読者です 読者をやめる 読者になる 読者になる

Rを通じて統計学を学ぶ備忘録ブログ

SPSSからRに移行したい私のような人向けのR解説ブログ兼学習用備忘録。

補足の回 重回帰分析の交互作用の検討 R実践

R 多変量解析

【お詫びのお知らせ】20160628

作図2:オンラインツール(Preacher, Curran, & Bauer, 2006)の箇所で入力に誤りがありましたので訂正いたしました。

 

今回は重回帰分析の交互作用の検討をRで実践してみます。前回も同様の記事を書きましたが、今回は地味に地味に色んなことをチェックしながら交互作用をみていきます。パッケージpequod、関数「lmres」はつかいません。そっちの方が簡単に結果は得られますが、過程を追っていくので今回はパス。重回帰分析で交互作用を検討する意味は、「zがどのような値のときに、xの効果がどのようになるかを検討する」ということになると思いますが、そもそも交互作用って??という方は、前回の記事を見て頂ければ幸いです。 

□変数の中心化(センタリング)

交互作用の分析では、まず交互作用を変数を掛け合わせることで作成します。しかし、掛け合わせただけでは、変数同士の相関係数が高くなってしまい、その結果「マルチコリニアリティ(multicollinearity)」が発生してしまいます。多重共線性、通称マルチコです。そのため、相関係数を下げるために、変数の平均から偏差をとる作業、つまり「中心化(センタリング)」を行う必要があるのです。中心化を行うことで、変数の平均値を0にできます。それにより、変数の効果を正しく評価できます。また、中心化した交互作用項、中心化していない交互作用項、どちらを使っても、回帰分析の係数は変化しません。

 f:id:teruaki-sugiura:20150927142720p:plain

demo <- read.csv("inter.csv",header=TRUE)
demo$xz <- demo$x*demo$z
demo$centered_x <- demo$x-(mean(demo$x))
demo$centered_z <- demo$z-mean(demo$z)
demo$centered_xz <- (demo$centered_x)*(demo$centered_z)

f:id:teruaki-sugiura:20150927134750p:plain

summary(demo)

中心化(センタリング)を行った、変数centerd_xとcenterd_zの平均値は0になっていることが確認できます。

f:id:teruaki-sugiura:20150927140402p:plain

round(cor(demo), 2)

また、相関係数を確認してみると、中心化(センタリング)を行っていない変数xzとx, zでは0.58, 0.74となっており、相関係数が高い状態です。一方、中心化(センタリング)を行った変数数centerd_x:centerd_zと数centerd_x, centerd_zの相関係数は-0.07, 0.24となっており、先ほどの数値と比べても下がっていることが確認できます。

あくまで、中心化(センタリング)は交互作用項xzと変数x, zの相関を下げるためのもので、変数x, zの相関係数が下がるわけではありません。f:id:teruaki-sugiura:20150927135511p:plain

result_normal <- lm(y ~ x + z + x:z, data = demo)
result_centered <- lm(y ~ centered_x + centered_z + centered_x:centered_z, data = demo)
summary(result_normal)
summary(result_centered)

f:id:teruaki-sugiura:20150927140909p:plain

中心化した交互作用項(result_centered)、中心化していない交互作用項(result_normal)の回帰係数を確認してみると、どちらの場合でも、回帰係数は変化していません。

中心化(センタリング)ではなく、標準化(平均0、標準偏差1)にした場合も、後々、xの主効果を検討する際に用いる単純傾斜分析の有意性は変化しません。しかし、標準化した値というのは、平均0、標準偏差1に規格化されているので、実際の観測値を計算するのに手間がかかります。

□単純傾斜分析(Simple Slope Analysis)

交互作用項が有意な場合、単純傾斜分析を行います。単純傾斜分析は変数xの主効果が有意なのかどうかを検討するためのものです。方法は、群分けする変数の±1sdで、xの傾きが0と異なるかを検定します。手順は、z高群(z低群)の変数を作成し、xとz高群(xとz低群)の交互作用項xz_high(xz_low)を作成し、重回帰分析を行い、xの回帰係数が有意かどうかをみることになります。
+1sd変数は(平均+1sd)、-1sd変数は(平均-1sd)することになります。
demo$z_high <- demo$z - (mean(demo$z) + sd(demo$z))
demo$z_low <- demo$z - (mean(demo$z) - sd(demo$z))
result_high <- lm(y ~ centered_x + z_high + centered_x:z_high, data = demo)
result_low <- lm(y ~ centered_x + z_low + centered_x:z_low, data = demo)
summary(result_high)
summary(result_low)

