绘制两个连续变量与 lme4 数据的交互答案

【问题标题】：Plotting interaction of two continuous variables with lme4 data绘制两个连续变量与 lme4 数据的交互
【发布时间】：2016-07-10 18:51:12
【问题描述】：

我试图在 R 中绘制两个连续变量之间的交互。但是，我的数据是多级的（人们嵌套在几天内），所以我在绘制图表时需要考虑数据的嵌套结构。我使用 lme4 库分析我的数据以解释嵌套结构，但我很难弄清楚如何绘制它。

## example data
spin = runif(600, 1, 24)
reg = runif(600, 1, 15)
ID = rep(c("1","2","3","4","5", "6", "7", "8", "9", "10"))
day = rep(1:30, each = 10)
testdata <- data.frame(
  spin, reg, ID, day)
testdata$fatigue <- testdata$spin*testdata$reg/10*rnorm(30, mean=3, sd=2)

在这里，我有自变量旋转和调节、疲劳因变量以及嵌套在几天内的人 (ID)。我在下面运行我的模型。

## running my multilevel model with lme4
library(lme4)
m1 <- lmer(fatigue ~ spin * reg + ( 1 | ID), data = testdata, REML = T)
(m1)
confint(m1, test = "Chisq")

假设我在 spin 和 reg 之间有交互。我需要将我的连续变量放入一个分类变量中才能绘制它。

所以我根据我的一个连续变量创建分类变量。在这里我选择旋转。注意：不确定下面的代码是否完全适合我想要的。可能必须做标准错误？也没有考虑我的嵌套数据结构，但不知道该怎么做。

x <- mean(testdata$spin, na.rm = T)
print(x)
y <- sd(testdata$spin, na.rm = T)
print(y)

testdata$SpinLevel[testdata$spin > x+y] <- "High"
testdata$SpinLevel[testdata$spin > x-y & testdata$spin <= x+y] <- "Mean"
testdata$SpinLevel[testdata$spin <= x-y] <- "Low"

rm(x,y)

根据我在网上找到的内容，我可以创建一个基本的情节来展示效果。但是没有考虑嵌套结构（人——变量ID——在几天内嵌套）。

library(ggplot2)
ggplot(testdata,aes(reg,fatigue,linetype=SpinLevel))+
  geom_smooth(method="lm",se=FALSE)

这个 ggplot 可以很好地解释基本效果，但线条可能会倾斜，因为它们没有考虑到我的数据的嵌套结构（几天内的人）。

我还可以使用效果库绘制我的模型。这确实考虑了嵌套结构。除了图表不漂亮而且是四分位数，而且很难解释。我希望它是高、平均和低，并且都在同一个图表上。但我不知道该怎么做。

library(effects)
plot(effect("spin*reg", m1), grid=TRUE, labels = T,
  xlevels=list(spin=quantile(testdata$spin, seq(0, 1, 0.25))))

有什么想法吗？将不胜感激。

【问题讨论】：

我在模型规范中看不到任何内容表明嵌套在 day 中。
stackoverflow.com/questions/9447329/…
我会推荐 broom 包作为 coefplot2 的替代品 ...但我认为 OP 想要一个效果图，而不是系数图 ...
抱歉，@42，我想我在统计意义上使用“嵌套”一词——当您的数据来自几天内响应的同一个人时，您不能假设响应是独立的。通过使用样本人，我打破了独立观察的统计假设，这就是我使用多级建模的原因。我在我的代码中包含了day，只是为了表明我在多天内使用了多个人，因为这就是我的数据的格式。 @BenBolker——你说得对，我确实想要效果图。
因此您忽略了相邻或附近日期之间潜在的自相关。您的“日子”是否如此相隔太远，以至于这在生理上具有很好的意义？

标签： r plot interaction lme4 multi-level

【解决方案1】：

我稍微更改了模型以使其同时反映ID 和day。

这个怎么样：

## example data
spin = runif(600, 1, 24)
reg = runif(600, 1, 15)
ID = rep(c("1","2","3","4","5", "6", "7", "8", "9", "10"))
day = rep(1:30, each = 10)
testdata <- data.frame(
spin, reg, ID, day)
testdata$fatigue <- testdata$spin*testdata$reg/10*rnorm(30, mean=3, sd=2)

