2つの連続変数のlme4データとの対話のプロット

Question

Rの2つの連続変数間の相互作用をプロットしようとしています。しかし、私のデータはマルチレベルです（人は数日以内に入れ子になります）。私はそれをグラフにしています。私はネストされた構造を説明するために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) 

私は日内にネストされた私のスピンとREGの独立変数、疲労の従属変数、および人（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") 

私はスピンとregとの間に相互作用があるとします。私はそれをプロットするために、連続変数をカテゴリ変数に入れる必要があります。

私は連続変数の1つに基づいてカテゴリ変数を作成します。ここで私はスピンを選ぶ。 注：以下のコードが私の望むところに適しているかどうかはわかりません。標準エラーが発生する可能性がありますか？また、ネストされたデータ構造を考慮に入れていませんが、何をすべきかはわかりません。

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

大変感謝しています。

Answer 1

IDとdayの両方を反映するようにモデルを少し変更しました。また、いくつかの非常に間がある

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

：

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

も：これはどのように

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) 

は本当にday内ネストされたIDですか？技術的には、それは、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がモデルに含まれていないのはなぜですか？

g0 + geom_point(data=testdata) 

ここから必要なデータを取得するための最初の試みだ...私はなぜ、これは因子レベルを反転さかなりわからないんだけど

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

：予測データを設定し

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)