交互作用を表現する線形モデルをRで実施する
はじめに
前回の記事で交互作用の概念について説明したが、その続きとして、シミュレーションデータを使ってRでの実行例を見ていきたいと思う。
mstour.hatenablog.com
なお今回は2値のカテゴリー変数どうしの交互作用だけを扱うことにする。
問題設定は次の通りである。
目的変数（結果側の変数）として期末テストの得点、説明変数として補習を受けたかどうか（補習あり or 補習なし）と性別（男子学生 or 女子学生）のデータがあり、テストの得点に対する補習の効果を見たいものとする。
補習の効果と性別との交互作用がないケース・あるケースそれぞれのサンプルデータを作成し、線形回帰モデルを実行して結果を見ていくことにしよう。
交互作用のないケース
最初に、補習の効果と性別とには交互作用がないようなデータを作成する。
#----- 前準備 ----- library(tidyverse) library(truncnorm) set.seed(992) n = 300
ここで、テストの得点の生成に「切断正規分布」を利用するためにtruncnormパッケージを読みこんでいる。切断正規分布とは、ランダムサンプルのとりうる値をある範囲に制限した正規分布のことで、以下で0から100までの値しか生成されないように設定している。
#----- 各変数の作成 -----
# ID
ID = 1:n
# Group
Hosyu = c(rep("Y", n/2), rep("N", n/2))
# Hosyu dummy
Hosyu_Y = rep(0, n)
Hosyu_Y[Hosyu == "Y"] = 1
# Gender
Gender = c(rep("M", n/4), rep("F", n/4), rep("M", n/4), rep("F", n/4))
# Male dummy
Gender_Male = rep(0, n)
Gender_Male[Gender == "M"] = 1
# Final exam. score (No interaction)
Score1 = round(rtruncnorm(n, a = 0, b = 100, mean = 30, sd = 10)
+ Hosyu_Y * 20 + Gender_Male * (-10), 0)
#----- データフレーム作成 -----
ad1 = data.frame(ID = ID, Hosyu = Hosyu, Gender = Gender,
Hosyu_Y = Hosyu_Y, Gender_Male = Gender_Male, Score = Score1)補習については補習ありの場合に1・なしの場合に0をとるダミー変数を、性別については男子学生の場合に1・女子学生の場合に0をとるダミー変数をそれぞれ用意している。
テストの得点は、補習の有無と性別の両方から影響を受けるが、両者の交互作用はないようにデータを生成している。
最後に各変数をデータフレームにまとめる。
以下のコードで、補習あり・なしそれぞれのグループのテスト得点の平均値を可視化してみる。
#----- プロット作成用に、補習の有無および性別ごとの平均値を算出 -----
ad1summary = summarise(group_by(ad1, Gender, Hosyu), m = mean(Score))
ad1summary$Hosyu_Y = c(0, 1, 0, 1)
#----- 平均値をプロットする -----
ggplot(ad1summary, aes(x = Hosyu_Y, y = m, color = Gender)) +
geom_line() + geom_point() +
scale_x_continuous(breaks = c(0, 1), labels = c("Nashi", "Ari")) +
xlab("Hosyu") + ylab("Final exam score")通常、2値の変数を使ってこのような折れ線グラフを描くことはあまりないと思うが、前回の話の流れと繋がりやすくするためにあえて作成した。横軸の「補習あり・なし」にはダミー変数を使っているが、ラベルを細工してカテゴリー型変数っぽくしている。
プロットは以下のようになる。
緑色の直線の傾きが男子学生の補習の効果、赤色の直線の傾きが女子学生の補習の効果を示している。直線はほとんど平行であり、補習の効果と性別とに交互作用はなさそうなことが分かる。
さて、まずはこのデータに「交互作用項を含めた線形回帰モデル」（補習と性別とに交互作用があると仮定したモデル）を当てはめることにする。
lm1_int = lm(Score ~ Hosyu_Y + Gender_Male + Hosyu_Y*Gender_Male, data = ad1) summary(lm1_int)
結果は以下のようになった。
赤線で示す部分が交互作用項の係数の推定を表しているが、p値（「係数が0である」という帰無仮説に対応する）を見る限り、交互作用があることを積極的に示すような結果にはなっていない*1。
次に、「交互作用項を除いた線形回帰モデル」（交互作用がないと仮定したモデル）を当てはめてみよう。
lm1_noint = lm(Score ~ Hosyu_Y + Gender_Male, data = ad1) summary(lm1_noint)
補習の効果の推定値（Estimate）は、交互作用項を含めたモデルとほぼ同じになった（交互作用項を含んだモデルでは、係数の推定値が0に近かったため）。
交互作用項を含めるべき理由は特に見当たらないため、今回のデータに関しては、交互作用項のないこちらのモデルで十分だろう。
続いて、交互作用のあるようなデータを作ってみよう。
交互作用のあるケース
説明変数については先ほど作成したものを流用して、テストの得点のみ、補習と性別とに交互作用が生じるような設定で新たに作成する。
# Final exam. score (Interaction Gender*Hosyu)
Score2 = round(rtruncnorm(n, a = 0, b = 100, mean = 30, sd = 10)
+ Hosyu_Y * 15 + Gender_Male * (-10) + Gender_Male * Hosyu_Y * 10, 0)
#----- データフレーム作成 -----
ad2 = data.frame(ID = ID, Hosyu = Hosyu, Gender = Gender,
Hosyu_Y = Hosyu_Y, Gender_Male = Gender_Male, Score = Score2)先ほどと同様にプロットしてみよう。
#----- プロット作成用に、補習の有無および性別ごとの平均値を算出 -----
ad2summary = summarise(group_by(ad2, Gender, Hosyu), m = mean(Score))
ad2summary$Hosyu_Y = c(0, 1, 0, 1)
#----- 平均値をプロットする -----
ggplot(ad2summary, aes(x = Hosyu_Y, y = m, color = Gender)) +
geom_line() + geom_point() +
scale_x_continuous(breaks = c(0, 1), labels = c("Nashi", "Ari")) +
xlab("Hosyu") + ylab("Final exam score")
直線は平行になっておらず、男子学生と女子学生とで補習の効果が大きく異なっていることが分かる。
交互作用のないシミュレーションデータと同じように、まずは「交互作用項を含めた線形回帰モデル」によるモデリングを行ってみよう。
lm2_int = lm(Score ~ Hosyu_Y + Gender_Male + Hosyu_Y*Gender_Male, data = ad2) summary(lm2_int)
補習と性別との交互作用項に関する推定結果を見ると、「Estimate」が7.96となっている。この項は男子学生と女子学生との補習効果の差を示しており、男子学生は女子学生よりも補習の効果がおよそ8点大きいと推定されている。
また係数についてのp値は非常に小さく、交互作用があることが強く示唆される結果となっている。
同様に「交互作用項を除いた線形回帰モデル」も当てはめてみる。
lm2_noint = lm(Score ~ Hosyu_Y + Gender_Male, data = ad2) summary(lm2_noint)
交互作用項を含むモデルと比較すると、補習の効果が4点程度大きくなっている。
これは、本来は男女で異なっているはずの補習の効果（交互作用ありのモデルからは、男子学生の効果は15.68 + 7.96 = 23.64、女子学生の効果は15.68と推定されている）を、男女で平均した結果となっている（(23.64 + 15.68) / 2 = 19.66）*2。
対象母集団全体での補習の効果を見積もりたい場合には交互作用項のないモデルで問題はないが、データを正確にモデリングすることが目的であれば、このようなケースでは交互作用項を含めたモデルを用いるほうがよいだろう。
まとめ
交互作用項を含めた線形回帰モデルのRでの実行について、交互作用のある場合・ない場合それぞれのシミュレーションデータを用いて説明した。
今回のデータ例において、交互作用項のあるモデルは男女別に回帰モデルを当てはめるのと同じようなことを行っている（厳密には異なる）。
Rでのプロット作成に関しては、以下の資料を参考にした。
[1] グラフ頂上決戦 もしも、SASのsgplotとRのggplot2を比較したら・・・
http://nfunao.web.fc2.com/files/sgplot_vs_ggplot2.pdf
[2] Mukku John Blog：ggplot2を使って、軸を制御する-3
http://mukkujohn.hatenablog.com/entry/2016/10/11/220722
