クロスレベル交互作用 - 井出草平の研究ノート

クロスレベル交互作用とは、階層データ (マルチレベルデータ) において、異なるレベルの変数間の交互作用を指す。つまり、あるレベルの変数 (例：生徒レベル) の効果が、別のレベルの変数 (例：学校レベル) の水準によって異なる場合、クロスレベル交互作用が存在すると言える。

例:

生徒の「勉強時間」の効果 (成績への影響) が、学校の「教育の質」によって異なる場合、これはクロスレベル交互作用の一例となる。勉強時間の効果は生徒レベル (レベル1) の変数であり、教育の質は学校レベル (レベル2) の変数である。教育の質が高い学校では、勉強時間の効果が大きくなる (勉強時間を増やすほど成績が上がりやすい) かもしれないし、逆に教育の質が低い学校では、勉強時間の効果が小さくなる (いくら勉強時間を増やしても成績が上がりにくい) かもしれない。

1. データの読み込みと確認:

# パッケージの読み込み
library(mlmRev)
library(lme4)

# データの読み込み
data(Exam)

# schgend を因子型に変換
#Exam$schgend <- factor(Exam$schgend)

# データの確認
head(Exam)
str(Exam)

データ。

  school   normexam schgend    schavg      vr     intake   standLRT sex type student
1      1  0.2613242   mixed 0.1661752 mid 50% bottom 25%  0.6190592   F  Mxd     143
2      1  0.1340672   mixed 0.1661752 mid 50%    mid 50%  0.2058022   F  Mxd     145
3      1 -1.7238820   mixed 0.1661752 mid 50%    top 25% -1.3645760   M  Mxd     142
4      1  0.9675862   mixed 0.1661752 mid 50%    mid 50%  0.2058022   F  Mxd     141
5      1  0.5443412   mixed 0.1661752 mid 50%    mid 50%  0.3711052   F  Mxd     138
6      1  1.7348992   mixed 0.1661752 mid 50% bottom 25%  2.1894372   M  Mxd     155

データの説明:

Exam データは、ロンドンの中等学校の生徒の試験成績データ。
4,059人の生徒が65の学校に所属している (階層データ)。
主な変数:
- normexam: 標準化された試験得点 (連続変数) - 従属変数
- standLRT: 標準化された以前のテスト得点 (連続変数) - 生徒レベル (レベル1) の変数
- schgend: 学校の性別 (1: 男子校, 2: 女子校, 3: 共学) - 学校レベル (レベル2) の変数
- avslrt: 学校の平均 standLRT - 学校レベル (レベル2) の変数
- school: 学校 ID - レベル2変数
- student: 生徒 ID - レベル1変数

数式的な表現:

最もシンプルなモデルは以下。

$\displaystyle \text{レベル 1（例：生徒レベル）:} \quad y_{ij} = \beta_{0j} + \beta_{1j} x_{ij} + e_{ij}$

$\displaystyle \text{レベル 2（例：学校レベル）:} \quad \beta_{0j} = \gamma_{00} + \gamma_{01} z_{j} + u_{0j}$

$\displaystyle \beta_{1j} = \gamma_{10} + \gamma_{11} z_{j} + u_{1j}$

ここで:

$y_{ij}$ ：生徒 $i$ 、学校 $j$ の従属変数（例：試験成績）
$x_{ij}$ ：生徒 $i$ 、学校 $j$ のレベル1変数（例：勉強時間）
$z_{j}$ ：学校 $j$ のレベル2変数（例：教育の質）
$\beta_{0j}$ ：学校 $j$ における切片
$\beta_{1j}$ ：学校 $j$ における $x_{ij}$ の傾き（効果）
$\gamma_{00}$ ：固定効果の係数
$\gamma_{01}$ ：固定効果の係数
$\gamma_{10}$ ：固定効果の係数
$\gamma_{11}$ ：固定効果の係数
$u_{0j}$ ：学校 $j$ のランダム切片
$u_{1j}$ ：学校 $j$ の $x_{ij}$ のランダム傾き
$e_{ij}$ ：生徒 $i$ 、学校 $j$ の残差

