使用扫帚从 glmmTMB 零 inflation 模型获取系数时出现错误消息
Error message when using broom to get coefficients from glmmTMB zero inflation models
使用 Owls 数据和 glmmTMB
包,我想直观地比较不同零膨胀模型的回归系数,这些模型在所使用的系列(ZIPOISS、ZINB1、ZINB2)和 with/out 偏移量 (logBroodSize)。
但是我的第一个问题是获取系数。包 broom
中的 tidy
函数应该为您提供稍后用 ggplot
绘制它们的系数,但是当我尝试获取它们时出现以下错误:
modList= list(zipoiss, zinb1, zinb2, zinb1_bs, zinb2_bs)
coefs= ldply(modList,tidy,effect="fixed",conf.int=TRUE,
.id="model") %>%
mutate(term=abbfun(term)) %>%
select(model,term,estimate,conf.low,conf.high) %>%
filter(!term %in% c("Intercept","Intercept.1","NCalls","zi_NCalls"))
Error in as.data.frame.default(x) :
cannot coerce class ""glmmTMB"" to a data.frame
Also: Warning message:
In tidy.default(X[[i]], ...) :
No method for tidying an S3 object of class glmmTMB , using as.data.frame
知道哪里出了问题吗?我已经被告知没有正确版本的 broom
可能是原因,但是我已经安装了正确版本的它......接下来提供了可重现示例的代码:
# Packages and dataset
library(glmmTMB)
library(broom) # devtools::install_github("bbolker/broom")
library(plyr)
library(dplyr)
data(Owls,package="glmmTMB")
Owls = plyr::rename(Owls, c(SiblingNegotiation="NCalls"))
Owls = transform(Owls,
ArrivalTime=scale(ArrivalTime, center=TRUE, scale=FALSE),
obs=factor(seq(nrow(Owls))))
# Models
zipoiss<-glmmTMB(NCalls~(FoodTreatment+ArrivalTime)*SexParent+
offset(logBroodSize)+(1|Nest),
data=Owls,
ziformula = ~ 1,
family="poisson")
zinb2<- glmmTMB(NCalls~(FoodTreatment+ArrivalTime)*SexParent+
offset(logBroodSize)+(1|Nest),
data=Owls,
ziformula = ~ 1,
family="nbinom2")
zinb1 <- glmmTMB(NCalls~(FoodTreatment+ArrivalTime)*SexParent+
offset(logBroodSize)+(1|Nest),
data=Owls,
ziformula = ~ 1,
family="nbinom1")
zinb1_bs<- glmmTMB(NCalls~(FoodTreatment+ArrivalTime)*SexParent+
BroodSize+(1|Nest),
data=Owls,
ziformula = ~ 1,
family="nbinom1")
zinb2_bs<- glmmTMB(NCalls~(FoodTreatment+ArrivalTime)*SexParent+
BroodSize+(1|Nest),
data=Owls,
ziformula = ~ 1,family="nbinom2")
# Get coefficients ("coefs" does not work yet...)
modList = list(zipoiss, zinb1, zinb2, zinb1_bs, zinb2_bs)
coefs = ldply(modList,tidy,effect="fixed",conf.int=TRUE,
.id="model") %>%
mutate(term=abbfun(term)) %>%
select(model,term,estimate,conf.low,conf.high) %>%
filter(!term %in% c("Intercept","Intercept.1","NCalls","zi_NCalls"))
这现在(截至今天)适用于 new/under-development broom.mixed 包,例如
devtools::install_github("bbolker/broom.mixed")
希望这会很快出现在 CRAN 上,但它只是中等优先级
对我来说,我不想在至少 90% 烘烤之前将其发布到 CRAN。欢迎请求请求!
最后一步稍微改了一下(一方面,我好像没有abbfun
):
modList = lme4:::namedList(zipoiss, zinb1, zinb2, zinb1_bs, zinb2_bs)
coefs = ldply(modList,tidy,effect="fixed",conf.int=TRUE,
.id="model") %>%
select(model,term,component,estimate,conf.low,conf.high) %>%
filter(!term %in% c("(Intercept)","NCalls"))
警告:为glmmTMB
模型开发这些工具非常棒new/experimental;你应该 (1) sanity-check 仔细检查你的结果,并且 (2) report an issue 如果某些事情没有按预期工作。
使用 Owls 数据和 glmmTMB
包,我想直观地比较不同零膨胀模型的回归系数,这些模型在所使用的系列(ZIPOISS、ZINB1、ZINB2)和 with/out 偏移量 (logBroodSize)。
但是我的第一个问题是获取系数。包 broom
中的 tidy
函数应该为您提供稍后用 ggplot
绘制它们的系数,但是当我尝试获取它们时出现以下错误:
modList= list(zipoiss, zinb1, zinb2, zinb1_bs, zinb2_bs)
coefs= ldply(modList,tidy,effect="fixed",conf.int=TRUE,
.id="model") %>%
mutate(term=abbfun(term)) %>%
select(model,term,estimate,conf.low,conf.high) %>%
filter(!term %in% c("Intercept","Intercept.1","NCalls","zi_NCalls"))
Error in as.data.frame.default(x) :
cannot coerce class ""glmmTMB"" to a data.frame
Also: Warning message:
In tidy.default(X[[i]], ...) :
No method for tidying an S3 object of class glmmTMB , using as.data.frame
知道哪里出了问题吗?我已经被告知没有正确版本的 broom
可能是原因,但是我已经安装了正确版本的它......接下来提供了可重现示例的代码:
# Packages and dataset
library(glmmTMB)
library(broom) # devtools::install_github("bbolker/broom")
library(plyr)
library(dplyr)
data(Owls,package="glmmTMB")
Owls = plyr::rename(Owls, c(SiblingNegotiation="NCalls"))
Owls = transform(Owls,
ArrivalTime=scale(ArrivalTime, center=TRUE, scale=FALSE),
obs=factor(seq(nrow(Owls))))
# Models
zipoiss<-glmmTMB(NCalls~(FoodTreatment+ArrivalTime)*SexParent+
offset(logBroodSize)+(1|Nest),
data=Owls,
ziformula = ~ 1,
family="poisson")
zinb2<- glmmTMB(NCalls~(FoodTreatment+ArrivalTime)*SexParent+
offset(logBroodSize)+(1|Nest),
data=Owls,
ziformula = ~ 1,
family="nbinom2")
zinb1 <- glmmTMB(NCalls~(FoodTreatment+ArrivalTime)*SexParent+
offset(logBroodSize)+(1|Nest),
data=Owls,
ziformula = ~ 1,
family="nbinom1")
zinb1_bs<- glmmTMB(NCalls~(FoodTreatment+ArrivalTime)*SexParent+
BroodSize+(1|Nest),
data=Owls,
ziformula = ~ 1,
family="nbinom1")
zinb2_bs<- glmmTMB(NCalls~(FoodTreatment+ArrivalTime)*SexParent+
BroodSize+(1|Nest),
data=Owls,
ziformula = ~ 1,family="nbinom2")
# Get coefficients ("coefs" does not work yet...)
modList = list(zipoiss, zinb1, zinb2, zinb1_bs, zinb2_bs)
coefs = ldply(modList,tidy,effect="fixed",conf.int=TRUE,
.id="model") %>%
mutate(term=abbfun(term)) %>%
select(model,term,estimate,conf.low,conf.high) %>%
filter(!term %in% c("Intercept","Intercept.1","NCalls","zi_NCalls"))
这现在(截至今天)适用于 new/under-development broom.mixed 包,例如
devtools::install_github("bbolker/broom.mixed")
希望这会很快出现在 CRAN 上,但它只是中等优先级 对我来说,我不想在至少 90% 烘烤之前将其发布到 CRAN。欢迎请求请求!
最后一步稍微改了一下(一方面,我好像没有abbfun
):
modList = lme4:::namedList(zipoiss, zinb1, zinb2, zinb1_bs, zinb2_bs)
coefs = ldply(modList,tidy,effect="fixed",conf.int=TRUE,
.id="model") %>%
select(model,term,component,estimate,conf.low,conf.high) %>%
filter(!term %in% c("(Intercept)","NCalls"))
警告:为glmmTMB
模型开发这些工具非常棒new/experimental;你应该 (1) sanity-check 仔细检查你的结果,并且 (2) report an issue 如果某些事情没有按预期工作。