## running my multilevel model with lme4
library(lme4)
m1 <- lmer(fatigue ~ spin * reg + ( 1 | ID/day), data = testdata, REML = T)
(m1)
confint(m1, test = "Chisq")

x <- mean(testdata$spin, na.rm = T)
print(x)
y <- sd(testdata$spin, na.rm = T)
print(y)

testdata$SpinLevel[testdata$spin > x+y] <- "High"
testdata$SpinLevel[testdata$spin > x-y & testdata$spin <= x+y] <- "Mean"
testdata$SpinLevel[testdata$spin <= x-y] <- "Low"

rm(x,y)

require(multicomp)
mp <- as.data.frame(confint(glht(m1))$confint)
tmp$Comparison <- rownames(tmp)
ggplot(tmp, aes(x = Comparison, y = Estimate, ymin = lwr, ymax = upr)) + geom_errorbar() + geom_point()

# or
library(multcomp)
tmp <- as.data.frame(confint(glht(m1))$confint)
tmp$Comparison <- rownames(tmp)
ggplot(tmp, aes(x = Comparison, y = Estimate, ymin = lwr, ymax = upr)) + geom_errorbar() + geom_point()

还有：

install.packages("coefplot2", # from this crackpot R coder named Bolker
                 repos="http://www.math.mcmaster.ca/bolker/R",
                 type="source") # I think he died a few years back
                                # jk Ben
library(coefplot2)
coefplot2(m1)
ggplot(tmp, aes(x = Comparison, y = Estimate, ymin = lwr, ymax = upr)) + geom_errorbar() + geom_point()

一个叫Wes here的人的回答中也有一些非常有趣的彩色图。

【讨论】：

【解决方案2】：

数据设置：

set.seed(101)
spin = runif(600, 1, 24)
reg = runif(600, 1, 15)
ID = rep(c("1","2","3","4","5", "6", "7", "8", "9", "10"))
day = rep(1:30, each = 10)
testdata <- data.frame(spin, reg, ID, day)
testdata$fatigue <- testdata$spin*testdata$reg/10*rnorm(30, mean=3, sd=2)

ID 真的嵌套在day 内吗？从技术上讲，这表明在第 1 天测量的个体 1 (ID=1) 代表了与ID=1 在第 2 天测量的不同人...？

library(lme4)
m1 <- lmer(fatigue ~ spin * reg + ( 1 | ID),
           data = testdata, REML = TRUE)
confint(m1, method = "Wald", parm="beta_")
## instead of test="Chisq", which doesn't work
##                    2.5 %    97.5 %
## (Intercept) -13.44726318 7.4959080
## spin         -0.04751327 1.2328254
## reg          -0.86763792 1.1550787
## spin:reg      0.11263238 0.2541709

为什么模型中没有day...？

设置预测数据：

## midpoints of bin
 spinvals <- quantile(testdata$spin,seq(0,1,length=5))[2:4]
 pframe <- with(testdata,
           expand.grid(ID=unique(ID),
                       reg=seq(min(reg),max(reg),length.out=51),
                       spin=spinvals))
 pframe$fatigue <- predict(m1,newdata=pframe)
 pframe$spinFac <- factor(pframe$spin,levels=spinvals)
 ## explicit factor() to prevent alphabetization of levels

 library(ggplot2); theme_set(theme_bw())
 g0 <- ggplot(pframe,aes(reg,fatigue,colour=spinFac))+
     geom_line(aes(group=interaction(spinFac,ID)))

 ## bins for cutting testdata into 3 levels (min, 0.33,0.66, max)
 ## label bins by midpoints
 spincuts <- quantile(testdata$spin,seq(0,1,length=4))
 testdata$spinFac <- cut(testdata$spin,
            spincuts,labels=spinvals)

我不太清楚为什么这会翻转因子水平......

 g0 + geom_point(data=testdata)

这是从effects 对象中提取所需数据的初步尝试：

library(effects)
ee <- effect("spin*reg", m1,
   xlevels=list(spin=spinvals))
eedat <- with(ee,data.frame(x,fatigue=fit,lwr=lower,upr=upper))
ggplot(eedat,aes(x=reg,y=fatigue,colour=factor(spin)))+
    geom_line()+
    geom_ribbon(aes(group=spin,ymin=lwr,ymax=upr),colour=NA,
                            alpha=0.4)

【讨论】：