上記の式で、 $\gamma_{11}$ はクロスレベル交互作用の係数を表す。 $\gamma_{11}$ が有意であれば、レベル1変数 $x_{ij}$ の効果（ $\beta_{1j}$ ）がレベル2変数 $z_{j}$ の値によって異なることを意味する。

検討可能なクロスレベル交互作用の例:

standLRT (生徒の以前のテスト得点) の効果が、schgend (学校の性別) によって異なるか？
- つまり、生徒の以前のテスト得点が試験成績に与える影響が、共学、男子校、女子校で異なるかどうかを検討する。
- 例えば、男子校では、以前のテスト得点が高い生徒ほど試験成績が高くなる傾向が強い一方、共学ではその傾向が弱い、といった可能性があるかもしれない。
- モデル: normexam ~ standLRT * schgend + schavg + (1 | school)
- 交互作用項: standLRT:schgend
standLRT (生徒の以前のテスト得点) の効果が、schavg (学校の平均的な以前のテスト得点) によって異なるか？
- つまり、生徒の以前のテスト得点が試験成績に与える影響が、学校全体の学力レベルによって異なるかどうかを検討する。
- 例えば、学校の平均 standLRT が高い学校では、生徒の standLRT の効果が小さくなる (生徒間の差が小さくなる) 一方、学校の平均 standLRT が低い学校では、生徒の standLRT の効果が大きくなる (生徒間の差が大きくなる) といった可能性があるかもしれない。
- モデル: normexam ~ standLRT * schavg + schgend + (1 | school)
- 交互作用項: standLRT:schavg

3. クロスレベル交互作用の解釈:

クロスレベル交互作用が有意である場合、レベル1変数の効果はレベル2変数の水準によって異なることを意味する。解釈の際には、以下の点に注意する必要がある。

交互作用項の係数の符号と大きさ: 係数の符号は、レベル2変数の値が大きくなるにつれて、レベル1変数の効果がどのように変化するかを示か。係数の大きさは、その変化の程度を表す。
レベル2変数の水準: レベル2変数がカテゴリ変数の場合、各水準におけるレベル1変数の効果を比較する必要がある。
グラフの活用: 交互作用の効果を視覚的に理解するために、グラフを作成することがのぞましい。レベル2変数の各水準における、レベル1変数と従属変数の関係をプロットするなど。

例：standLRT と schgend のクロスレベル交互作用の解釈:

仮に、standLRT:schgend の交互作用項が有意であり、男子校 (schgend = "boys") における standLRT の係数が、共学 (schgend = "mixed") よりも大きいとする。この場合、以下のように解釈できる。

男子校では、共学に比べて、生徒の以前のテスト得点 (standLRT) が高いほど、試験成績 (normexam) が高くなる傾向が強い。
言い換えると、男子校では、以前のテスト得点に基づく生徒間の成績差が、共学よりも大きくなる可能性がある。

2. 分析の目的とモデル:

ここでは、以下の2つのクロスレベル交互作用に注目して分析を行う。

生徒の以前のテスト得点 (standLRT) の効果が、学校の性別 (schgend) によって異なるか？
生徒の以前のテスト得点 (standLRT) の効果が、学校の平均的な以前のテスト得点 (avslrt) によって異なるか？

これらの交互作用を検討するために、以下のモデルを推定する。

モデル 1: ベースラインモデル (ランダム切片モデル)

model_base <- lmer(normexam ~ standLRT + schgend + schavg + (1 | school), data = Exam)
summary(model_base)

モデル 2: standLRT と schgend のクロスレベル交互作用を追加

model_int_gend <- lmer(normexam ~ standLRT * schgend + schavg + (1 | school), data = Exam)
summary(model_int_gend)

モデル 3: standLRT と schavg のクロスレベル交互作用を追加

model_int_schavg <- lmer(normexam ~ standLRT * schavg + schgend + (1 | school), data = Exam)
summary(model_int_schavg)

モデル 4: 両方のクロスレベル交互作用を追加

model_int_both <- lmer(normexam ~ standLRT * schgend + standLRT * schavg + (1 | school), data = Exam)
summary(model_int_both)

