用 lme4 数据绘制两个连续变量的交互作用

我正在尝试绘制 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))))

有任何想法吗？将不胜感激。

数据设置：

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)

Is ID really nested within day？从技术上讲，这表明个体 1 (ID=1）在第 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)