1 この資料の目的

この資料は，重回帰分析の単純傾斜を検討するときにする変数処理の理屈について，個人的な理解を，統計ソフトのRを使ってデータと照らし合わせながらまとめたものです。以下，なんでこんなものを作ったのかの理由です。

重回帰分析で交互作用の影響を検討するとき，交互作用項を作って重回帰モデルに入れると思います。多くの場合，交互作用項の影響が有意であれば，次に下位検定である単純傾斜（simple slope）の影響を検討をしていきます。

単純主効果を検討するときに，量的変数の調整変数が高い時と低い時の変数を作って，それぞれ検討します。その際，もとの値から-1SDした調整変数の値を調整変数が高い時，もとの値から+1SDした調整変数の値を調整変数が低い時とするのが通例と思います（Rを使った具体的な実行方法だと，例えば広島大学の平川さんの資料がシンプルにまとまっていてわかりやすいかと思います）。

このとき，素朴に感じる疑問は，なんで調整変数の高い時を検討するのに-1SDした変数を高い時の値として扱うのか，ということです。今までグラフの平行移動をイメージして，調整変数の高い場合を検討したいときにはプラス側に移動するから-1SD(実際?は-(+1SD))して，調整変数の低い場合を検討したいときにはマイナス側に移動するから+1SD(実際?は-(-1SD))するという理解でした。でもこれが何でなのか，友人に聞かれた際に明確に答えられなかったので，自分の理解の仕方をまとめてみようと思って作ってみた次第です。

私はあまり数学的な素養がないので，実際のデータと照らし合わせながら理解することを目指します。以下では，Rを使ってやっている処理を可視化しながら説明をしてみます。

library(knitr)
library(rmdformats)
library(tidyverse)
library(MASS)
library(lattice)
library(jtools)

## Global options
options(max.print="75")
opts_knit$set(width=75)

2 使うデータについて

データは関西学院大学の清水先生のブログからもってきました。詳細も書かれています。以下の3つの変数を関係を検討します。

• $y$(従属変数): 満足度($sat$)

• $x$(説明変数): 話す頻度($talk$)

• $z$(調整変数): パフォーマンス($per$)

dat01 <- read_csv("sampledata_reg.csv")

## Parsed with column specification:
## cols(
##   group = col_double(),
##   sat = col_double(),
##   talk = col_double(),
##   per = col_double(),
##   con = col_double()
## )

dat01としてサンプルデータを読み込みます。データフレームの上からいくつかのデータを見てみるとこんな感じです。

kable(head(dat01))

とりあえず，可視化してみます。

g <- ggplot(dat01, aes(x = talk, y = sat, group = per, colour = per)) +
  geom_jitter(height = 0.1, width = 0.4)
g

見た感じ，色の濃いのが右肩下がりで色の薄いのが右肩上がりな様子。 もう少しわかりやすくするために，次は，パフォーマンス(per)の高群と低群のみのデータで描いてみる。なお，今回は高群を1, 低群を0, 中間を2とコーディングしています。

per_h <- mean(dat01$per) + sd(dat01$per)
per_l <- mean(dat01$per) - sd(dat01$per)

dat01$per_hl <- ifelse(dat01$per >= per_h, 1, 
                   ifelse(dat01$per <= per_l, 0, 2))

群分けしたときのそれぞれの人数の確認をしてみます。Var1っていうのが群分けの変数です。Freqは人数です。

kable(table(dat01$per_hl)) 

高低群のみをグラフにして，それぞれの群の回帰直線をかいてやるとわかりやすいかもしれません。

dat02 <- dplyr::filter(dat01, per_hl != 2)
dat02$per_hl <- as.factor(dat02$per_hl) 
g01 <- ggplot(dat02, aes(x = talk, y = sat, group = per_hl, colour = per_hl)) +
  geom_jitter(height = 0.1, width = 0.4) +
  stat_smooth(method = "lm")
g01

交互作用有りっぽい図が描けました。こういう結果(高群は回帰係数が正，低群は回帰係数が0か負)に近い分析結果ならokだと考えられます。

4 交互作用項を含む重回帰分析における単純主効果の検討方法の意味

sd(dat01$per)

## [1] 1.756184

dat01m <- dplyr::filter(dat01, per_hl == 2)
g03m <- ggplot(dat01m, aes(x = talk_c, y = sat, colour = per_hl)) +
  geom_jitter(height = 0.1, width = 0.4) +
  stat_smooth(method = "lm")
g03m

lattice::cloud(sat ~ talk_c * per_c, data = dat01, groups = per_hl)

lattice::cloud(sat ~ talk * per, data = dat02, groups = per_hl)

dat02 <- mutate(dat02, talk_c = talk - mean(talk))
g01 <- ggplot(dat02, aes(x = talk_c, y = sat, group = per_hl, colour = per_hl)) +
  geom_jitter(height = 0.1, width = 0.4) +
  stat_smooth(method = "lm")
g01

lattice::cloud(sat ~ talk * per, data = dat01, groups = per_hl)

g04 <- ggplot(dat01, aes(x = per_c, y = sat)) +
  geom_jitter(height = 0.1, width = 0.4, colour = "black", alpha = 0.5) +
  stat_smooth(method = "lm", colour = "black") + geom_jitter(aes(x = per_high, y = sat), height = 0.1, width = 0.4, colour = "red", alpha = 0.5) +
  stat_smooth(aes(x = per_high, y = sat), method = "lm", colour = "red") + 
  geom_jitter(aes(x = per_low, y = sat), height = 0.1, width = 0.4, colour = "blue", alpha = 0.5) +
  stat_smooth(aes(x = per_low, y = sat), method = "lm", colour = "blue") +
  xlab("per_high(red), per_c(black), per_low(blue)")