f:id:teruaki-sugiura:20150927141227p:plain

---------------------------------------------------------------------------------------

f:id:teruaki-sugiura:20150927141236p:plain

上記の単純傾斜分析の結果をみる限り、z低群ではxの効果はないが、z高群ではxの効果があることが確認できます。あとは、交互作用という数値だけでは理解しにくいので、図を作成して結果を解釈していきましょう。

□作図1

以下のコードは関西学院清水先生のものをベースに少しだけいじりました。下のコードを自分の分析に合わせて名前を変えていき、読み込めば作図完了です。

x_sd <- sd(demo$x)
y_sd <- sd(demo$y)
y_mean <- mean(demo$y)
x_name <- "x"
y_name <- "y"
z_name <- "z"
centered_x_name <-"centered_x"
y11 <- coef(result_low)[1] - coef(result_low)[centered_x_name] * x_sd
y12 <- coef(result_low)[1] + coef(result_low)[centered_x_name] * x_sd
y21 <- coef(result_high)[1] - coef(result_high)[centered_x_name] * x_sd
y22 <- coef(result_high)[1] + coef(result_high)[centered_x_name] * x_sd
yy1 <- y_mean - y_sd
yy2 <- y_mean + y_sd
xx <- c(-x_sd*1.5,x_sd*1.5)
yy <- c(yy1,yy2)
x_domain <- c(-x_sd,x_sd)
y1 <- c(y11,y12)
y2 <- c(y21,y22)
plot(xx, yy, type = 'n', font = 2, font.lab = 2, xlab = x_name, ylab = y_name, main = 'Interaction_Plot', las = 1, tcl = 0.3)
lines(x_domain, y1, lwd = 3, lty = 1, col = 1)
lines(x_domain, y2, lwd = 3, lty = 6, col = 2)
points(x_domain, y1, col =1, pch = 15)
points(x_domain, y2, col = 2, pch = 15)
legend("bottomleft", legend = c(paste(z_name, "+1SD"), paste(z_name, "-1SD")), lwd = c(3, 3), lty = c(6, 1), col = c(2, 1), bty = "n")

f:id:teruaki-sugiura:20160628020218p:plain

□作図 2

オンラインツール(Preacher, Curran, & Bauer, 2006)のSimple intercepts, simple slopes, and regions of significance in MLR 2-way interactionsを使用すれば、一瞬で作図ができますが、共分散などが必要になりますが、vcov関数で入手します。

f:id:teruaki-sugiura:20160628015634p:plain

f:id:teruaki-sugiura:20150927152550p:plain 

f:id:teruaki-sugiura:20160628015710p:plain

※この表の入力に誤りがありました。

※zが0, 1の2値変数であれば、チェックボックスにチェックを入れ、cv_z1に0、cv_z3に1を入力してください。 

xx <- c(-1.5,1.5) # <-- change to alter plot dims
yy <- c(2.3267,3.8061) # <-- change to alter plot dims
leg <- c(-1.5,2.5116) # <-- change to alter legend location
x <- c(-1.5,1.5) # <-- x-coords for lines
y1 <- c(2.774,3.8061)
y2 <- c(2.6922,3.1897)
y3 <- c(2.6105,2.5733)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='X',ylab='Y',main='MLR 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('CVz1(1)','CVz1(2)','CVz1(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))

f:id:teruaki-sugiura:20160628020111p:plain

確認用

TWO-WAY INTERACTION SIMPLE SLOPES OUTPUT

