R как заставить функцию итога работать на моделях пакетов censReg - PullRequest
Новая
       34

R как заставить функцию итога работать на моделях пакетов censReg

0 голосов
/ 14 марта 2020

Я пытаюсь использовать функцию censReg в R, которая является частью пакета censReg. Я пытаюсь смоделировать данные о биомассе, которые являются непрерывными по природе и содержат нули. Ноль - это наименьшее возможное значение, потому что нет такой вещи, как отрицательная биомасса. Я решил использовать пакет censReg как способ справиться с непрерывным распределением переменной total_biomass с сильно завышенным нулем. Мне удалось запустить приведенную ниже модель, однако, когда я пытаюсь запустить функцию summary () на модели, я получаю следующую ошибку: 

Error in printCoefmat(coef(x, logSigma = logSigma), digits = digits) : 
  'x' must be coefficient matrix/data frame

Глядя на это сообщение об ошибке, я не понимаю, что это значит, ни в чем проблема с моим кодом или базой данных. Может ли кто-нибудь предоставить мне какой-либо дополнительный код или настройки, которые мне нужно сделать, чтобы успешно получить сводку модели?

База данных 

library(tidyverse)
library(censReg)

mean_fish_totals <- structure(list(Date = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 2L, 
2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 
4L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 
7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 9L, 10L, 
10L, 10L, 10L, 10L, 10L, 11L, 11L, 11L, 11L, 11L, 11L, 12L, 12L, 
12L, 12L, 12L, 12L, 13L, 13L, 13L, 13L, 13L, 13L, 14L, 14L, 14L, 
14L, 14L, 14L, 15L, 15L, 15L, 15L, 15L, 15L, 16L, 16L, 16L, 16L, 
16L, 16L, 17L, 17L, 17L, 17L, 17L, 17L, 18L, 18L, 18L, 18L, 18L, 
18L, 19L, 19L, 19L, 19L, 19L, 19L, 20L, 20L, 20L, 20L, 20L, 20L, 
21L, 21L, 21L, 21L, 22L, 22L, 22L, 23L, 23L, 23L, 24L, 24L, 25L, 
25L, 25L, 25L, 25L, 25L, 26L, 26L, 26L, 26L, 26L, 26L), .Label = c("11/28/17", 
"12/10/19", "12/14/19", "12/3/17", "12/4/18", "12/8/18", "2/25/17", 
"3/19/19", "3/22/19", "4/18/17", "5/15/18", "5/20/18", "5/25/17", 
"6/3/19", "6/4/19", "6/6/17", "8/28/18", "9/1/18", "9/10/19", 
"9/15/19", "9/20/16", "9/22/16", "9/25/16", "9/27/16", "9/5/17", 
"9/7/17"), class = "factor"), `Module #` = c(211L, 212L, 213L, 
214L, 215L, 216L, 211L, 212L, 213L, 214L, 215L, 216L, 111L, 112L, 
113L, 114L, 115L, 116L, 111L, 112L, 113L, 114L, 115L, 116L, 211L, 
212L, 213L, 214L, 215L, 216L, 111L, 112L, 113L, 114L, 115L, 116L, 
111L, 112L, 113L, 114L, 115L, 116L, 211L, 212L, 213L, 214L, 215L, 
216L, 111L, 112L, 113L, 114L, 115L, 116L, 211L, 212L, 213L, 214L, 
215L, 216L, 211L, 212L, 213L, 214L, 215L, 216L, 111L, 112L, 113L, 
114L, 115L, 116L, 111L, 112L, 113L, 114L, 115L, 116L, 111L, 112L, 
113L, 114L, 115L, 116L, 211L, 212L, 213L, 214L, 215L, 216L, 211L, 
212L, 213L, 214L, 215L, 216L, 211L, 212L, 213L, 214L, 215L, 216L, 
111L, 112L, 113L, 114L, 115L, 116L, 211L, 212L, 213L, 214L, 215L, 
216L, 111L, 112L, 113L, 114L, 115L, 116L, 213L, 214L, 215L, 216L, 
114L, 115L, 116L, 111L, 112L, 113L, 211L, 212L, 211L, 212L, 213L, 
214L, 215L, 216L, 111L, 112L, 113L, 114L, 115L, 116L), Site_long = structure(c(1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 
2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 
1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 
2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("Hanauma Bay", 
"Waikiki"), class = "factor"), Treatment_long = structure(c(2L, 
1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 
1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 
1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 
1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 
1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 
1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 
1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 
1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 
2L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L), .Label = c("Closed", 
"Open"), class = "factor"), Shelter = structure(c(1L, 2L, 1L, 
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 2L, 1L, 1L, 
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L), .Label = c("High", 
"Low"), class = "factor"), TimeStep = c(5, 5, 5, 5, 5, 5, 15, 
15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 5, 5, 5, 5, 5, 5, 
11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 0.75, 0.75, 0.75, 
0.75, 0.75, 0.75, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 
12, 0.75, 2, 0.75, 0.75, 0.75, 2, 7, 7, 7, 7, 7, 7, 7, 7, 7, 
7, 7, 7, 1, 1, 1, 1, 1, 1, 13, 13, 13, 13, 13, 13, 13, 13, 13, 
13, 13, 13, 1, 1, 1, 1, 1, 1, 10, 10, 10, 10, 10, 10, 10, 10, 
10, 10, 10, 10, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 
0.5, 0.5, 0.5, 1, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 1, 4, 4, 
4, 4, 4, 4, 4, 4, 4, 4, 4, 4), total_biomass = c(0.0526394784788566, 
0.00650991088549517, 0.180596698411345, 0, 0.526717015131238, 
0, 0.206894204748519, 0.0164498212264971, 0, 0, 0, 0, 0.201603325166693, 
0, 0.283172330030785, 0, 0, 0, 0.209199851062596, 0, 0.283172330030785, 
0, 0.394899281159044, 0.176129136061979, 0.169586003898729, 0, 
0.120034320302602, 0, 0.996727938559748, 0, 0, 0, 0, 0, 0.257380940140258, 
0.443609909402316, 0.308392987176445, 0, 0.748305018557033, 0, 
0.169586003898729, 0, 0.120034320302602, 0, 0.120034320302602, 
0, 1.08474493759439, 0, 0, 0.0745413774887557, 0, 0.0467403151010407, 
0.233352036920048, 0, 0.0899664423257818, 0, 0.308392987176445, 
0, 0, 0, 10.6461511880577, 0, 26.4504921170652, 0, 0.526717015131238, 
0, 0, 0, 0.403061371653634, 0.209695276260751, 0.120034320302602, 
0.0206199419933497, 0.078489026854395, 0, 0.165344302422082, 
0, 0.0317487117533543, 0, 0, 0, 0, 0, 0.94950027744447, 0, 0, 
0, 0, 0, 1.03519325399826, 0, 0.169586003898729, 0.125258604363503, 
0.310810458426215, 0, 10.6461511880577, 0, 0, 0, 0.976669817356845, 
0, 0.996727938559748, 0, 0, 0.0202327772014947, 0.0403651893214743, 
0, 0.168776380464422, 0, 0.133454408606973, 0, 0, 0.621724957549784, 
0.0164498212264971, 0, 0.110237738607749, 0, 0.116136901985565, 
0, 0, 0.0135959326713389, 0.00889824575015321, 0.078489026854395, 
0.16627064788403, 0, 0.028053154008064, 0.0526394784788566, 0.0419621766803908, 
0, 0, 0, 0, 0, 0.308392987176445, 0.0520989800170256, 0.222619650542138, 
0, 21.2935111117403, 15.7227719241434, 0.232861857335555, 0, 
0.0634974235067086, 0, 0.0492074004365164, 0), new_date = structure(c(17498, 
17498, 17498, 17498, 17498, 17498, 18240, 18240, 18240, 18240, 
18240, 18240, 18244, 18244, 18244, 18244, 18244, 18244, 17503, 
17503, 17503, 17503, 17503, 17503, 17869, 17869, 17869, 17869, 
17869, 17869, 17873, 17873, 17873, 17873, 17873, 17873, 17222, 
17222, 17222, 17222, 17222, 17222, 17974, 17974, 17974, 17974, 
17974, 17974, 17977, 17977, 17977, 17977, 17977, 17977, 17274, 
17274, 17274, 17274, 17274, 17274, 17666, 17666, 17666, 17666, 
17666, 17666, 17671, 17671, 17671, 17671, 17671, 17671, 17311, 
17311, 17311, 17311, 17311, 17311, 18050, 18050, 18050, 18050, 
18050, 18050, 18051, 18051, 18051, 18051, 18051, 18051, 17323, 
17323, 17323, 17323, 17323, 17323, 17771, 17771, 17771, 17771, 
17771, 17771, 17775, 17775, 17775, 17775, 17775, 17775, 18149, 
18149, 18149, 18149, 18149, 18149, 18154, 18154, 18154, 18154, 
18154, 18154, 17064, 17064, 17064, 17064, 17066, 17066, 17066, 
17069, 17069, 17069, 17071, 17071, 17414, 17414, 17414, 17414, 
17414, 17414, 17416, 17416, 17416, 17416, 17416, 17416), class = "Date")), row.names = c(NA, 
-144L), groups = structure(list(Date = structure(c(1L, 1L, 1L, 
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 
4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 
6L, 7L, 7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 
9L, 9L, 9L, 10L, 10L, 10L, 10L, 10L, 10L, 11L, 11L, 11L, 11L, 
11L, 11L, 12L, 12L, 12L, 12L, 12L, 12L, 13L, 13L, 13L, 13L, 13L, 
13L, 14L, 14L, 14L, 14L, 14L, 14L, 15L, 15L, 15L, 15L, 15L, 15L, 
16L, 16L, 16L, 16L, 16L, 16L, 17L, 17L, 17L, 17L, 17L, 17L, 18L, 
18L, 18L, 18L, 18L, 18L, 19L, 19L, 19L, 19L, 19L, 19L, 20L, 20L, 
20L, 20L, 20L, 20L, 21L, 21L, 21L, 21L, 22L, 22L, 22L, 23L, 23L, 
23L, 24L, 24L, 25L, 25L, 25L, 25L, 25L, 25L, 26L, 26L, 26L, 26L, 
26L, 26L), .Label = c("11/28/17", "12/10/19", "12/14/19", "12/3/17", 
"12/4/18", "12/8/18", "2/25/17", "3/19/19", "3/22/19", "4/18/17", 
"5/15/18", "5/20/18", "5/25/17", "6/3/19", "6/4/19", "6/6/17", 
"8/28/18", "9/1/18", "9/10/19", "9/15/19", "9/20/16", "9/22/16", 
"9/25/16", "9/27/16", "9/5/17", "9/7/17"), class = "factor"), 
    `Module #` = c(211L, 212L, 213L, 214L, 215L, 216L, 211L, 
    212L, 213L, 214L, 215L, 216L, 111L, 112L, 113L, 114L, 115L, 
    116L, 111L, 112L, 113L, 114L, 115L, 116L, 211L, 212L, 213L, 
    214L, 215L, 216L, 111L, 112L, 113L, 114L, 115L, 116L, 111L, 
    112L, 113L, 114L, 115L, 116L, 211L, 212L, 213L, 214L, 215L, 
    216L, 111L, 112L, 113L, 114L, 115L, 116L, 211L, 212L, 213L, 
    214L, 215L, 216L, 211L, 212L, 213L, 214L, 215L, 216L, 111L, 
    112L, 113L, 114L, 115L, 116L, 111L, 112L, 113L, 114L, 115L, 
    116L, 111L, 112L, 113L, 114L, 115L, 116L, 211L, 212L, 213L, 
    214L, 215L, 216L, 211L, 212L, 213L, 214L, 215L, 216L, 211L, 
    212L, 213L, 214L, 215L, 216L, 111L, 112L, 113L, 114L, 115L, 
    116L, 211L, 212L, 213L, 214L, 215L, 216L, 111L, 112L, 113L, 
    114L, 115L, 116L, 213L, 214L, 215L, 216L, 114L, 115L, 116L, 
    111L, 112L, 113L, 211L, 212L, 211L, 212L, 213L, 214L, 215L, 
    216L, 111L, 112L, 113L, 114L, 115L, 116L), Site_long = c("Hanauma Bay", 
    "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", 
    "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", 
    "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", "Waikiki", "Waikiki", 
    "Waikiki", "Waikiki", "Waikiki", "Waikiki", "Waikiki", "Waikiki", 
    "Waikiki", "Waikiki", "Waikiki", "Waikiki", "Hanauma Bay", 
    "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", 
    "Hanauma Bay", "Waikiki", "Waikiki", "Waikiki", "Waikiki", 
    "Waikiki", "Waikiki", "Waikiki", "Waikiki", "Waikiki", "Waikiki", 
    "Waikiki", "Waikiki", "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", 
    "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", "Waikiki", "Waikiki", 
    "Waikiki", "Waikiki", "Waikiki", "Waikiki", "Hanauma Bay", 
    "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", 
    "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", 
    "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", "Waikiki", "Waikiki", 
    "Waikiki", "Waikiki", "Waikiki", "Waikiki", "Waikiki", "Waikiki", 
    "Waikiki", "Waikiki", "Waikiki", "Waikiki", "Waikiki", "Waikiki", 
    "Waikiki", "Waikiki", "Waikiki", "Waikiki", "Hanauma Bay", 
    "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", 
    "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", 
    "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", 
    "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", 
    "Hanauma Bay", "Waikiki", "Waikiki", "Waikiki", "Waikiki", 
    "Waikiki", "Waikiki", "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", 
    "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", "Waikiki", "Waikiki", 
    "Waikiki", "Waikiki", "Waikiki", "Waikiki", "Hanauma Bay", 
    "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", "Waikiki", "Waikiki", 
    "Waikiki", "Waikiki", "Waikiki", "Waikiki", "Hanauma Bay", 
    "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", 
    "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", "Waikiki", "Waikiki", 
    "Waikiki", "Waikiki", "Waikiki", "Waikiki"), Treatment_long = c("Open", 
    "Closed", "Open", "Closed", "Open", "Closed", "Open", "Closed", 
    "Open", "Closed", "Open", "Closed", "Open", "Closed", "Open", 
    "Closed", "Open", "Closed", "Open", "Closed", "Open", "Closed", 
    "Open", "Closed", "Open", "Closed", "Open", "Closed", "Open", 
    "Closed", "Open", "Closed", "Open", "Closed", "Open", "Closed", 
    "Open", "Closed", "Open", "Closed", "Open", "Closed", "Open", 
    "Closed", "Open", "Closed", "Open", "Closed", "Open", "Closed", 
    "Open", "Closed", "Open", "Closed", "Open", "Closed", "Open", 
    "Closed", "Open", "Closed", "Open", "Closed", "Open", "Closed", 
    "Open", "Closed", "Open", "Closed", "Open", "Closed", "Open", 
    "Closed", "Open", "Closed", "Open", "Closed", "Open", "Closed", 
    "Open", "Closed", "Open", "Closed", "Open", "Closed", "Open", 
    "Closed", "Open", "Closed", "Open", "Closed", "Open", "Closed", 
    "Open", "Closed", "Open", "Closed", "Open", "Closed", "Open", 
    "Closed", "Open", "Closed", "Open", "Closed", "Open", "Closed", 
    "Open", "Closed", "Open", "Closed", "Open", "Closed", "Open", 
    "Closed", "Open", "Closed", "Open", "Closed", "Open", "Closed", 
    "Open", "Closed", "Open", "Closed", "Closed", "Open", "Closed", 
    "Open", "Closed", "Open", "Open", "Closed", "Open", "Closed", 
    "Open", "Closed", "Open", "Closed", "Open", "Closed", "Open", 
    "Closed", "Open", "Closed"), Shelter = c("High", "Low", "High", 
    "Low", "High", "Low", "High", "Low", "High", "Low", "High", 
    "Low", "High", "Low", "High", "Low", "High", "Low", "High", 
    "Low", "High", "Low", "High", "Low", "High", "Low", "High", 
    "Low", "High", "Low", "High", "Low", "High", "Low", "High", 
    "Low", "High", "Low", "High", "Low", "High", "Low", "High", 
    "Low", "High", "Low", "High", "Low", "High", "Low", "High", 
    "Low", "High", "Low", "High", "Low", "High", "Low", "High", 
    "Low", "High", "Low", "High", "Low", "High", "Low", "High", 
    "Low", "High", "Low", "High", "Low", "High", "Low", "High", 
    "Low", "High", "Low", "High", "Low", "High", "Low", "High", 
    "Low", "High", "Low", "High", "Low", "High", "Low", "High", 
    "Low", "High", "Low", "High", "Low", "High", "Low", "High", 
    "Low", "High", "Low", "High", "Low", "High", "Low", "High", 
    "Low", "High", "Low", "High", "Low", "High", "Low", "High", 
    "Low", "High", "Low", "High", "Low", "High", "Low", "High", 
    "Low", "Low", "High", "Low", "High", "Low", "High", "High", 
    "Low", "High", "Low", "High", "Low", "High", "Low", "High", 
    "Low", "High", "Low", "High", "Low"), .rows = list(1L, 2L, 
        3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 
        15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L, 25L, 
        26L, 27L, 28L, 29L, 30L, 31L, 32L, 33L, 34L, 35L, 36L, 
        37L, 38L, 39L, 40L, 41L, 42L, 43L, 44L, 45L, 46L, 47L, 
        48L, 49L, 50L, 51L, 52L, 53L, 54L, 55L, 56L, 57L, 58L, 
        59L, 60L, 61L, 62L, 63L, 64L, 65L, 66L, 67L, 68L, 69L, 
        70L, 71L, 72L, 73L, 74L, 75L, 76L, 77L, 78L, 79L, 80L, 
        81L, 82L, 83L, 84L, 85L, 86L, 87L, 88L, 89L, 90L, 91L, 
        92L, 93L, 94L, 95L, 96L, 97L, 98L, 99L, 100L, 101L, 102L, 
        103L, 104L, 105L, 106L, 107L, 108L, 109L, 110L, 111L, 
        112L, 113L, 114L, 115L, 116L, 117L, 118L, 119L, 120L, 
        121L, 122L, 123L, 124L, 125L, 126L, 127L, 128L, 129L, 
        130L, 131L, 132L, 133L, 134L, 135L, 136L, 137L, 138L, 
        139L, 140L, 141L, 142L, 143L, 144L)), row.names = c(NA, 
-144L), class = c("tbl_df", "tbl", "data.frame"), .drop = TRUE), class = c("grouped_df", 
"tbl_df", "tbl", "data.frame"))