モデル 5: standLRTのランダムスロープを追加

model_int_both_slope <- lmer(normexam ~ standLRT * schgend + standLRT * schavg + (standLRT | school), data = Exam)
summary(model_int_both_slope)

3. モデルの推定と結果の解釈:

上記のコードを実行すると、各モデルの推定結果が得られる。結果の解釈のポイントは以下の通り。

固定効果の係数:
- standLRT の係数は、生徒の以前のテスト得点が1単位高いと、試験得点がどのくらい高くなるかを表す。
- schgend の係数は、学校の性別が共学と比較して男子校や女子校である場合に試験得点がどう変化するかを表す。(数値ではなくカテゴリーとして扱う必要があるので、schgendをfactorに変換)
- avslrt の係数は、学校の平均的な以前のテスト得点が1単位高いと、生徒の試験得点がどのくらい高くなるかを表す。
- 交互作用項 (standLRT:schgend, standLRT:avslrt) の係数は、standLRT の効果が、schgend や avslrt の水準によってどのくらい異なるかを表す。
ランダム効果の分散:
- (1 | school) は学校ごとの切片のばらつき (分散) を表す。
- (standLRT | school) は学校ごとの standLRT の傾きのばらつき (分散) を表す。
モデルの比較:
- anova() 関数を用いて、異なるモデルを比較することができる。例えば、anova(model_base, model_int_gend) は、ベースラインモデルと standLRT と schgend の交互作用を含むモデルを比較する。
- AIC や BIC などの情報量基準を用いて、モデルの適合度を比較する。

4. クロスレベル交互作用の解釈の例:

model_int_gend において、standLRT:schgend の交互作用項が有意である場合:
- 生徒の以前のテスト得点 (standLRT) が試験得点に与える影響は、学校の性別 (schgend) によって異なることを意味する。
- 例えば、男子校では standLRT の効果が共学よりも大きい (つまり、以前のテスト得点が高い生徒ほど、試験得点が高くなる傾向が強い) 可能性がある。
model_int_avslrt において、standLRT:avslrt の交互作用項が有意である場合:
- 生徒の以前のテスト得点 (standLRT) が試験得点に与える影響は、学校の平均的な以前のテスト得点 (avslrt) によって異なることを意味する。
- 例えば、学校の平均 avslrt が高いほど、standLRT の効果が小さくなる (つまり、以前のテスト得点が高い生徒と低い生徒の差が小さくなる) 可能性がある。

モデルの要約

summary(model_base)
summary(model_int_gend)
summary(model_int_schavg)
summary(model_int_both)
summary(model_int_both_slope)

モデルの比較

anova(model_base, model_int_gend)
anova(model_base, model_int_schavg)
anova(model_base, model_int_both)
anova(model_int_gend, model_int_both)
anova(model_int_schavg, model_int_both)
anova(model_int_both, model_int_both_slope)

結果。

> anova(model_base, model_int_gend)
 
refitting model(s) with ML (instead of REML)
Data: Exam
Models:
model_base: normexam ~ standLRT + schgend + schavg + (1 | school)
model_int_gend: normexam ~ standLRT * schgend + schavg + (1 | school)
               npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
model_base        7 9352.8 9397.0 -4669.4   9338.8                     
model_int_gend    9 9355.9 9412.7 -4669.0   9337.9 0.9365  2     0.6261

結果。

> anova(model_base, model_int_schavg)
refitting model(s) with ML (instead of REML)
Data: Exam
Models:
model_base: normexam ~ standLRT + schgend + schavg + (1 | school)
model_int_schavg: normexam ~ standLRT * schavg + schgend + (1 | school)
                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
model_base          7 9352.8 9397.0 -4669.4   9338.8                         
model_int_schavg    8 9335.8 9386.2 -4659.9   9319.8 19.082  1  1.252e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

結果。

> anova(model_base, model_int_both)
refitting model(s) with ML (instead of REML)
Data: Exam
Models:
model_base: normexam ~ standLRT + schgend + schavg + (1 | school)
model_int_both: normexam ~ standLRT * schgend + standLRT * schavg + (1 | school)
               npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)    
