如何在 lme 中混合模型的预测函数中包含异方差性答案

【问题标题】：How to include heteroscedasticity in the predict function of a mixed model in lme如何在 lme 中混合模型的预测函数中包含异方差性
【发布时间】：2015-12-16 22:02:58
【问题描述】：

大家好，

我从“nlme”和“ggplot2”绘制了一个混合模型，我想考虑模型中我的类别变量的不同方差（weights=varIdent(form = ~1 | category)）。当我检查情节时，我意识到情节中的模型预测不包括异方差性，尽管模型确实如此。如何调整预测函数，使其包含异方差？

创建一些数据：

no <- 1:20
type <- factor("A","B")
set.seed(1)
baseline <- rnorm(20, 30, 10)
set.seed(2)
event <- rnorm(20, 60, 30)
set.seed(3)
postevent1 <- rnorm(20, 40, 20)
set.seed(4)
postevent2 <- rnorm(20, 30, 10)
set.seed(5)
postevent3 <- rnorm(20, 20, 5)
data <- data.frame(no, baseline, event, postevent1, postevent2, postevent3)
data$type[data$no %% 2 == 0] <- "A"
data$type[data$no %% 2 != 0] <- "B"
data$event[data$type == "A"] <- data$event[data$type == "A"] -2
data$no <- factor(data$no)
data$type <- factor(data$type)

将数据转换为长格式

library(dplyr)

long <- data %>% gather(key = category, value = measure, -no, -type)

创建模型

library(nlme)
model <- lme(measure ~ category*type, random = ~ 1|no, data= long,   
weights=varIdent(form = ~1 | category), method = "REML")
new <- long %>% select(-measure)
long$pred <- predict(model, newdata = new)

画出来

library(ggplot2)
ggplot(long, aes(x=category, y =  measure, colour = type,    
group=interaction(type,category))) +
geom_point() +
facet_grid(~ type) +
geom_line(aes(y=pred), size=0.8, colour = "black") +
theme(axis.text.x = element_text(angle=30, hjust=1, vjust=1))

【问题讨论】：

如果您愿意切换到lme4 而不是nlme，sim 函数可能会满足您的需求。
A 和 B 的方差相同，但在类别变量中方差不同。 lme 模型考虑了这一点，但 lme 模型的预测函数没有。
亲爱的格雷戈尔，你是对的！我已经改了标题。
我不知道如何解释 lme4 中的异方差性。我在互联网上找到的评论说重量参数在 lmer link 中没有正常工作
再一次你的权利。我刚刚修好了。

标签： r ggplot2 predict mixed-models nlme

【解决方案1】：

这是我使用“AICcmodavg”包的解决方案。如果有人找到更好的解决方案，请告诉我。

预测值

library(AICcmodavg)
pred <- predictSE.lme(model, newdata = new, print.matrix = TRUE)
pred <- data.frame(pred)
newlong <- data.frame(long, pred)

用预测的 95% 置信区间绘制平均值

ggplot(newlong, aes(x= category, y= measure)) + 
geom_point(aes(x = category, y = fit)) + 
facet_grid(~ type) + 
geom_point(aes(colour = type), alpha = 0.5, size = 3)  + 
geom_errorbar(aes(ymax=fit + se.fit*1.96, ymin=fit - se.fit*1.96), 
position=position_dodge(0.9), width=0.25) +
theme(axis.text.x = element_text(angle=30, hjust=1, vjust=1))

【讨论】：