Анализ 

## Herbivorous fish with interaction
# variables
module_fish <- mean_fish_totals$`Module #`

# Data distribution #
plotNormalHistogram(mean_fish_totals$total_biomass)

## Generalized Mixed Effects Model Fishes ##
# glmer new database #
mean_fish_totals$Shelter <- factor(mean_fish_totals$Shelter, levels = c("Low", "High"), ordered = TRUE)
mean_fish_totals$Site_long <- as.factor(mean_fish_totals$Site_long)

# Censreg Model #
fish_mixed_effects_censreg <- censReg(total_biomass ~ Site_long*Shelter + (1|module_fish), left = 0, right = Inf, data = mean_fish_totals)
summary(fish_mixed_effects_censreg)

Заранее благодарен за ваш вклад!

1 Ответ

2 голосов
/ 15 марта 2020

Проблема в части "(1 | module_fi sh)", потому что censReg () не "знает" этот способ задания формулы регрессии. Если у вас есть данные панели, cenReg () может учитывать случайные эффекты в перехвате (см. «Виньетка» в пакете censReg), но не учитывает случайные эффекты в параметрах наклона. Вы можете оценить модель без случайных эффектов следующим образом: 

mean_fish_totals$module_fish <- mean_fish_totals$`Module #`

fish_censreg <- censReg(
  total_biomass ~ Site_long*Shelter + module_fish,
  left = 0, right = Inf, data = mean_fish_totals )

summary( fish_censreg )

Вы также можете оценить эту спецификацию с помощью функции tobit () пакета "AER": 

library( "AER" )
fish_tobit <- tobit(total_biomass ~ Site_long*Shelter + module_fish,
  data = mean_fish_totals)
summary( fish_tobit )

Это дает идентичные результаты : 

toMatrix <- function(x){ class( x ) <- "matrix"; x }
all.equal( toMatrix( coef( summary( fish_tobit ) ) ), 
  coef( summary( fish_censreg ) ), check.attributes = FALSE )

Остатки можно получить из моделей, оцененных с помощью tobit (): 

qqnorm( resid( fish_tobit ) ) 
qqline( resid( fish_tobit ) )

Обратите внимание, что существуют разные способы определения (и, таким образом, вычисления) остатков цензуры регрессионные модели (см., например, документацию residuals.survreg () в пакете «выживание», который используется внутри системы при получении остатков моделей, оцениваемых с помощью tobit ()).

Пожалуйста, не стесняйтесь свяжитесь со мной, если вы (или кто-либо еще) заинтересованы в реализации отсутствующих в настоящее время функций в пакете censReg: https://r-forge.r-project.org/projects/sampleselection/

...