R v 3.6.1 中的 lmer 错误和预测函数/从给定自变量的 lmer 模型中获取置信区间答案

【问题标题】：lmer errors and predict function in R v 3.6.1 / obtain confidence intervals from lmer model for a given independent variableR v 3.6.1 中的 lmer 错误和预测函数/从给定自变量的 lmer 模型中获取置信区间
【发布时间】：2019-11-27 03:27:34
【问题描述】：

这是两个问题合二为一，希望可以。

首先，我试图从lme4 包中的lmer 对象获取置信区间值。我以前使用过 R v 3.4.4，模型运行良好，我可以找到并显示图上均值拟合的置信区间。我最近升级到 R v 3.6.1，现在在使用 predict 函数时收到错误消息。我已将其作为代码的一部分在下面显示

我的主要问题是，如何计算给定值的置信上限和下限？对于传统的lm，我会使用：

new.dat <- data.frame(variable = ##)
predict(lm_object, newdata = new.dat, interval = 'confidence')

但这不适用于lmer 对象。

这是数据：

TYPE <- c(rep("A", 100), rep("B", 31), rep("C", 18))
MAX<-c(NA,32.6,19.5,23.5,0,17.3,0,31,35.3,23.9,20.8,18.3,10.6,19.4,0,9,14.5,
       27.1,0,27.5,21,0,14.7,23.7,17.4,13.7,30.7,25.3,NA,0,16.5,0,NA,18.5,23.9,
       8.6,11.9,21.5,0,20.3,10.1,0,20.2,33.6,40.6,21.9,16.6,18.3,0,28.3,36.4,0,
       29.4,25.7,24.8,25,0,36.9,19,19.3,27.8,20.4,19.2,0,25.5,26.3,30.6,0,27.8,
       5.7,0,21,19.7,15.3,0,16.5,14.5,17.2,31.7,13,21.5,20,32.5,0,6.8,26.2,0,
       24.6,21.2,0,32.3,17.3,29.2,43.1,26.2,0,29.5,26.1,36.8,10.9,45.76,17.41,
       62.475,0,11.82,57.12,0,41.35,52.935,13.01,0,56.095,60.345,56.645,78.775,
       69.565,47.98,15.28,16.46,12.91,0,14.76,29.185,29.26,0,77.72,78.25,0,45.875,
       0,40.27,16.43,27.065,45.44,71.38,21.875,0,33.625,45.825,51.79,39.705,27.46,
       36.61,44.21,62.38,0,120.295,26.61,0)
STEM_PERCENT<-c(-8.708496564,-8.708496564,-8.708496564,-8.708496564,0,-8.708496564,
                0,-8.708496564,-8.708496564,-8.708496564,-8.708496564,-8.708496564,
                -8.708496564,-8.708496564,-8.708496564,-8.708496564,-8.708496564,
                -8.708496564,0,-8.708496564,-8.708496564,-8.708496564,-8.708496564,
                -8.708496564,-8.708496564,-8.708496564,-8.708496564,-8.708496564,
                -8.708496564,-8.708496564,-8.708496564,0,-8.708496564,-8.708496564,
                -8.708496564,-8.708496564,-8.708496564,-8.708496564,-8.708496564,
                -8.708496564,-8.708496564,0,-8.708496564,-8.708496564,-8.708496564,
                -8.708496564,-8.708496564,-8.708496564,0,-8.708496564,-8.708496564,0,
                -8.708496564,-8.708496564,-8.708496564,-8.708496564,-8.708496564,
                -8.708496564,-8.708496564,-8.708496564,-8.708496564,-8.708496564,
                -8.708496564,0,-8.708496564,-8.708496564,-8.708496564,0,-8.708496564,
                -8.708496564,-8.708496564,-8.708496564,-8.708496564,-8.708496564,0,
                -8.708496564,-8.708496564,-8.708496564,-8.708496564,-8.708496564,
                -8.708496564,-8.708496564,-8.708496564,-8.708496564,-8.708496564,
                -8.708496564,0,-8.708496564,-8.708496564,0,-8.708496564,-8.708496564,
                -8.708496564,-8.708496564,-8.708496564,0,-8.708496564,-8.708496564,
                -8.708496564,-8.708496564,-7.541459043,0,-7.541459043,0,0,0,0,
                -7.541459043,-7.541459043,0,0,-7.541459043,-7.541459043,-7.541459043,
                -7.541459043,-7.541459043,0,-7.541459043,8.156584524,-7.541459043,0,
                0,8.156584524,-7.541459043,8.156584524,0,-7.541459043,0,-7.541459043,
                0,0,0,0,-7.541459043,-15.08291809,0,0,0,-7.541459043,-7.541459043,
                -8.156584524,8.156584524,0,0,-16.31316905,0,-16.31316905,0,8.156584524)
val<- data.frame(TYPE, MAX, STEM_PERCENT)
na.strings=c("",NA)
val<-subset(val,MAX!= "NA")

这是我用来运行混合效果分析的代码，它在 R v 3.4.4 中运行良好

library(lme4)
library(merTools)
mod3<-lmer(STEM_PERCENT~1+MAX+(1|TYPE)+(0+MAX|TYPE),data=val)

在这里，我收到以下错误消息，我从未在 R v 3.4.4 中收到过

Warning message:
In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model failed to converge with max|grad| = 0.107451 (tol = 0.002, component 1)

代码继续

mod4<-lmer(STEM_PERCENT~1+MAX+(1+MAX|TYPE),data=val) # also produces error not seen before


val2<-expand.grid(MAX=seq(0,120,length=1000),
                  TYPE=levels(val$TYPE))

val2$STEM_PERCENT<-predict(mod3,newdata=val2)

val3<-data.frame(MAX=seq(0,120,0.1))
val3$STEM_PERCENT<-predict(mod3,newdata=val3,re.form=~0)

CI <- cbind(val2, predictInterval(mod4, val2))


plot(val$MAX,val$STEM_PERCENT,pch=19,col=(1+as.integer((val$TYPE))),
     bty="l",cex=0.5,xaxt="n",las=1,ylim=c(-30,10),
     xlim=c(0,120),yaxt="n",xlab=NA,ylab=NA)

axis(side=2,tck=-0.02,at=seq(-30,10,10),cex.axis=0.5,
     font.axis=1,las=2,mgp=c(0,.5,0),labels=T)
axis(side=1,tck=-0.02,at=seq(0,120,20),
     cex.axis=0.5,labels=T,mgp=c(0,-0.1,0))


xv<-seq(0,120,0.01)
typea<-rep("A",length(xv))
yv<-predict(mod3,list(MAX=xv,TYPE=typea),type="response")

这是我在 R v 3.6.1 中收到错误消息的地方，而在 v 3.4.4 中没有。 yv的每个实例也是如此

Error in rep(0, nobs) : invalid 'times' argument

代码继续

lines(xv,yv,col="red",lwd=1.5)

typeb<-rep("B",length(xv))
yv<-predict(mod3,list(MAX=xv,TYPE=typeb),type="response")
lines(xv,yv,col="green",lwd=1.5)

typec<-rep("C",length(xv))
yv<-predict(mod3,list(MAX=xv,TYPE=typec),type="response")
lines(xv,yv,col="blue",lwd=1.5)

lines(val3$MAX,val3$STEM_PERCENT,lwd=2)

ssmin<-smooth.spline(CI$MAX,CI$lwr,df=4)
lines(ssmin,lty=2,col="black",lwd=1)

ssmax<-smooth.spline(CI$MAX,CI$upr,df=4)
lines(ssmax,lty=2,col="black",lwd=1)

这使得 R v 3.4.4 中的以下情节没有问题。

我希望能够在给定的 MAX 值附近找到上下置信区间值，例如当 MAX = 50 时，fit 的 upr 和 lwr 限制。对于upr 和 lwr 限制，它看起来应该分别在 -2 和 -12 左右。我还想知道以前没有收到的错误消息是怎么回事。

感谢您的帮助和建议，谢谢

【问题讨论】：

相关的很可能不是不同的R版本，而是lme4的不同版本。
newdata 参数被记录为期望一个 data.frame，实际上，yv<-predict(mod3, data.frame(MAX=xv,TYPE=typea), type="response") 工作得很好。如果您还没有，您可能应该研究 help("convergence") 关于“不收敛”警告（这不是错误）。
谢谢，@Roland。我通过在算法无法收敛之前用参数估计重新启动模型解决了收敛问题，因此：ss <- getME(mod3,c("theta","fixef")) new_mod <- update(mod3,start=ss)对于其他人，关于收敛问题有更多帮助here您可以就CI查询提供建议吗？跨度>

标签： r predict lme4 mixed-models confidence-interval

【解决方案1】：

我想出了如何确定 CI 值。希望其他人可以从解决方案中受益。像往常一样，这很简单。

首先，将 CI 舍入为 0 d.p

CI$MAX<-as.numeric(round(CI$MAX,0))

然后，确定 MAX = 50 处的平均上限和下限

low<-mean(CI$lwr[which(CI$MAX==50)])
high<-mean(CI$upr[which(CI$MAX==50)])

在绘图中添加要检查的线

abline(v = 50)
abline(h = low)
abline(h = high)

【讨论】：