0

我正在尝试使用 glmulti 和 glmer 进行模型平均并获得模型平均预测。我遵循了 glmulti 文档中的示例（'将 glmulti 与任何类型的统计模型一起使用，带有示例'，包含在包中）和本网站上提供的更新（glmulti 和 liner 混合模型）和包维护者的博客（https： //vcalcagnoresearch.wordpress.com/package-glmulti/）。我设法为 glmer 函数创建了一个包装器：

glmer.glmulti <- function (formula, data, family, random = "") {
glmer(paste(deparse(formula), random), data = data, family = binomial)
}

我已经设法为 glmer 分配了一个 getfit 方法，这样我就可以获得模型的平均系数：

setMethod('getfit', 'merMod', function(object, ...) {
summ=summary(object)$coef
summ1=summ[,1:2]
if (length(dimnames(summ)[[1]])==1) {
summ1=matrix(summ1, nr=1, dimnames=list(c("(Intercept)"),c("Estimate","Std. Error")))
}
cbind(summ1, df=rep(10000,length(fixef(object))))
})

下一步是为 glmer 分配一个预测函数。这是包文档中提供的示例：

predict.mer=function(objectmer,random=random, newdata, withRandom=F,se.fit=F, ...){
if (missing(newdata) || is.null(newdata)) {
DesignMat <- model.matrix(objectmer) }
else {
DesignMat=model.matrix(delete.response(terms(objectmer)),newdata)
}
output=DesignMat %*% fixef(objectmer)
if(withRandom){
z=unlist(ranef(objectmer))
if (missing(newdata) || is.null(newdata)) {
Zt<- objectmer@Zt
} else {
Zt<-as(as.factor(newdata[,names(ranef(objectmer))]),"sparseMatrix")
}
output = as.matrix(output + t(Zt) %*% z)
}
if(se.fit){
pvar <- diag(DesignMat %*% tcrossprod(vcov(objectmer),DesignMat))
if(withRandom){
pvar <- pvar+ VarCorr(objectmer)[[1]]
}
output=list(fit=output,se.fit=sqrt(pvar))
}
return(output)
}

然后得到模型平均预测（bab在示例中是一个拟合的 glmulti 对象）：

> predict(bab, se.fit=T, withR=T)
Error in predict.merMod(coffee[[i]], se.fit = se.fit, ...) :
cannot calculate predictions with both standard errors and random effects

我也试过：

> predict(bab, se.fit=T, withR=F)
Error in predict.merMod(coffee[[i]], se.fit = se.fit, ...) :
cannot calculate predictions with both standard errors and random effects

和：

> predict(bab, se.fit=F, withR=T)
Error in waou %*% t(matrix(unlist(preds), nrow = nbpo)) :
non-conformable arguments
In addition: There were 50 or more warnings (use warnings() to see the first 50)

我不太确定出了什么问题，尽管它可能很明显。我收集到 lme4 包自编写此示例以来已更新和更改，因此可能与此有关（？）。

另一种可能性是文档说这个函数只会处理一个随机变量。我的模型有两个嵌套的随机变量：x ~ y + z + w + (1|u/v).

我需要 a) 让它工作，b) 更新函数以处理两个随机变量。任何建议将不胜感激。

4

0 回答 0