使用 ggplot2 函数的混合效应模型残差图答案

【问题标题】：Graphs of the mixed effects model residuals using the ggplot2 function使用 ggplot2 函数的混合效应模型残差图
【发布时间】：2021-04-29 15:56:27
【问题描述】：

我正在尝试使用 ggplot2 函数绘制混合效应模型的残余效应。但是，在执行搜索后，我发现了一些可用的功能，但在我看来，对于功能 nlme 它们不起作用。

我打算制作的图表是以下示例的图表：

我最初尝试的计算例程如下，查看在ggplot2中执行函数时出现的错误。

library(splines)
library(ggplot2)
library(nlme)
library(gridExtra)

datanew1$DummyVariable = as.factor(datanew1$DummyVariable)
datanew1$Variable2 = as.factor(datanew1$Variable2)
datanew1$Variable3 = as.factor(datanew1$Variable3)
#############################################################################
############################## Model ########################################
#############################################################################
model <-  lme(Response~(bs(Variable1, df=3)) + DummyVariable,
                         random=~1|Variable2/Variable3, datanew1, method="REML")
completemodel <- update(model, weights = varIdent(form=~1|DummyVariable))

p1 <- qplot(.fitted, .resid, data = completemodel) +
  geom_hline(yintercept = 0) +
  geom_smooth(se = FALSE)

Erro: `data` must be a data frame, or other object coercible by `fortify()`, not an S3 object with class lme
Run `rlang::last_error()` to see where the error occurred.

p2 <- qplot(sample =.stdresid, data = completemodel, stat = "qq") + geom_abline()
grid.arrange(p1,p2)

Erro: `data` must be a data frame, or other object coercible by `fortify()`, not an S3 object with class lme
Run `rlang::last_error()` to see where the error occurred.
Além disso: Warning message:
`stat` is deprecated

我尝试制作图表的另一种方法是使用下面的函数，但我没有成功。

ggplot(completemodel, aes(.fitted, .resid)) + geom_point()

Erro: `data` must be a data frame, or other object coercible by `fortify()`, not an S3 object with class lme
Run `rlang::last_error()` to see where the error occurred.

【问题讨论】：

您似乎在数据导入/清理步骤中遗漏了一些代码，因为我无法理解该数据集中的列是什么（单个列有一个标题，带有内容“Variable2;Variable3;Response;Variable1;DummyVariable”，其他两个未标记）。但是，我怀疑您的问题是 nmle 正在创建一个自定义 S3 对象来保存其信息，您必须查看该对象以查看要绘制的数据在哪里。 ggplot2 似乎正在努力将您的数据强制转换为矩形格式 - 检查 str(completemodel) 是什么。
@Dubukay 在那里，我用必要的信息和缺少的样条包更新了代码。看到“bs”、“weights”、“form”、“random”是包 splines 和 nlme 的函数，而不是变量！
啊，欧洲的csv格式！我应该想尝试一下。看起来我的怀疑是正确的——数据存储在 S3 模型对象中。您的列名（.fitted、.resid）表明数据应该首先通过fortify 运行，但这不是为 NMLE 对象定义的，因此会出现错误。但是，如果您想自己计算拟合和残差，可以使用 completemodel$fitted 和 completemodel$residuals。
别担心。噢，那太糟糕了！！您是在告诉我，这些图表不可能在 ggplot2 中使用通过调整 nlme 函数获得的残差来制作吗？感谢您的信息！
@Dubukay 我刚刚发布了一个解决方案。

标签： r ggplot2 regression mixed-models nlme

【解决方案1】：

[解决方案]

library(splines)
library(ggplot2)
library(nlme)
library(gridExtra)

datanew1$DummyVariable = as.factor(datanew1$DummyVariable)
datanew1$Variable2 = as.factor(datanew1$Variable2)
datanew1$Variable3 = as.factor(datanew1$Variable3)


model <-  lme(Response~(bs(Variable1, df=3)) + DummyVariable,
              random=~1|Variable2/Variable3, datanew1, method="REML")
completemodel <- update(model, weights = varIdent(form=~1|DummyVariable))

df_model <- broom.mixed::augment(completemodel)
#> Registered S3 method overwritten by 'broom.mixed':
#>   method      from 
#>   tidy.gamlss broom
df_model[".stdresid"] <- resid(completemodel, type = "pearson")

p1 <- ggplot(df_model, aes(.fitted, .resid)) + 
  geom_point() +
  geom_hline(yintercept = 0) +
  geom_smooth(se=FALSE)

p2 <- ggplot(df_model, aes(sample = .stdresid)) +
  geom_qq() +
  geom_qq_line()

grid.arrange(p1,p2)
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

【讨论】：