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

Question

这是两个合二为一的问题，希望没问题。

首先，我试图从 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 左右分别。我还想知道我以前没有收到的错误消息是怎么回事。

感谢收到的任何帮助和建议，谢谢

Answer 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)

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 variable

r

predict

lme4

confidence-interval

mixed-models