从 lme 拟合中提取预测带

Posted 2023-02-16

【中文标题】从 lme 拟合中提取预测带

【英文标题】：Extract prediction band from lme fit

【发布时间】：2012-12-30 19:05:44

【问题描述】：

我有以下型号

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 我可以绘制平均人口响应。如何从整个人口的 nlme 对象中提取和绘制 95% 置信区间/预测带？

【问题讨论】：

好问题！如果理解，您尝试拥有与此 curve(predict(model.lm, data.frame(x=x),interval ='confidence'),add=T) 等效的东西，其中 model.lm 例如是 lm(y~x)

是的。具有较低和较高的 CI。

我认为即使是 lm 也是一件苦差事。有intervals .lme 功能，但它并没有给乐队信心一分。

intervals 获取拟合的估计值/系数的 CI。对于任何给定的 x，我需要 y 的 CI。

其实@ECII你有没有努力过...我的意思是自己计算乐队..？

【参考方案1】：

很抱歉回到这么老的话题，但这可能会解决这里的评论：

如果某个包可以提供这个功能就好了

此功能包含在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)

^{由reprex package (v0.3.0) 于 2019 年 6 月 18 日创建}

【参考方案2】：

警告：在执行此操作之前请阅读this thread on r-sig-mixed models。解释结果预测带时要非常小心。

来自r-sig-mixed models 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

【讨论】：

+1 为我省去了这样做的麻烦或为不这样做而感到内疚。请注意（如常见问题解答中所述），这仅说明了以随机效应方差和 BLUP/条件模式的估计为条件的固定效应的不确定性

最好能从您的常见问题解答到此处进行交叉引用。我记得我曾多次重新发明过那个***。

我认为 Bates 教授不会将其包含在 predict.lme 中，但如果某些软件包可以提供此功能会很好（当然，在帮助文件中有关于统计限制的明确警告)。

大家介意解释一下newdat 中发生了什么吗？此代码是否适用于任意数量的回归器？

@Roland：我担心道格拉斯·贝茨会切换到他的批判模式并问“为什么？你确定你的非统计客户正确地阅读了图表吗？”。 “以估计为条件”是一个向临床研究人员解释的丑陋的概念。