Извлечь полосу предсказания из lme fit

У меня следующая модель

x  <- rep(seq(0, 100, by=1), 10)
y  <- 15 + 2*rnorm(1010, 10, 4)*x + rnorm(1010, 20, 100)
id <- NULL
for(i in 1:10){ id <- c(id, rep(i,101)) }
dtfr <- data.frame(x=x,y=y, id=id)
library(nlme)
with(dtfr, summary(     lme(y~x, random=~1+x|id, na.action=na.omit)))
model.mx <- with(dtfr, (lme(y~x, random=~1+x|id, na.action=na.omit)))
pd       <- predict( model.mx, newdata=data.frame(x=0:100), level=0)
with(dtfr, plot(x, y))
lines(0:100, predict(model.mx, newdata=data.frame(x=0:100), level=0), col="darkred", lwd=7)

введите описание изображения здесь

с predict и level=0 я могу построить средний ответ населения. Как я могу извлечь и построить 95% доверительные интервалы / диапазоны прогнозов из объекта nlme для всей популяции?


r prediction confidence-interval mixed-models
person ECII    schedule 16.01.2013    source источник
comment
хороший вопрос! Если вы понимаете, вы пытаетесь получить эквивалент этого curve(predict(model.lm, data.frame(x=x),interval ='confidence'),add=T), где model.lm, например, lm (y ~ x)   -  person agstudy    schedule 16.01.2013
comment
да. С нижним и верхним КИ.   -  person ECII    schedule 16.01.2013
comment
Я думаю, что даже для фильма получить его - это утомительно. есть intervals .lme функция, но она не дает группе уверенности ни на один балл.   -  person agstudy    schedule 16.01.2013
comment
intervals получает CI оценок / коэффициентов совпадений. Что мне нужно CI y для любого заданного x.   -  person ECII    schedule 16.01.2013
comment
на самом деле, @ECII, ты старался на собственном опыте ... Я имею в виду, что сам рассчитал полосу ..?   -  person agstudy    schedule 17.01.2013
comment
Что ты имеешь в виду под собой?   -  person ECII    schedule 21.01.2013
comment
Прочтите часто задаваемые вопросы о моделях r-sig-mixed.   -  person Roland    schedule 21.01.2013

Ответы (2)


arrow_upward
24
arrow_downward

Предупреждение. Прочтите эту ветку на моделях r-sig-mixed, прежде чем делать это. Будьте очень осторожны при интерпретации результирующего диапазона прогнозов.

Из r-sig - FAQ по смешанным моделям, адаптированный к вашему примеру:

set.seed(42)
x <- rep(0:100,10)
y <- 15 + 2*rnorm(1010,10,4)*x + rnorm(1010,20,100)
id<-rep(1:10,each=101)

dtfr <- data.frame(x=x ,y=y, id=id)

library(nlme)

model.mx <- lme(y~x,random=~1+x|id,data=dtfr)

#create data.frame with new values for predictors
#more than one predictor is possible
new.dat <- data.frame(x=0:100)
#predict response
new.dat$pred <- predict(model.mx, newdata=new.dat,level=0)

#create design matrix
Designmat <- model.matrix(eval(eval(model.mx$call$fixed)[-2]), new.dat[-ncol(new.dat)])

#compute standard error for predictions
predvar <- diag(Designmat %*% model.mx$varFix %*% t(Designmat))
new.dat$SE <- sqrt(predvar) 
new.dat$SE2 <- sqrt(predvar+model.mx$sigma^2)

library(ggplot2) 
p1 <- ggplot(new.dat,aes(x=x,y=pred)) + 
geom_line() +
geom_ribbon(aes(ymin=pred-2*SE2,ymax=pred+2*SE2),alpha=0.2,fill="red") +
geom_ribbon(aes(ymin=pred-2*SE,ymax=pred+2*SE),alpha=0.2,fill="blue") +
geom_point(data=dtfr,aes(x=x,y=y)) +
scale_y_continuous("y")
p1

сюжет

person Roland    schedule 21.01.2013
comment
+1 за то, что избавил меня от необходимости делать это или чувствовать себя виноватым за то, что не делаю этого. Обратите внимание (как указано в FAQ), что это учитывает только неопределенность фиксированных эффектов, обусловленных оценками дисперсии случайных эффектов и BLUP / условных режимов. - person Ben Bolker; 21.01.2013
comment
Было бы неплохо получить здесь ссылку из FAQ. Я помню, что довольно много раз изобретал это колесо заново. - person Dieter Menne; 21.01.2013
comment
Я не думаю, что профессор Бейтс включил бы это в predict.lme, но было бы неплохо, если бы какой-нибудь пакет мог предоставить эту функцию (конечно, с четкими предупреждениями о статистических ограничениях в файле справки). - person Roland; 22.01.2013
comment
Не могли бы вы немного объяснить, что происходит в newdat? Применим ли этот код для любого количества регрессоров? - person ECII; 22.01.2013
comment
@ Роланд: Я боюсь, что Дуглас Бейтс переключится в свой критический режим и спросит, зачем? Вы уверены, что ваши нестатистические клиенты правильно читают график? Обусловленность оценок - это чертовски уродливая концепция, которую нужно объяснять клиническим исследователям. - person Dieter Menne; 22.01.2013

arrow_upward
2
arrow_downward

Извините за возвращение к такой старой теме, но это может быть связано с комментарием здесь:

было бы неплохо, если бы какой-нибудь пакет мог предоставить эту функциональность

Эта функция включена в ggeffects-package, когда вы используете type = "re" (который затем будет включать дисперсии случайного эффекта, а не только остаточной дисперсии, которая, однако, одинакова в этом конкретном примере).

library(nlme)
library(ggeffects)

x  <- rep(seq(0, 100, by = 1), 10)
y  <- 15 + 2 * rnorm(1010, 10, 4) * x + rnorm(1010, 20, 100)
id <- NULL
for (i in 1:10) {
  id <- c(id, rep(i, 101))
}
dtfr <- data.frame(x = x, y = y, id = id)
m <- lme(y ~ x,
         random =  ~ 1 + x | id,
         data = dtfr,
         na.action = na.omit)

ggpredict(m, "x") %>% plot(rawdata = T, dot.alpha = 0.2)


ggpredict(m, "x", type = "re") %>% plot(rawdata = T, dot.alpha = 0.2)

Создано 18 июня 2019 г. пакетом REPEX (v0.3.0)

person Daniel    schedule 18.06.2019