我希望为一组线性模型制作一个汇总表。模型是 lmer()、glmer() 或 gamlss()。我正在尝试将其中 6 个模型的结果打印到一个表中。但是,当我尝试使用 sjPlot::tab_model 执行此操作时,我收到以下错误消息,“错误:$ 运算符对原子向量无效”。
我删除了两个 gamlss() 模型以查看表格是否会打印并且确实如此。在我看来,问题在于 sjPlot::tab_model 不能与 gamlss() 模型一起使用。
数据库
recruitment_data <- structure(list(Site_long = structure(c(1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 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, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("Hanauma Bay",
"Waikiki"), class = "factor"), Shelter = structure(c(1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 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, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 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, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("High",
"Low"), class = "factor"), `Module #` = structure(c(7L, 7L, 7L,
7L, 7L, 7L, 7L, 7L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 11L, 11L,
11L, 11L, 11L, 11L, 11L, 11L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L,
10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 12L, 12L, 12L, 12L, 12L,
12L, 12L, 12L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 6L, 6L, 6L, 6L,
6L, 6L, 6L, 6L), .Label = c("111", "112", "113", "114", "115",
"116", "211", "212", "213", "214", "215", "216"), class = "factor"),
TimeStep = c(4L, 4L, 5L, 5L, 6L, 6L, 7L, 7L, 4L, 4L, 5L,
5L, 6L, 6L, 7L, 7L, 4L, 4L, 5L, 5L, 6L, 6L, 7L, 7L, 4L, 4L,
5L, 5L, 6L, 6L, 7L, 7L, 4L, 4L, 5L, 5L, 6L, 6L, 7L, 7L, 4L,
4L, 5L, 5L, 6L, 6L, 7L, 7L, 4L, 4L, 5L, 5L, 6L, 6L, 7L, 7L,
4L, 4L, 5L, 5L, 6L, 6L, 7L, 7L, 4L, 4L, 5L, 5L, 6L, 6L, 7L,
7L, 4L, 4L, 5L, 5L, 6L, 6L, 7L, 7L, 4L, 4L, 5L, 5L, 6L, 6L,
7L, 7L, 4L, 4L, 5L, 5L, 6L, 6L, 7L, 7L), Side = 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), .Label = c("N", "S"), class = "factor"),
recruits = c(0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 2, 2, 1,
0, 0, 3, 0, 0, 1, 1, 0, 2, 0, 1, 0, 0, 0, 0, 0, 0, 0, 2,
3, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0,
0, 2, 2, 0, 1, 1, 0, 1, 3, 1, 0, 1, 1, 0, 1, 5, 2, 0, 1,
1, 2, 0, 3, 1, 2, 2, 3, 6, 5, 2, 0, 1, 2, 0, 4, 1, 4, 1,
0, 0, 4, 0, 1)), row.names = c(NA, -96L), class = c("grouped_df",
"tbl_df", "tbl", "data.frame"), vars = c("Site_long", "Shelter",
"Module #", "TimeStep"), drop = TRUE)
survival_data <- structure(list(Date = structure(c(17288, 17288, 17288, 17288,
17288, 17288, 17292, 17299, 17299, 17304, 17306, 17386, 17386,
17386, 17386, 17386, 17386, 17387, 17387, 17387, 17389, 17389,
17389, 17390, 17398, 17404, 17475, 17475, 17477, 17480, 17482,
17484, 17484, 17484, 17484, 17484, 17484, 17489, 17575, 17575,
17575, 17575, 17575, 17582, 17586, 17594, 17600, 17601, 17603,
17603), class = "Date"), Year = structure(c(1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("17",
"18"), class = "factor"), Site = structure(c(1L, 1L, 1L, 1L,
1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L,
1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("HAN",
"WAI"), class = "factor"), Treatment = c("CLO", "CLO", "CLO",
"OPE", "OPE", "OPE", "CLO", "CLO", "OPE", "OPE", "CLO", "CLO",
"CLO", "CLO", "OPE", "OPE", "OPE", "CLO", "OPE", "OPE", "CLO",
"CLO", "OPE", "OPE", "CLO", "CLO", "CLO", "OPE", "CLO", "OPE",
"CLO", "CLO", "CLO", "CLO", "OPE", "OPE", "OPE", "OPE", "CLO",
"CLO", "CLO", "OPE", "OPE", "OPE", "CLO", "OPE", "OPE", "CLO",
"CLO", "OPE"), `Module #` = c(212L, 214L, 216L, 211L, 213L, 215L,
116L, 114L, 115L, 113L, 112L, 212L, 214L, 216L, 211L, 213L, 215L,
116L, 111L, 115L, 112L, 114L, 115L, 113L, 114L, 114L, 112L, 115L,
116L, 113L, 114L, 212L, 214L, 216L, 211L, 213L, 215L, 111L, 212L,
214L, 216L, 213L, 215L, 211L, 116L, 115L, 113L, 114L, 112L, 111L
), n.x = c(1L, 1L, 2L, 2L, 3L, 3L, 6L, 7L, 5L, 4L, 2L, 2L, 2L,
1L, 2L, 3L, 5L, 10L, 1L, 4L, 10L, 13L, 6L, 5L, 2L, 2L, 8L, 6L,
10L, 8L, 12L, 2L, 6L, 2L, 2L, 3L, 5L, 2L, 2L, 5L, 1L, 8L, 9L,
2L, 10L, 15L, 10L, 16L, 12L, 4L), n.y = c(1, 1, 0, 2, 3, 3, 6,
7, 4, 4, 2, 2, 2, 1, 2, 3, 5, 9, 1, 3, 10, 11, 5, 5, 1, 1, 7,
6, 7, 8, 11, 2, 5, 2, 2, 3, 5, 2, 2, 5, 1, 7, 7, 0, 8, 14, 9,
9, 9, 4), `%_Survival` = c(100, 100, 0, 100, 100, 100, 100, 100,
80, 100, 100, 100, 100, 100, 100, 100, 100, 90, 100, 75, 100,
84.6153846153846, 83.3333333333333, 100, 50, 50, 87.5, 100, 70,
100, 91.6666666666667, 100, 83.3333333333333, 100, 100, 100,
100, 100, 100, 100, 100, 87.5, 77.7777777777778, 0, 80, 93.3333333333333,
90, 56.25, 75, 100), Quarter = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3,
3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4), Mortality = c(0,
0, 2, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 2,
1, 0, 1, 1, 1, 0, 3, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2,
2, 2, 1, 1, 7, 3, 0), Survival = c(100, 100, 0, 100, 100, 100,
100, 100, 80, 100, 100, 100, 100, 100, 100, 100, 100, 90, 100,
75, 100, 84.6153846153846, 83.3333333333333, 100, 50, 50, 87.5,
100, 70, 100, 91.6666666666667, 100, 83.3333333333333, 100, 100,
100, 100, 100, 100, 100, 100, 87.5, 77.7777777777778, 0, 80,
93.3333333333333, 90, 56.25, 75, 100), Module = c(212L, 214L,
216L, 211L, 213L, 215L, 116L, 114L, 115L, 113L, 112L, 212L, 214L,
216L, 211L, 213L, 215L, 116L, 111L, 115L, 112L, 114L, 115L, 113L,
114L, 114L, 112L, 115L, 116L, 113L, 114L, 212L, 214L, 216L, 211L,
213L, 215L, 111L, 212L, 214L, 216L, 213L, 215L, 211L, 116L, 115L,
113L, 114L, 112L, 111L), Survival_prop = c(1, 1, 0, 1, 1, 1,
1, 1, 0.8, 1, 1, 1, 1, 1, 1, 1, 1, 0.9, 1, 0.75, 1, 0.846153846153846,
0.833333333333333, 1, 0.5, 0.5, 0.875, 1, 0.7, 1, 0.916666666666667,
1, 0.833333333333333, 1, 1, 1, 1, 1, 1, 1, 1, 0.875, 0.777777777777778,
0, 0.8, 0.933333333333333, 0.9, 0.5625, 0.75, 1), Date_new = structure(c(17378,
17378, 17378, 17378, 17378, 17378, 17382, 17389, 17389, 17394,
17396, 17476, 17476, 17476, 17476, 17476, 17476, 17477, 17477,
17477, 17479, 17479, 17479, 17480, 17488, 17494, 17565, 17565,
17567, 17570, 17572, 17574, 17574, 17574, 17574, 17574, 17574,
17579, 17665, 17665, 17665, 17665, 17665, 17672, 17676, 17684,
17690, 17691, 17693, 17693), class = "Date"), Site_long = structure(c(1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L,
1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,
2L), .Label = c("Hanauma Bay", "Waikiki"), class = "factor"),
Treatment_long = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 1L,
1L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 1L, 1L,
2L, 2L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L,
2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 1L, 1L, 2L), .Label = c("Closed",
"Open"), class = "factor"), Shelter = c("Low", "Low", "Low",
"High", "High", "High", "Low", "Low", "High", "High", "Low",
"Low", "Low", "Low", "High", "High", "High", "Low", "High",
"High", "Low", "Low", "High", "High", "Low", "Low", "Low",
"High", "Low", "High", "Low", "Low", "Low", "Low", "High",
"High", "High", "High", "Low", "Low", "Low", "High", "High",
"High", "Low", "High", "High", "Low", "Low", "High")), row.names = c(NA,
-50L), vars = c("Date", "Year", "Site", "Treatment", "Module #"
), labels = structure(list(Date = structure(c(17288, 17288, 17288,
17288, 17288, 17288, 17292, 17299, 17299, 17304, 17306, 17386,
17386, 17386, 17386, 17386, 17386, 17387, 17387, 17387, 17389,
17389, 17389, 17390, 17398, 17404, 17475, 17475, 17477, 17480,
17482, 17484, 17484, 17484, 17484, 17484, 17484, 17489, 17575,
17575, 17575, 17575, 17575, 17582, 17586, 17594, 17600, 17601,
17603, 17603), class = "Date"), Year = c(17, 17, 17, 17, 17,
17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17,
17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17,
17, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18), Site = structure(c(1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L,
1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,
2L), .Label = c("HAN", "WAI"), class = "factor"), Treatment = c("CLO",
"CLO", "CLO", "OPE", "OPE", "OPE", "CLO", "CLO", "OPE", "OPE",
"CLO", "CLO", "CLO", "CLO", "OPE", "OPE", "OPE", "CLO", "OPE",
"OPE", "CLO", "CLO", "OPE", "OPE", "CLO", "CLO", "CLO", "OPE",
"CLO", "OPE", "CLO", "CLO", "CLO", "CLO", "OPE", "OPE", "OPE",
"OPE", "CLO", "CLO", "CLO", "OPE", "OPE", "OPE", "CLO", "OPE",
"OPE", "CLO", "CLO", "OPE"), `Module #` = c(212L, 214L, 216L,
211L, 213L, 215L, 116L, 114L, 115L, 113L, 112L, 212L, 214L, 216L,
211L, 213L, 215L, 116L, 111L, 115L, 112L, 114L, 115L, 113L, 114L,
114L, 112L, 115L, 116L, 113L, 114L, 212L, 214L, 216L, 211L, 213L,
215L, 111L, 212L, 214L, 216L, 213L, 215L, 211L, 116L, 115L, 113L,
114L, 112L, 111L)), row.names = c(NA, -50L), class = "data.frame", vars = c("Date",
"Year", "Site", "Treatment", "Module #"), drop = TRUE), indices = list(
0L, 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), drop = TRUE, group_sizes = c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L), biggest_group_size = 1L, class = c("grouped_df", "tbl_df",
"tbl", "data.frame"))
分析
library(lme4)
library(gamlss)
recruitment_glmer_n3 <- glmer(recruits ~ Site_long*Shelter + (1|module_recruit), data = n3, family = poisson, na.action = "na.fail")
summary(recruitment_glmer_n3)
survival_gamlss <- gamlss(Survival_prop ~ Site_long*Treatment_long + (1|module_survival), data = survival_results_long_2, family = BEINF())
summary(survival_gamlss)
表代码
library(sjPlot)
tab_model(recruitment_glmer_n3, survival_gamlss)
有没有办法使用 sjPlot::tab_model 来获取带有 gamlss 模型对象的 html 表输出,或者是否有其他软件包可以推荐用于为 gamlss 以及 lmer 和 glmer 的线性模型对象制作出版质量表?谢谢!