Построение графика взаимодействия двух непрерывных переменных с данными 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)

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

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

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


r interaction multi-level plot lme4
person sk17    schedule 10.07.2016    source источник
comment
Я не вижу в спецификации модели ничего, что указывало бы на вложение в day.   -  person IRTFM    schedule 10.07.2016
comment
stackoverflow .com / questions / 9447329 /   -  person Hack-R    schedule 10.07.2016
comment
Я бы порекомендовал пакет broom в качестве замены coefplot2 ... но я думаю, что OP хочет график эффектов, а не график коэффициентов ...   -  person Ben Bolker    schedule 10.07.2016
comment
Извините, @ 42, я полагаю, что я использую термин «вложенность» в статистическом смысле - вы не можете предполагать, что ответы независимы, если ваши данные принадлежат тем же людям, которые отвечают в течение нескольких дней. Используя выборку людей, я нарушаю статистическое предположение о независимых наблюдениях, поэтому я использую многоуровневое моделирование. Я включил day в свой код, чтобы просто показать, что я использую несколько человек в течение нескольких дней, потому что именно так форматируются мои данные. @ BenBolker - вы правы, мне нужны графики эффектов.   -  person sk17    schedule 11.07.2016
comment
Таким образом, вы игнорируете потенциальную автокорреляцию по соседним или близлежащим дням. Ваши дни настолько разделены, что это имеет физиологический смысл?   -  person IRTFM    schedule 11.07.2016
comment
@ 42, в этом примере модели я учитываю только отсутствие независимости в отношении ID. Я мог бы запустить модель и учесть day, если я считаю, что это имеет значение для моих гипотез, но меня здесь нет.   -  person sk17    schedule 13.07.2016

Ответы (2)


arrow_upward
1
arrow_downward

Я немного изменил модель, чтобы она отражала и 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()

В ответе человека по имени Уэс здесь.

person Hack-R    schedule 10.07.2016

arrow_upward
1
arrow_downward

Настройка данных:

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 измеряли во второй день ...?

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)
person Ben Bolker    schedule 10.07.2016