Your Input
=======================================================
X1 = -1.5
X2 = 1.5
cv1 = 1.843251
cv2 = 0
cv3 = -1.843251
Intercept = 2.94096
X Slope = 0.16582
Z Slope = 0.18939
XZ Slope = 0.09669
df = 96
alpha = 0.05 Asymptotic (Co)variances
=======================================================
var(b0) 0.00565453
var(b1) 0.00526563
var(b2) 0.00177854
var(b3) 0.00161122
cov(b2,b0) -0.00004498
cov(b3,b1) 0.00016283 Region of Significance
=======================================================
Z at lower bound of region = -10.0355
Z at upper bound of region = -0.2241
(simple slopes are significant *outside* this region.) Simple Intercepts and Slopes at Conditional Values of Z
=======================================================
At Z = cv1...
simple intercept = 3.2901(0.1074), t=30.6381, p=0
simple slope = 0.344(0.1065), t=3.2308, p=0.0017
At Z = cv2...
simple intercept = 2.941(0.0752), t=39.1103, p=0
simple slope = 0.1658(0.0726), t=2.2851, p=0.0245
At Z = cv3...
simple intercept = 2.5919(0.1089), t=23.7966, p=0
simple slope = -0.0124(0.1007), t=-0.1232, p=0.9022 Simple Intercepts and Slopes at Region Boundaries
=======================================================
Lower Bound...
simple intercept = 1.0403(0.4309), t=2.4143, p=0.0177
simple slope = -0.8045(0.4053), t=-1.985, p=0.05
Upper Bound...
simple intercept = 2.8985(0.0759), t=38.178, p=0
simple slope = 0.1442(0.0726), t=1.985, p=0.05 Points to Plot
=======================================================
Line for cv1: From {X=-1.5, Y=2.774} to {X=1.5, Y=3.8061}
Line for cv2: From {X=-1.5, Y=2.6922} to {X=1.5, Y=3.1897}
Line for cv3: From {X=-1.5, Y=2.6105} to {X=1.5, Y=2.5733}

以上で今回はおしまい。

ーーーーーーーーーーーーーーーーーーーーー

□交互作用が連続変数とカテゴリ変数の場合

fit<- lm(y ~ class * x, data =data)
## 予測値のデータフレームを作る
pred <- with(data, expand.grid(x = seq(min(x, na.rm = TRUE), max(x, na.rm = TRUE), length =100),class = c(0,1)))
pred <- mutate(pred, y = predict(fit, newdata = pred))
err <- predict(fit, newdata = pred, se.fit = TRUE)
pred$lower <- err$fit + qt(0.025, df = err$df) * err$se.fit
pred$upper <- err$fit + qt(0.975, df = err$df) * err$se.fit
p4 <- ggplot(data, aes(x = x, y = y, col = as.factor(class))) + geom_point() + geom_line(data = pred)
p4 <- p4 + geom_smooth(data = pred, aes(ymin = lower, ymax = upper), stat ="identity")
p4 <- p4 + scale_color_discrete(name = "class", labels = c("0","1")) + guides(color = guide_legend(reverse = TRUE))
p4 <- p4 + ggtitle("interaction x with class") + labs(x = "x", y = "y")
print(p4)

f:id:teruaki-sugiura:20160111214623p:plain

パッケージpequodを使う場合

パッケージpequodのインストールは忘れずに。関数は、「lmres」を使います。引数は以下です。

lmres(目的変数 ~ 説明変数 * 調整変数, centered = c(“目的変数”, “調整変数”), data = データ)

回帰式が通常の重回帰分析を含む書き方と同じです。centerdでは中心化する変数を指定します。

demo <- read.csv("inter.csv")
head(demo)

f:id:teruaki-sugiura:20160628033820p:plain

lm.inter <- lmres(y ~ x*z, centered = c("x", "z"), data = demo)
summary(lm.inter)

f:id:teruaki-sugiura:20160628033842p:plain

lm.slope<- simpleSlope(lm.inter, pred = "x", mod1 = "z")
summary(lm.slope)

f:id:teruaki-sugiura:20160628033901p:plain

PlotSlope(lm.slope)

f:id:teruaki-sugiura:20160628033603p:plain

 

zが高群(+1SD)の場合に、xの効果は有意です。

しかし、zが低群(-1SD)の場合、非有意であることがわかります。

 

参考になるリンク集

平川真(2014)重回帰分析で交互作用効果

深谷達史(2012)重回帰分析による調整効果(交互作用)の検討

前田和寛(2008)重回帰分析の応用的手法-交互作用項ならびに統制変数を含む分析-

 

 

広告を非表示にする