ダミー変数・交互作用
佐藤圭一
2023-09-12
1 ダミー変数
1.1 例題
ある政治家の好感度を0~100点で測定したとする。政治家の好感度に対する性別の影響について調べたい。ただし年齢を統制変数とする。
模擬データの作成
age <- c(57, 49, 50, 74, 48, 63, 66, 32, 18, 41, 34, 64, 75, 48, 41, 51, 54, 55, 58, 66)
gender <- c(1, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1)
gender <- factor(gender,
levels = c(0, 1),
labels = c("Female", "Male")) # 0=female, 1=maleと定義
like <- c(93, 47, 83, 42, 35, 66, 95, 43, 63, 80, 23, 97, 85, 67, 81, 66, 92, 89, 16, 63)
dat <- data.frame(
ID = 1:20,
age = age, # 年齢
gender = gender, # 性別
like = like # 政治家好感度
)
# データを見てみる(最初の6行分)
head(dat)## ID age gender like
## 1 1 57 Male 93
## 2 2 49 Female 47
## 3 3 50 Male 83
## 4 4 74 Female 42
## 5 5 48 Female 35
## 6 6 63 Male 66分析
- 従属変数：政治家の好感度
- 独立変数：性別（０＝女性；１＝男性）、年齢（統制変数）
# モデルに従って計算をし、outというobjectに付ける
out <- lm(formula = like ~ gender + age, # object名で式を書く
data = dat) # データ元
summary(out)##
## Call:
## lm(formula = like ~ gender + age, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -29.962 -13.184 2.093 11.235 35.974
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 35.5091 14.9047 2.382 0.029144 *
## genderMale 35.6386 8.0142 4.447 0.000354 ***
## age 0.1802 0.2808 0.642 0.529479
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 17.14 on 17 degrees of freedom
## Multiple R-squared: 0.5693, Adjusted R-squared: 0.5186
## F-statistic: 11.23 on 2 and 17 DF, p-value: 0.0007774# 多重共線性の確認
library(car)
vif(out) # いずれもVIF < 2 で多重共線性の心配はない## gender age
## 1.048954 1.048954# 統計表の出力
library(texreg)
screenreg(out,
custom.coef.names = c("切片", "男性（ダミー）", "年齢"))##
## =====================
## Model 1
## ---------------------
## 切片 35.51 *
## (14.90)
## 男性（ダミー） 35.64 ***
## (8.01)
## 年齢 0.18
## (0.28)
## ---------------------
## R^2 0.57
## Adj. R^2 0.52
## Num. obs. 20
## =====================
## *** p < 0.001; ** p < 0.01; * p < 0.05# もしもhtmlファイルに出力したいならば
library(texreg)
htmlreg(out,
custom.coef.names = c("切片", "男性（ダミー）", "年齢"),
file = "dummy.html")1.2 カテゴリー変数の基準値の変更
基準値を女性から男性に変更してみる。
# 基準値を変えた新しいダミー変数を作成する
dat$gender2 <- dat$gender # gender変数をコピーし新たにgender2を作成
dat$gender2 <- relevel(dat$gender2, ref = "Male")
# きちんとできているか必ず確認する
table(dat$gender) # もとの変数##
## Female Male
## 8 12table(dat$gender2) # 新しい変数（男女の数は変わらないが順番が逆になる）##
## Male Female
## 12 8回帰分析
out2 <- lm(formula = like ~ gender2 + age, # 男性を基準カテゴリーにした変数gender2をgenderの代わりに入れる
data = dat)
summary(out2) # Femaleの効果が表に出る。##
## Call:
## lm(formula = like ~ gender2 + age, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -29.962 -13.184 2.093 11.235 35.974
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 71.1477 16.1261 4.412 0.000381 ***
## gender2Female -35.6386 8.0142 -4.447 0.000354 ***
## age 0.1802 0.2808 0.642 0.529479
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 17.14 on 17 degrees of freedom
## Multiple R-squared: 0.5693, Adjusted R-squared: 0.5186
## F-statistic: 11.23 on 2 and 17 DF, p-value: 0.0007774# genderとgender2をそれぞれ入れた結果を並べてみる
library(texreg)
screenreg(list(out, out2)) # 切片も変わっていることに注意##
## =====================================
## Model 1 Model 2
## -------------------------------------
## (Intercept) 35.51 * 71.15 ***
## (14.90) (16.13)
## genderMale 35.64 ***
## (8.01)
## age 0.18 0.18
## (0.28) (0.28)
## gender2Female -35.64 ***
## (8.01)
## -------------------------------------
## R^2 0.57 0.57
## Adj. R^2 0.52 0.52
## Num. obs. 20 20
## =====================================
## *** p < 0.001; ** p < 0.01; * p < 0.051.3 三つ以上のカテゴリーからなる変数のダミー変数化
模擬データの作成
# 学歴（中卒5人、高卒10人、大卒5人の仮想変数を作成する）
educ <- c(rep("Junior", 5), rep("High", 10), rep("Univ", 5))
educ <- factor(educ) # 因子型として定義するダミー変数の作成
# ライブラリを読み込む
library(makedummies)
# ダミー変数を作成
educ.dat <- makedummies(data.frame(educ), # このパッケージはデータフレーム型しか受け付けないため、データフレームとして括る
basal_level = TRUE) # 基準変数として因子型で定義されているものもアウトプットに含めるか
# 中身を見てみる
head(educ.dat) # それぞれの学歴を二値で置いた新たな変数が作成される。## educ_High educ_Junior educ_Univ
## 1 0 1 0
## 2 0 1 0
## 3 0 1 0
## 4 0 1 0
## 5 0 1 0
## 6 1 0 0回帰分析
# ダミー変数化したデータを既存のデータに張り付ける
dat <- data.frame(dat, educ.dat)
# 中卒が基準カテゴリーとする場合⇒高卒・大卒ダミーのみ投入する。
out3 <- lm(like ~ gender + age + educ_High + educ_Univ,
data = dat)
summary(out3)##
## Call:
## lm(formula = like ~ gender + age + educ_High + educ_Univ, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -22.920 -9.189 -1.072 8.437 29.919
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 29.8719 17.5311 1.704 0.109013
## genderMale 37.1360 8.3252 4.461 0.000458 ***
## age 0.2747 0.2916 0.942 0.361011
## educ_High 4.6056 9.7910 0.470 0.644841
## educ_Univ -9.9840 11.2444 -0.888 0.388603
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 17.01 on 15 degrees of freedom
## Multiple R-squared: 0.6258, Adjusted R-squared: 0.526
## F-statistic: 6.272 on 4 and 15 DF, p-value: 0.0035772 交互作用
2.1 例題
ある政治家の好感度を0~100点で測定したとする。政治家の好感度に対して、性別と年齢の交互作用について調べたい。
データの読み込み
age <- c(57, 49, 50, 74, 48, 63, 66, 32, 18, 41, 34, 64, 75, 48, 41, 51, 54, 55, 58, 66)
gender <- c(1, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1)
gender <- factor(gender,
levels = c(0, 1),
labels = c("Female", "Male")) # 0=female, 1=maleと定義
tv <- c(2, 2, 2, 2, 3, 3, 2, 3, 2, 2, 4, 2, 1, 3, 2, 3, 2, 2, 3, 3)
like <- c(93, 47, 83, 42, 35, 66, 95, 43, 63, 80, 23, 97, 85, 67, 81, 66, 92, 89, 16, 63)
# データフレームにまとめる
dat2 <- data.frame(
ID = 1:20,
age = age, # 年齢
gender = gender, # 性別
tv = tv, # テレビの視聴頻度
like = like # 政治家好感度
)
# データを見てみる(最初の6行分)
head(dat2)## ID age gender tv like
## 1 1 57 Male 2 93
## 2 2 49 Female 2 47
## 3 3 50 Male 2 83
## 4 4 74 Female 2 42
## 5 5 48 Female 3 35
## 6 6 63 Male 3 66分析
- 従属変数：政治家の好感度
- 独立変数：年齢 × テレビの視聴頻度（交互作用）、なお年齢、テレビの視聴頻度の主効果についても含める
- 統制変数：性別
# 交互作用用のパッケージを読み込む
library(pequod)
# モデルを組む（テレビの年齢、視聴頻度をそれぞれ中心化する）
out.intr <- lmres(like ~ age + tv + gender + age*tv,
centered = c("age", "tv"),
data = dat2)
# 結果
summary(out.intr)## Formula:
## like ~ age + tv + genderMale + age.XX.tv
## <environment: 0x000002c6860e0498>
##
## Models
## R R^2 Adj. R^2 F df1 df2 p.value
## Model 0.969 0.940 0.923 58.301 4.000 15 5.8e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residuals
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -10.2359 -4.3669 0.7971 0.0000 4.1978 10.3843
##
## Coefficients
## Estimate StdErr t.value beta p.value
## (Intercept) 41.34991 2.87138 14.40071 <2e-16 ***
## age -0.26730 0.12460 -2.14534 -0.1552 0.0487 *
## tv -22.61931 2.43936 -9.27263 -0.6230 <2e-16 ***
## genderMale 39.11191 3.61638 10.81522 0.7956 <2e-16 ***
## age.XX.tv -0.50612 0.16390 -3.08802 -0.2274 0.0075 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Collinearity
## VIF Tolerance
## age 1.2992 0.7697
## tv 1.1204 0.8925
## genderMale 1.3432 0.7445
## age.XX.tv 1.3461 0.7429交互作用のプロット
# スロープの計算
slope <- simpleSlope(out.intr, pred = "age", mod1 = "tv") # テレビ視聴量が±1SDごとにスロープを計算
summary(slope) # 計算結果をみる##
## ** Estimated points of like **
##
## Low age (-1 SD) High age (+1 SD)
## Low tv (-1 SD) 55.637 57.850
## High tv (+1 SD) 34.733 17.179
##
##
##
## ** Simple Slopes analysis ( df= 15 ) **
##
## simple slope standard error t-value p.value
## Low tv (-1 SD) 0.0771 0.1364 0.57 0.5800
## High tv (+1 SD) -0.6118 0.1932 -3.17 0.0064 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
##
## ** Bauer & Curran 95% CI **
##
## lower CI upper CI
## tv -1.3666 -0.0565# スロープのプロット
PlotSlope(slope,
namemod = c("TV時間短", "TV時間長"),
namex = "年齢", # X軸名
namey = "好感度") # y軸名
2.2 カテゴリー変数の交互作用
カテゴリー変数（ここでは性別）を調整変数とした場合の交互作用の扱いについて取り上げる。以下のようにモデルを組むとする。 - 従属変数：好感度 - 独立変数：年齢×性別（調整変数）、ただし、年齢（主効果）、性別（主効果）もふくめる。
# 年齢を中心化する
dat2$age_c <- dat2$age - mean(dat2$age)
# モデルを組む
out.intr2 <- lm(like ~ age_c + gender + age_c*gender, # カテゴリー変数genderを含むので、通常のlm関数の方が後の処理がしやすい
data = dat2)
summary(out.intr2) # age×genderは有意な効果はない##
## Call:
## lm(formula = like ~ age_c + gender + age_c * gender, data = dat2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -29.970 -13.199 2.082 11.235 35.953
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 44.919783 6.366403 7.056 2.72e-06 ***
## age_c 0.181022 0.330665 0.547 0.591626
## genderMale 35.642096 8.290425 4.299 0.000552 ***
## age_c:genderMale -0.003406 0.683504 -0.005 0.996086
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 17.67 on 16 degrees of freedom
## Multiple R-squared: 0.5693, Adjusted R-squared: 0.4885
## F-statistic: 7.049 on 3 and 16 DF, p-value: 0.003105# 多重共線性の確認
library(car)
vif(out.intr2, type = "predictor")## GVIFs computed for predictors## GVIF Df GVIF^(1/(2*Df)) Interacts With Other Predictors
## age_c 1 3 1 gender --
## gender 1 3 1 age_c --交互作用のプロット 上述の結果のように、有意な効果はないがプロットの例として示す。
# パッケージの読み込み
library(rockchalk) # prequodパッケージは連続変数にしか使えないため
# プロット
# 以下のプロットでは、実際男性と女性の直線はほぼ並行であることが見える。
plotSlopes(out.intr2,
plotx = "age_c", modx = "gender",
col = c("black", "gray")) 
2.3 交互作用のプロット：より一般的かつ柔軟な方法
ここまでの方法はパッケージの機能に頼った方法で、一般性には欠ける。調整変数がどのような型かに関わらずより一般的にやろうとするならば、次のような方法もある。まずは、例題と同じように、次のようなモデルを組むとする。
- 従属変数：政治家の好感度
- 独立変数：年齢 × テレビの視聴頻度（交互作用） 年齢、テレビの視聴頻度の主効果についても含める
- 統制変数：性別
lm関数を用いた交互作用の計算
# 交互作用に含める変数を中心化する
dat2$tv_c <- dat2$tv - mean(dat2$tv)
dat2$age_c <- dat2$age - mean(dat2$age)
# モデルを組む（通常のlm関数を使う。）
out.intr3 <- lm(like ~ age_c + tv_c + gender + age_c*tv_c,
data = dat2)
summary(out.intr3) このように、交互作用を統計的に計算するだけならば、パッケージに頼らずとも、baseパッケージにある関数（lm）でそのまま計算できる。
交互作用のプロット
プロットをする場合には、effectsパッケージを用いて、回帰式から予測される各予測値を計算し、それをもとにプロットを行う。
library(effects)
# TVの頻度の-1SD, +SD それぞれの値を計算する
tv.L <- mean(dat2$tv_c) - sd(dat2$tv_c)
tv.H <- mean(dat2$tv_c) + sd(dat2$tv_c)
age.L <- mean(dat2$age_c) - sd(dat2$age_c)
age.H <- mean(dat2$age_c) + sd(dat2$age_c)
tv.L <- round(tv.L, 1) #小数第二位以下を丸め、定点にする
tv.H <- round(tv.H, 1)
age.L <- round(age.L, 1)
age.H <- round(age.H, 1)
# モデルに基づきそれぞれの独立変数の値での従属変数の値を推定、
out.intr3.effect <- effect(term = "age_c*tv_c",
mod = out.intr3,
xlevels = list(tv_c = c(tv.L, tv.H),
age_c= c(age.L, age.H)
))
out.intr3.effect ###
## age_c*tv_c effect
## tv_c
## age_c -0.7 0.7
## -14.3 79.40673 57.87229
## 14.3 81.89442 40.09478それぞれのポイントでの予測される好感度の値が記載されている。これを用いて線を引く。
# データ・フレームに変換する
effect.dat <- as.data.frame(out.intr3.effect)
# テレビの視聴度合い（1～2行目；３～４行目）ごとにのポイントを結ぶ直線を引き、
# 図を重ね書きする
plot(x = effect.dat$age_c[1:2], y = effect.dat$fit[1:2], type = "l",
xlim = c(min(effect.dat$age_c), max(effect.dat$age_c)),
ylim = c(min(effect.dat$fit), max(effect.dat$fit)),
lty = 1, # 線の種類
ann = F) # 座標名を描かない
par(new = T) # 重ね書き
plot(x = effect.dat$age_c[3:4], y = effect.dat$fit[3:4], type = "l",
xlim = c(min(effect.dat$age_c), max(effect.dat$age_c)),
ylim = c(min(effect.dat$fit), max(effect.dat$fit)),
lty = 2,
xlab = "年齢",
ylab = "好感度")
legend("right",
legend = c(c("TV時間短", # tv_c = -0.7の方なので
"TV時間長")), # tv_c = 0.7の方なので
lty = c(1,2) # 線種
) 
なお、simpleSlopeでは第三変数の代入値が異なるらしく、上のものと予測値が同じにならないことに注意。
(線の傾きは同じ) #
effect関数で具体的にどのような値が代入されているかは、out.intr3.effectオブジェクト内のmodel、model.matrixをみると分かる。
統計表の出力
library(texreg)
screenreg(out.intr3,
custom.coef.names = c("切片", "年齢", "TV時間","男性（ダミー）","年齢 × TV時間"))##
## =====================
## Model 1
## ---------------------
## 切片 41.35 ***
## (2.87)
## 年齢 -0.27 *
## (0.12)
## TV時間 -22.62 ***
## (2.44)
## 男性（ダミー） 39.11 ***
## (3.62)
## 年齢 × TV時間 -0.51 **
## (0.16)
## ---------------------
## R^2 0.94
## Adj. R^2 0.92
## Num. obs. 20
## =====================
## *** p < 0.001; ** p < 0.01; * p < 0.05# 中心化なし
out.no.centering <- lm(like ~ age + tv + gender, data = dat2)
# 中心化あり
out.with.centering <- lm(like ~ age_c + tv_c + gender, data = dat2)
# 結果の比較
library(texreg)
screenreg(list(
no.centering = out.no.centering,
with.centering = out.with.centering))##
## =========================================
## no.centering with.centering
## -----------------------------------------
## (Intercept) 106.13 *** 45.95 ***
## (12.13) (3.04)
## age -0.14
## (0.15)
## tv -22.08 ***
## (3.01)
## genderMale 33.92 *** 33.92 ***
## (3.96) (3.96)
## age_c -0.14
## (0.15)
## tv_c -22.08 ***
## (3.01)
## -----------------------------------------
## R^2 0.90 0.90
## Adj. R^2 0.88 0.88
## Num. obs. 20 20
## =========================================
## *** p < 0.001; ** p < 0.01; * p < 0.05# 注意：交互作用を入れた場合には、交互作用項が他の変数との強い相関を持つため、この相関に上記のように中心化ありとなしの係数が同じにはならない。