model_base        7 9352.8 9397.0 -4669.4   9338.8                        
model_int_both   10 9338.2 9401.3 -4659.1   9318.2  20.6  3  0.0001275 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

結果。

> anova(model_int_gend, model_int_both)
refitting model(s) with ML (instead of REML)
Data: Exam
Models:
model_int_gend: normexam ~ standLRT * schgend + schavg + (1 | school)
model_int_both: normexam ~ standLRT * schgend + standLRT * schavg + (1 | school)
               npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
model_int_gend    9 9355.9 9412.7 -4669.0   9337.9                         
model_int_both   10 9338.2 9401.3 -4659.1   9318.2 19.663  1  9.236e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

結果。

> anova(model_int_schavg, model_int_both)
refitting model(s) with ML (instead of REML)
Data: Exam
Models:
model_int_schavg: normexam ~ standLRT * schavg + schgend + (1 | school)
model_int_both: normexam ~ standLRT * schgend + standLRT * schavg + (1 | school)
                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
model_int_schavg    8 9335.8 9386.2 -4659.9   9319.8                     
model_int_both     10 9338.2 9401.3 -4659.1   9318.2 1.5172  2     0.4683

結果。

> anova(model_int_both, model_int_both_slope)
refitting model(s) with ML (instead of REML)
Data: Exam
Models:
model_int_both: normexam ~ standLRT * schgend + standLRT * schavg + (1 | school)
model_int_both_slope: normexam ~ standLRT * schgend + standLRT * schavg + (standLRT | school)
                     npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
model_int_both         10 9338.2 9401.3 -4659.1   9318.2                         
model_int_both_slope   12 9314.9 9390.6 -4645.4   9290.9 27.341  2  1.156e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

モデル比較の結果のまとめ

比較対象モデル	npar	AIC	BIC	logLik	deviance	Chisq	Df	Pr(>Chisq)
`model_base` vs `model_int_gend`	7 vs 9	9352.8 vs 9355.9	9397.0 vs 9412.7	-4669.4 vs -4669.0	9338.8 vs 9337.9	0.9365	2	0.6261
`model_base` vs `model_int_schavg`	7 vs 8	9352.8 vs 9335.8	9397.0 vs 9386.2	-4669.4 vs -4659.9	9338.8 vs 9319.8	19.082	1	1.252e-05 ***
`model_base` vs `model_int_both`	7 vs 10	9352.8 vs 9338.2	9397.0 vs 9401.3	-4669.4 vs -4659.1	9338.8 vs 9318.2	20.6	3	0.0001275 ***
`model_int_gend` vs `model_int_both`	9 vs 10	9355.9 vs 9338.2	9412.7 vs 9401.3	-4669.0 vs -4659.1	9337.9 vs 9318.2	19.663	1	9.236e-06 ***
`model_int_schavg` vs `model_int_both`	8 vs 10	9335.8 vs 9338.2	9386.2 vs 9401.3	-4659.9 vs -4659.1	9319.8 vs 9318.2	1.5172	2	0.4683
`model_int_both` vs `model_int_both_slope`	10 vs 12	9338.2 vs 9314.9	9401.3 vs 9390.6	-4659.1 vs -4645.4	9318.2 vs 9290.9	27.341	2	1.156e-06 ***

結果の解釈