g04

dat01 %>% 
  dplyr::select(per, per_c, per_high, per_low) %>% 
  summarise_all(funs(mean)) %>% 
  kable(digits =2)

g05 <- ggplot(dat01, aes(x = per_high, y = sat)) +
  geom_jitter(height = 0.1, width = 0.4, colour = "red", alpha = 0.5) +
  stat_smooth(method = "lm", colour = "red") +
  geom_vline(xintercept = 0)
g05

g06 <- ggplot(dat01, aes(x = per_low, y = sat)) +
  geom_jitter(height = 0.1, width = 0.4, colour = "blue", alpha = 0.5) +
  stat_smooth(method = "lm", colour = "blue") +
  geom_vline(xintercept = 0)
g06

func <- function(x){2*x}
func_h <- function(x){2*x - 4}
func_l <- function(x){2*x + 4}
plot(func, -5, 5, col = "green")
plot(func_h, -5, 5, col = "blue", add = TRUE)
plot(func_l, -5, 5, col = "red", add = TRUE)
abline(v = 0, lty = 2, add =TRUE)
abline(h = 0, lty = 2, add =TRUE)
text(4, labels = "y = 2x + 4", col = "red")
text(0, labels = "y = 2x", col = "green")
text(-5, labels = "y = 2x - 4", col = "blue")

x = seq(-5, 5, 1)
y = seq(-10, 10, 2)
d <- matrix(c(x, y), ncol = 2)
dat_test <- data.frame(d) 
colnames(dat_test) <- c("x", "y")
dat_test$x_plus_1sd <- dat_test$x + 2
g <- ggplot(dat_test) +
  geom_point(aes(x = x, y = y), colour = "green") +
  geom_point(aes(x = x_plus_1sd, y = y), colour = "blue")
g

5 まとめ

調整変数の高群，低群を作って再度重回帰分析で説明変数の変化を検討する場合，

• 高群にはセンタリングした調整変数に$-1SD$した変数を作成して，重回帰分析を行う。

• 低群にはセンタリングした調整変数に$+1SD$した変数を作成して，重回帰分析を行う。

• 符号が逆転する理由は，説明変数の回帰係数の解釈のさせ方を変えられるようにするためである。 *ie., 高群であれば，調整変数の全体の平均値を1SD低くしているため，調整変数が0のときの説明変数の回帰係数の解釈は，調整変数が1SD分の高いときの解釈となる。低群であれば，調整変数の全体の平均値を1SD高くしているため，調整変数が0のときの説明変数の回帰係数の解釈は，調整変数が1SD分の低いときの解釈となる。

ついでに，

• 重回帰式に代入して，説明変数の回帰係数を算出するとき：センタリングした変数を用いて推定された回帰係数を代入し，調整変数の得点には，高群のときには+1SDの値を，低群のときには-1SDの値を代入したら算出できる。ただし，それにくわえてプールされた標準誤差を算出して，回帰係数の差分の有意性検定を行うことになるっぽい。

6 実行環境

sessionInfo()

## R version 3.5.2 (2018-12-20)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18362)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=Japanese_Japan.932  LC_CTYPE=Japanese_Japan.932   
## [3] LC_MONETARY=Japanese_Japan.932 LC_NUMERIC=C                  
## [5] LC_TIME=Japanese_Japan.932    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] stargazer_5.2.2  jtools_2.0.1     lattice_0.20-38  MASS_7.3-51.1   
##  [5] forcats_0.4.0    stringr_1.4.0    dplyr_0.8.0.1    purrr_0.3.0     
##  [9] readr_1.3.1      tidyr_0.8.2      tibble_2.0.1     ggplot2_3.1.0   
## [13] tidyverse_1.2.1  rmdformats_0.3.5 knitr_1.22      
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.0       lubridate_1.7.4  assertthat_0.2.0 digest_0.6.18   
##  [5] mime_0.6         R6_2.4.0         cellranger_1.1.0 plyr_1.8.4      
##  [9] backports_1.1.3  evaluate_0.13    ggstance_0.3.1   httr_1.4.0      
## [13] highr_0.7        pillar_1.3.1     rlang_0.3.1      lazyeval_0.2.1  
## [17] readxl_1.3.0     rstudioapi_0.9.0 miniUI_0.1.1.1   rmarkdown_1.12  
## [21] labeling_0.3     questionr_0.7.0  pander_0.6.3     munsell_0.5.0   
## [25] shiny_1.2.0      broom_0.5.1      compiler_3.5.2   httpuv_1.4.5.1  
## [29] modelr_0.1.4     xfun_0.5         pkgconfig_2.0.2  htmltools_0.3.6 
## [33] tidyselect_0.2.5 bookdown_0.9     crayon_1.3.4     withr_2.1.2     
## [37] later_0.8.0      grid_3.5.2       nlme_3.1-137     jsonlite_1.6    
## [41] xtable_1.8-3     gtable_0.2.0     magrittr_1.5     scales_1.0.0    
## [45] cli_1.0.1        stringi_1.3.1    promises_1.0.1   xml2_1.2.0      
## [49] generics_0.0.2   tools_3.5.2      glue_1.3.0       hms_0.4.2       
## [53] yaml_2.2.0       colorspace_1.4-0 rvest_0.3.2      haven_2.1.0