Ошибки Imer и функция прогнозирования в R v 3.6.1 / получение доверительных интервалов из модели Imer для данной независимой переменной

Это два вопроса в одном, я надеюсь, что это нормально.

Во-первых, я пытаюсь получить значения доверительного интервала из объекта 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 без проблем.

Ошибки Imer и функция прогнозирования в R v 3.6.1 / получение доверительных интервалов из модели Imer для данной независимой переменной

Я хотел бы иметь возможность находить значения верхнего и нижнего доверительного интервала вокруг среднего значения при заданном значении MAX, например. upr и lwr пределы fit при MAX = 50. Похоже, что это должно быть где-то около -2 и -12 для upr и lwr пределов соответственно. Я также хотел бы знать, что происходит с сообщениями об ошибках, которые я раньше не получал.

Любая помощь и совет с благодарностью принимаются, спасибо

Релевантны, скорее всего, не разные версии R, а разные версии lme4.

Roland 18.07.2019 11:05

Документировано, что параметр newdata ожидает data.frame, и действительно, yv<-predict(mod3, data.frame(MAX=xv,TYPE=typea), type = "response") работает просто отлично. Если вы еще этого не сделали, вам, вероятно, следует изучить help("convergence") предупреждения об отсутствии конвергенции (которые не являются ошибками).

Roland 18.07.2019 11:12

Спасибо, @Роланд. Я решил проблему сходимости, перезапустив модели с оценками параметров до того, как алгоритм не сойдется, таким образом: ss <- getME(mod3,c("theta","fixef")) new_mod <- update(mod3,start=ss) Для других есть дополнительная помощь по проблеме сходимости здесь Можете ли вы дать совет по запросу CI?

A.Benson 18.07.2019 21:44
Стоит ли изучать PHP в 2023-2024 годах?
Стоит ли изучать PHP в 2023-2024 годах?
Привет всем, сегодня я хочу высказать свои соображения по поводу вопроса, который я уже много раз получал в своем сообществе: "Стоит ли изучать PHP в...
Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией
Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией
В JavaScript одним из самых запутанных понятий является поведение ключевого слова "this" в стрелочной и обычной функциях.
Приемы CSS-макетирования - floats и Flexbox
Приемы CSS-макетирования - floats и Flexbox
Здравствуйте, друзья-студенты! Готовы совершенствовать свои навыки веб-дизайна? Сегодня в нашем путешествии мы рассмотрим приемы CSS-верстки - в...
Тестирование функциональных ngrx-эффектов в Angular 16 с помощью Jest
В системе управления состояниями ngrx, совместимой с Angular 16, появились функциональные эффекты. Это здорово и делает код определенно легче для...
Концепция локализации и ее применение в приложениях React ⚡️
Концепция локализации и ее применение в приложениях React ⚡️
Локализация - это процесс адаптации приложения к различным языкам и культурным требованиям. Это позволяет пользователям получить опыт, соответствующий...
Пользовательский скаляр GraphQL
Пользовательский скаляр GraphQL
Листовые узлы системы типов GraphQL называются скалярами. Достигнув скалярного типа, невозможно спуститься дальше по иерархии типов. Скалярный тип...
0
3
142
1
Перейти к ответу Данный вопрос помечен как решенный

Ответы 1

Ответ принят как подходящий

Я понял, как определить значения 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)

lmer_solution

Другие вопросы по теме