model_base vs model_int_gend:
- standLRT と schgend の交互作用を追加した model_int_gend は、model_base と比較して、統計的に有意な改善を示していない (p = 0.6261)。
- これは、生徒の以前のテスト得点 (standLRT) の効果が、学校の性別 (schgend) によって大きくは異ならない可能性を示唆している。
- ただし、AICとBICは、model_baseの方がわずかに低いため、より良いモデルである可能性がある。
model_base vs model_int_schavg:
- standLRT と schavg の交互作用を追加した model_int_schavg は、model_base と比較して、統計的に有意な改善を示している (p = 1.252e-05)。
- これは、生徒の以前のテスト得点 (standLRT) の効果が、学校の平均的な以前のテスト得点 (schavg) によって異なることを示唆している。
- AICとBICも、model_int_schavgの方が低いため、より良いモデルであると考えられる。
model_base vs model_int_both:
- 両方の交互作用を追加した model_int_both は、model_base と比較して、統計的に有意な改善を示している (p = 0.0001275)。
- これは、standLRT の効果が schgend と schavg の両方によって異なる可能性を示唆している。
- ただし、AICはmodel_int_schavgの方が低く、BICはmodel_int_bothの方が低いため、どちらが良いモデルかは判断が分かれる。
model_int_gend vs model_int_both:
- model_int_both は model_int_gend と比較して、統計的に有意な改善を示している (p = 9.236e-06)。
- これは、standLRT と schgend の交互作用に加えて、standLRT と schavg の交互作用もモデルに含めることで、適合度が向上することを示唆している。
- AICとBICをふまえると、model_int_bothの方が良いモデルであると考えられる。
model_int_schavg vs model_int_both:
- model_int_both は model_int_schavg と比較して、統計的に有意な改善を示していない (p = 0.4683)。
- これは、standLRT と schavg の交互作用に加えて、standLRT と schgend の交互作用をモデルに含めても、適合度はあまり向上しないことを示唆している。
- AICとBICをふまえると、model_int_schavgの方が良いモデルである可能性がある。
model_int_both vs model_int_both_slope:
- standLRT のランダム傾きを追加した model_int_both_slope は、model_int_both と比較して、統計的に有意な改善を示している (p = 1.156e-06)。
- これは、学校ごとに standLRT の効果が異なることを示唆している。
- AICとBICも、model_int_both_slopeの方が低いため、より良いモデルであると考えられる。

モデル評価まとめ

上記の比較結果から、model_int_both_slope が最も良いモデルである可能性が高い と考えられる。このモデルは、以下の特徴を持っている。

standLRT と schgend のクロスレベル交互作用を含む。
standLRT と schavg のクロスレベル交互作用を含む。
standLRT のランダム傾きを含む (学校ごとに standLRT の効果が異なる)。

プロット

傾きが異なれば、レベル1変数 (standLRT) の効果がレベル2変数 (schgend または schavg) によって変化していることを意味する。
回帰直線が交差している場合は、強い交互作用効果があることを示唆する。

model_base のプロット（交互作用なし）

# 予測値の取得
Exam$predicted_base <- predict(model_base)

# プロットの作成
ggplot(Exam, aes(x = standLRT, y = normexam, color = schgend)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = "lm", se = FALSE, aes(y = predicted_base)) +
  labs(
    title = "ベースラインモデル: standLRT の単純効果",
    x = "生徒の以前のテスト得点 (standLRT)",
    y = "標準化された試験得点 (normexam)",
    color = "学校の性別 (schgend)"
  ) +
  theme_minimal()

model_int_gend: standLRT × schgend の交互作用　

# パッケージの読み込み
library(ggplot2)
library(lme4)

# 交互作用モデルの適用
model_int_gend <- lmer(normexam ~ standLRT * schgend + schavg + (1 | school), data = Exam)

# 予測値の取得
Exam$predicted_gend <- predict(model_int_gend)

# プロットの作成
ggplot(Exam, aes(x = standLRT, y = normexam, color = schgend)) +
  geom_point(alpha = 0.4) +  # 生徒ごとの散布図
  geom_smooth(method = "lm", se = FALSE, aes(y = predicted_gend)) +  # 予測値の回帰直線
  labs(
    title = "クロスレベル交互作用: standLRT × schgend",
    x = "生徒の以前のテスト得点 (standLRT)",
    y = "標準化された試験得点 (normexam)",
    color = "学校の性別 (schgend)"
  ) +
  theme_minimal()

model_int_schavg: standLRT × schavg の交互作用

# 交互作用モデルの適用
model_int_schavg <- lmer(normexam ~ standLRT * schavg + schgend + (1 | school), data = Exam)

# schavgを3つのカテゴリに分割（低・中・高）
Exam$schavg_cat <- cut(Exam$schavg, breaks = quantile(Exam$schavg, probs = c(0, 0.33, 0.66, 1)),
                       labels = c("低", "中", "高"), include.lowest = TRUE)

# 予測値の取得
Exam$predicted_schavg <- predict(model_int_schavg)

