Плотное взаимодействие двух непрерывных переменных с данными lme4

Question

Я пытаюсь построить взаимодействие между двумя непрерывными переменными в 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") 

Предположим, что у меня есть взаимодействие между спином и регистром. Мне нужно включить переменную непрерывную в категориальную переменную для ее построения.

Поэтому я создаю категориальную переменную, основанную на одной из моих переменных. Здесь я выбираю спин. Примечание: не уверен, что этот код ниже подходит для того, что я хочу. Должна ли быть стандартная ошибка? Также не учитывает мою вложенную структуру данных, но не уверен, что делать иначе.

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) 

Основываясь на том, что я нашел в Интернете, я могу создать основной сюжет, чтобы показать эффекты. НО не учитывает вложенную структуру (люди - идентификаторы переменных - вложены в течение дней).

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)))) 

Любые идеи? Было бы очень признательно.

Answer 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.

Answer 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 (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)