Это два вопроса в одном, я надеюсь, что это нормально.
Во-первых, я пытаюсь получить значения доверительного интервала из объекта lmer
из пакета lme4
. Раньше я использовал 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, например. upr
и lwr
пределы fit
при MAX = 50. Похоже, что это должно быть где-то около -2 и -12 для upr
и lwr
пределов соответственно. Я также хотел бы знать, что происходит с сообщениями об ошибках, которые я раньше не получал.
Любая помощь и совет с благодарностью принимаются, спасибо
Документировано, что параметр newdata
ожидает data.frame, и действительно, yv<-predict(mod3, data.frame(MAX=xv,TYPE=typea), type = "response")
работает просто отлично. Если вы еще этого не сделали, вам, вероятно, следует изучить help("convergence")
предупреждения об отсутствии конвергенции (которые не являются ошибками).
Спасибо, @Роланд. Я решил проблему сходимости, перезапустив модели с оценками параметров до того, как алгоритм не сойдется, таким образом: ss <- getME(mod3,c("theta","fixef")) new_mod <- update(mod3,start=ss)
Для других есть дополнительная помощь по проблеме сходимости здесь Можете ли вы дать совет по запросу 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, а разные версии
lme4
.