# プロットの作成
ggplot(Exam, aes(x = standLRT, y = normexam, color = schavg_cat)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = "lm", se = FALSE, aes(y = predicted_schavg)) +
  labs(
    title = "クロスレベル交互作用: standLRT × schavg",
    x = "生徒の以前のテスト得点 (standLRT)",
    y = "標準化された試験得点 (normexam)",
    color = "学校の平均学力 (schavg)"
  ) +
  theme_minimal()

model_int_both: standLRT × schgend と standLRT × schavg の両方の交互作用を含むモデル

# 予測値の取得
Exam$predicted_both <- predict(model_int_both)

# プロットの作成（schgendによる色分け + schavgのカテゴリ別でファセット表示）
ggplot(Exam, aes(x = standLRT, y = normexam, color = schgend)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = "lm", se = FALSE, aes(y = predicted_both)) +
  facet_wrap(~ schavg_cat) +  # 学校の平均学力カテゴリで分割表示
  labs(
    title = "クロスレベル交互作用: standLRT × (schgend + schavg)",
    x = "生徒の以前のテスト得点 (standLRT)",
    y = "標準化された試験得点 (normexam)",
    color = "学校の性別 (schgend)"
  ) +
  theme_minimal()

model_int_both_slope: model_int_bothに加えて、standLRT のランダム傾きを含むモデル

# 予測値の取得（ランダムスロープを考慮）
Exam$predicted_both_slope <- predict(model_int_both_slope)

# プロットの作成
ggplot(Exam, aes(x = standLRT, y = normexam, color = schgend)) +
  geom_point(alpha = 0.4) +
  geom_line(aes(y = predicted_both_slope, group = school), alpha = 0.6) +  # 学校ごとの予測線
  facet_wrap(~ schavg_cat) +
  labs(
    title = "クロスレベル交互作用: standLRT × (schgend + schavg) + ランダムスロープ",
    x = "生徒の以前のテスト得点 (standLRT)",
    y = "標準化された試験得点 (normexam)",
    color = "学校の性別 (schgend)"
  ) +
  theme_minimal()

ごちゃごちゃしているので、少し工夫を加える。グループ平均の傾きを強調表示する。

# 平均的な傾きを強調
ggplot(Exam, aes(x = standLRT, y = normexam, color = schgend)) +
  geom_point(alpha = 0.4) +
  stat_smooth(method = "lm", se = FALSE, size = 2) +  # 太い平均線を表示
  facet_wrap(~ schavg_cat) +
  labs(
    title = "学校の性別ごとの平均的な傾き",
    x = "生徒の以前のテスト得点 (standLRT)",
    y = "標準化された試験得点 (normexam)",
    color = "学校の性別 (schgend)"
  ) +
  theme_minimal()

特定の学校だけを選んで表示する

ランダムに選んだ数校（例: 5校）だけの線を表示する。

# ランダムに5校だけ抽出
set.seed(123)
selected_schools <- sample(unique(Exam$school), 5)
subset_data <- Exam[Exam$school %in% selected_schools, ]

# プロット
ggplot(subset_data, aes(x = standLRT, y = normexam, color = factor(school))) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = "lm", se = FALSE) +
  facet_wrap(~ schavg_cat) +
  labs(
    title = "ランダムに選んだ5校の傾き",
    x = "生徒の以前のテスト得点 (standLRT)",
    y = "標準化された試験得点 (normexam)",
    color = "学校ID"
  ) +
  theme_minimal()

デルのランダム効果（ranef）を使い、切片と傾きのばらつきを直接可視化する

library(lattice)

# ランダム効果の可視化
ranef_plot <- ranef(model_int_both_slope, condVar = TRUE)
dotplot(ranef_plot, scales = list(relation = "free"))

交互作用プロットでstandLRT と schgend の効果を比較

interaction.plot(Exam$standLRT, Exam$schgend, Exam$normexam,
                 col = c("red", "green", "blue"), lwd = 2,
                 xlab = "生徒の以前のテスト得点 (standLRT)",
                 ylab = "標準化された試験得点 (normexam)",
                 legend = TRUE)