第11回講義資料

交互作用

スライド

セットアップ

　本日の実習用データ（LMSからダウンロード可能）と必要なパッケージ（{tidyverse}、{summarytools}、{marginaleffects}）を読み込む。ただし、{marginaleffects}がインストールされていない場合は、コンソール上にinstall.packages("marginaleffects")を入力し、インストールしておくこと。

Code 01

library(tidyverse)
library(summarytools)
library(marginaleffects)

jes_df <- read_csv("Data/JES6_W1_2.csv")

交互作用とは

　以下の図のように主な説明変数（\(X\)）と応答変数（\(Y\)）の関係において、\(X\)が\(Y\)に与える影響がその他の変数（\(Z\)）の影響を受ける場合を考えてみよう。

　これは\(X\)が\(Y\)に与える影響は一定ではないことを意味する（= \(X\)が\(Y\)に与える影響は\(Z\)に依存する）。ここで、\(Z\)は調整変数 (moderation variable; moderator)と呼ばれる。調整変数はダミー変数でも、連続変数でも可能だ。

　この交互作用を実際の回帰式として表すと以下のようになる。

\[ \hat{Y} = \alpha + \beta_1 X + \beta_2 Z + \beta_3 X \cdot Z \tag{1}\]

　ここで説明変数に調整変数をかけた変数（\(X \cdot Z\)）は交差項 (interaction term)と呼ばれる。これまでの（重）回帰分析では変数\(X\)が\(Y\)に与える効果は\(\beta_1\)であると解説したが、回帰式に交差項が含まれている場合は解釈に注意が必要だ。なぜなら回帰式において\(X\)は\(\beta_3\)にも登場するからである。これは式 1 を変形してみると分かりやすい。

\[ \hat{Y} = \alpha + (\beta_1 + \beta_3 Z) X + \beta_2 Z \tag{2}\]

　式 2 によると、変数\(X\)が\(Y\)に与える効果、つまり\(X\)の傾き係数は\(\beta_1\)でなく、\(\beta_1 + \beta_3 Z\)だということが分かる。交互作用を仮定したモデルに解釈はやや面倒ではないが、難しい数式が登場するわけではない。ここでは調整変数がダミー変数の場合と連続変数の場合の例を紹介する。

　まず、調整変数\(Z\)が0、または1の値のみをとるダミー変数の場合（\(Z \in \{0, 1\}\)）だ。回帰分析の結果、以下のような1次関数が得られたとしよう。

\[ \begin{align} \hat{Y} & = 3 + 2 X + 1 Z + 3 X \cdot Z \\ & = 3 + (2 + 3Z) X + 1 Z \end{align} \tag{3}\]

　解釈する場合は\(Z = 0\)の場合と\(Z = 1\)の場合に分けて解釈する。まず、\(Z = 0\)の場合、\(\hat{Y} = 3 + (2 + 3\cdot0) X + 1 Z = 3 + 2X + 1Z\)となり、\(X\)が\(Y\)に与える影響は2となる。一方、\(Z = 1\)の場合は\(\hat{Y} = 3 + (2 + 3\cdot1) X + 1 Z = 3 + 5X + 1Z\)となり、\(X\)が\(Y\)に与える影響は5となる。これを可視化すると以下のようになる。

　調整変数\(Z\)が0になると傾きが2の回帰直線（赤）が、1になると傾きが5の回帰直線（青）が得られる。

　続いて、調整変数\(Z\)が無数の値をとる連続変数の場合を考えてみよう。回帰分析の結果、以下のような1次関数が得られたとしよう。

\[ \begin{align} \hat{Y} & = 2 + 3 X + 2 Z - 1 X \cdot Z \\ & = 2 + (3 - 1Z) X + 2 Z \end{align} \tag{4}\]

　ここで\(Z = -1\)の場合、\(\hat{Y} = 2 + (3 - 1\cdot(-1)) X + 2 Z = 3 + 4X + 1Z\)となり、\(X\)が\(Y\)に与える影響は4となる。また、\(Z = 2\)の場合は\(\hat{Y} = 2 + (3 - 1\cdot2) X + 2 Z = 3 + 1X + 1Z\)となり、\(X\)が\(Y\)に与える影響は1となる。\(Z\)が3.5の場合は\(\hat{Y} = 2 + (3 - 1\cdot3.5) X + 2 Z = 3 - 0.5X + 1Z\)となり、\(X\)が\(Y\)に与える影響は-0.5となる。これを可視化すると以下のようになる。

　ただし、調整変数が連続変数の場合は、-1、2、3.5以外の値を取ることもできる。\(Z\)が取り得るすべての値に対して回帰直線を計算することは出来ないため、いくつかの代表的な値に絞って計算する必要があろう。

推定

　それでは実際に推定してみよう。まず、実習データの中身を確認してみる。

Code 02

jes_df

# A tibble: 3,000 × 6
   TempKyosan Female   Age Satisfaction Interest Ideology
        <dbl>  <dbl> <dbl>        <dbl>    <dbl>    <dbl>
 1         20      1    69            4        4        9
 2         20      1    47            1        1        7
 3          0      1    37            3        3       11
 4          0      0    51            4        3       11
 5         20      0    38            2        3        7
 6          0      0    71            5        4       11
 7         10      0    47            3        3        9
 8          0      0    71            4        4       11
 9         25      0    75            3        4        9
10         40      1    66            2        3        6
# … with 2,990 more rows

　各変数の説明は以下の通りだ。

変数	説明	備考
`TempKyosan`	日本共産党に対する感情温度	高いほど好感
`Female`	女性ダミー	0: 男性 / 1: 女性
`Age`	回答者の年齢
`Satisfaction`	政治満足度	高いほど満足
`Interest`	回答者の政治関心	高いほど関心あり
`Ideology`	回答者のイデオロギー	高いほど保守的

　データ分析の前にjes_dfの記述統計量を確認する。今回のデータはすべて連続変数扱いとなるため、前処理は不要だ。性別は名目変数であるが、既にダミー変数に変換済みである。ダミー変数の記述統計量は連続変数と同じ扱いで問題ないため、データをそのまま{summarytools}パッケージのdescr()関数に渡せば良い。

Code 03

jes_df %>%
  descr(stats = c("mean", "sd", "min", "max", "n.valid"),
        transpose = TRUE, order = "p")

Descriptive Statistics  
jes_df  
N: 3000  

                      Mean   Std.Dev     Min      Max   N.Valid
------------------ ------- --------- ------- -------- ---------
        TempKyosan   26.88     24.95    0.00   100.00   3000.00
            Female    0.50      0.50    0.00     1.00   3000.00
               Age   47.34     15.63   18.00    75.00   3000.00
      Satisfaction    2.45      1.08    1.00     5.00   3000.00
          Interest    2.74      0.83    1.00     4.00   3000.00
          Ideology    6.34      2.10    1.00    11.00   3000.00

　このデータを用い、以下の問いに答えるとする。

　政治満足度が共産党に対する感情温度に与える影響を調べたい。ただし、この影響は一定ではなく、性別や年齢によって異なるかも知れない。政治満足度が共産党に対する感情温度に与える影響の不均一性を調べるためにはどうすれば良いだろうか。

　この問いにおける応答変数は「共産党に対する感情温度（TempKyosan）」、主な説明変数は「政治満足度（Satisfaction）」だ。ただし、政治満足度が共産感情温度に与える影響は性別（Female）や年齢（Age）に依存する可能性がある。調整変数が2つであるため、ここでは2つのモデルを作成する。

モデル1

種類	変数	変数名
主な説明変数	政治満足度	`Satisfaction`
応答変数	共産党に対する感情温度	`TempKoysan`
調整変数	女性ダミー	`Female`
統制変数	政治関心、イデオロギー、年齢	`Interest`、`Ideology`、`Age`

モデル2

種類	変数	変数名
主な説明変数	政治満足度	`Satisfaction`
応答変数	共産党に対する感情温度	`TempKoysan`
調整変数	年齢	`Age`
統制変数	政治関心、イデオロギー、女性ダミー	`Interest`、`Ideology`、`Female`

　モデル1は政治満足度（= 説明変数）が共産党に対する感情温度（= 応答変数）に与える影響は性別（= 調整変数）に依存することを意味し、以下のような図として表現できる。

　これを回帰式にする場合、式 5 のような1次関数になる（以下の図参照）。

\[ \widehat{\mbox{TempKyosan}} = \alpha + \beta_1 \mbox{Satisfaction} + \beta_2 \mbox{Female} + \beta_3 \mbox{Interest} + \beta_4 \mbox{Ideology} + \beta_5 \mbox{Age} + \beta_6 (\mbox{Satisfaction} \cdot \mbox{Female}) \tag{5}\]

　つづいて、モデル2は政治満足度（= 説明変数）が共産党に対する感情温度（= 応答変数）に与える影響は年齢（= 調整変数）に依存することを意味し以下のような図として表現できる。

　これを回帰式にする場合、式 6 のような1次関数になる（以下の図参照）。

\[ \widehat{\mbox{TempKyosan}} = \alpha + \beta_1 \mbox{Satisfaction} + \beta_2 \mbox{Age} + \beta_3 \mbox{Interest} + \beta_4 \mbox{Ideology} + \beta_5 \mbox{Female} + \beta_6 (\mbox{Satisfaction} \cdot \mbox{Age}) \tag{6}\]

　推定にはこれまでと同様、lm()関数を使用する。注意するところは第1引数である回帰式（formula）であり、交互作用が存在すると考えられる2つの変数を+でなく、*でつなぐだけだ。ここでA * Bの意味はA\(\times\)bだけでなく、A、B、A\(\times\)Bが同時に投入することを意味する。

　モデル1と2をそれぞれfit1、fit2に格納する。

Code 04

fit1 <- lm(TempKyosan ~ Satisfaction * Female + Interest + Ideology + Age, 
           data = jes_df)
fit2 <- lm(TempKyosan ~ Satisfaction * Age + Interest + Ideology + Female + Age, 
           data = jes_df)

　2つのモデルの推定結果を横に並べるために{modelsummary}とmodelsummary()関数を使用する。2つのモデルオブジェクトはlist()関数でまとめることを忘れずに。

Code 05

modelsummary(list(fit1, fit2))

	Model 1	Model 2
(Intercept)	49.275	28.250
	(2.587)	(4.164)
Satisfaction	−4.732	3.308
	(0.561)	(1.354)
Female	2.752	5.061
	(2.186)	(0.888)
Interest	0.217	−0.118
	(0.571)	(0.570)
Ideology	−1.887	−1.690
	(0.214)	(0.215)
Age	−0.040	0.357
	(0.030)	(0.073)
Satisfaction × Female	0.902
	(0.817)
Satisfaction × Age		−0.157
		(0.026)
Num.Obs.	3000	3000
R2	0.086	0.096
R2 Adj.	0.084	0.094
AIC	27560.1	27526.4
BIC	27608.1	27574.4
Log.Lik.	−13772.027	−13755.176
F	46.763	52.927
RMSE	23.85	23.72

モデル1の解釈

　モデル1（fit1）から得られた回帰式は以下の通りである。

共産に対する感情温度の予測値 = 49.28 - 4.73 \(\times\) 政治満足度 + 2.75 \(\times\) 女性ダミー + 0.22 \(\times\) 政治関心 - 1.89 \(\times\) イデオロギー - 0.04 \(\times\) 年齢 + 0.90 \(\times\) 政治満足度 \(\times\) 女性ダミー

　これを政治満足度でまとめると、

共産に対する感情温度の予測値 = 49.28 - (4.73 - 0.90 \(\times\) 女性ダミー) \(\times\) 政治満足度 + 2.75 \(\times\) 女性ダミー + 0.22 \(\times\) 政治関心 - 1.89 \(\times\) イデオロギー - 0.04 \(\times\) 年齢

になる。つまり、政治満足度が共産感情温度に与える影響は「-(4.73 - 0.90 \(\times\) 女性ダミーの値)」だ。女性ダミーが取り得る値は0（男性）か1（女性）しかないので、ここでは2つのケースにおける政治満足度が共産感情温度に与える影響を計算してみよう。

　まず、男性の場合（Femaleの値 = 0）、政治満足度が共産感情温度に与える影響は約-4.73（= -(4.73 - 0.90 \(\times\) 0) ）である。これは男性の場合、政治満足度が上がると共産感情温度が下がることを意味し、具体的には「男性の場合、政治満足度が1単位上がると、共産に対する感情温度は約4.73度下がる」と解釈できる。また、女性の場合（Femaleの値 = 1）、政治満足度が共産感情温度に与える影響は約-3.83（= -(4.73 - 0.90 \(\times\) 1) ）である。これもまた、女性の場合、政治満足度が上がると共産感情温度が下がることを意味し、具体的には「女性の場合、政治満足度が1単位上がると、共産に対する感情温度は約3.83度下がる」と解釈できる。

　以上の結果を回帰直線で示したものが以下の図である。

モデル2の解釈

　モデル2（fit2）に対しても同じことが言える。まずは回帰式をまとめてみよう。

共産に対する感情温度の予測値 = 28.25 + 3.31 \(\times\) 政治満足度 + 0.36 \(\times\) 年齢 - 0.12 \(\times\) 政治関心 - 1.69 \(\times\) イデオロギー + 5.06 \(\times\) 女性ダミー - 0.16 \(\times\) 政治満足度 \(\times\) 年齢

　これを政治満足度でまとめると、

共産に対する感情温度の予測値 = 28.25 + (3.31 - 0.16 \(\times\) 年齢) \(\times\) 政治満足度 + 0.36 \(\times\) 年齢 - 0.12 \(\times\) 政治関心 - 1.69 \(\times\) イデオロギー + 5.06 \(\times\) 女性ダミー

になる。つまり、政治満足度が共産感情温度に与える影響は「(3.31 - 0.16 \(\times\) 年齢の値)」だ。年齢は18、19、20、…など様々な値を取り得る。これら全てに対して計算することは効率的ではないため、ここでは年齢が20歳、40歳、60歳の場合に絞って政治満足度が共産感情温度に与える影響を計算してみよう。

　まず、20歳の場合（Ageの値 = 20）、政治満足度が共産感情温度に与える影響は約0.11（= 3.31 - 0.16 \(\times\) 20）である（小数点7桁まで計算すると約0.1724813）。これは「20歳の場合、政治満足度が1単位上がると、共産に対する感情温度は約0.11度上がる」ことを意味する。つづいて、40歳の場合（Ageの値 = 40）、政治満足度が共産感情温度に与える影響は約-3.09（= 3.31 - 0.16 \(\times\) 40）である（小数点7桁まで計算すると約-2.962879）。これは「40歳の場合、政治満足度が1単位上がると、共産に対する感情温度は約3.09度下がる」ことを意味する。最後に60歳の場合（Ageの値 = 60）、政治満足度が共産感情温度に与える影響は約-6.29（= 3.31 - 0.16 \(\times\) 60）だ（小数点7桁まで計算すると約-6.098239）。これは「60歳の場合、政治満足度が1単位上がると、共産に対する感情温度は約6.29度下がる」ことを意味する。

　以上の結果を回帰直線で示したものが以下の図である。

　今回は計算を簡単にするために小数点2桁目まで丸めた係数を用いたが、より性格に計算するためにはcoef()関数から抽出した係数を用いた方が良い。fit2の係数を抽出してみよう。

Code 06

coef(fit2)

     (Intercept)     Satisfaction              Age         Interest 
      28.2497444        3.3078413        0.3572177       -0.1180184 
        Ideology           Female Satisfaction:Age 
      -1.6902682        5.0610601       -0.1567680

　ここで政治満足度（Satisfication）の係数は2番目であるため、coef(fit2)の後ろに[2]を付けると2番目の要素のみ抽出できる。

Code 07

coef(fit2)[2]

Satisfaction 
    3.307841

　同じく交差項の係数は7番目であるため、coef(fit2)[7]で抽出できる。

Code 08

coef(fit2)[7]

Satisfaction:Age 
       -0.156768

　この数値を使えば、政治満足度が共産感情温度に与える影響をより正確な数値で計算できる。たとえば、60歳における政治満足度が共産感情温度に与える影響は

Code 09

coef(fit2)[2] + coef(fit2)[7] * 60

Satisfaction 
   -6.098239

のように計算できる。

限界効果の話

　交互作用とは説明変数\(X\)が応答変数\(Y\)に与える影響が調整変数\(Z\)の値に依存することを意味する。この場合、「説明変数\(X\)が応答変数\(Y\)に与える影響」の統計的有意性はどうだろうか。\(X\)が\(Y\)に与える影響の統計的有意性を検定する際に用いられる検定統計量は「\(X\)の係数 / \(X\)の標準誤差」である。しかし、ここでの\(X\)の係数（と標準誤差）は\(Z\)の値によって変わる。これは\(Z\)の値によって\(X \rightarrow Y\)の統計的有意性は変わることを意味する。

　したがって、\(Z\)の値ごとに、\(X\)が\(Y\)に与える影響（= 限界効果; marginal effects）を計算するだけでなく、95%信頼区間 or \(p\)値も示す必要がある。95%信頼区間内に0が含まれないことは、\(p < 0.05\)を意味するため、統計的に有意な関連があると判断する。

　限界効果の計算方法については解説済みである。しかし、統計的有意性に関しては説明がやや難しくなるため、これはRパッケージにお任せするとしよう。もし、数学的な背景に関心のある履修者はBrambor (2006)を参照されたい。

　今回使用するパッケージ&関数は、冒頭で読み込んだ{marginaleffect}パッケージのmarginaleffects()関数だ。これは回帰モデルの限界効果と統計的有意性検定までやってくれる大変便利なパッケージである。使い方は以下の通り。

marginaleffects(回帰オブジェクト名, variables = "説明変数名", 
                newdata = datagrid(調整変数名 = 調整変数の値))

　fit1の場合、Femaleの値が0と1の場合のSatisfactonの限界効果を求めることになる。したがって、回帰オブジェクト名はfit1、variablesの実引数は"Satisfaction"（"で囲む）、datagrid()の中身はFemale = 0:1またはFemale = c(0, 1)だ。計算結果をfit1_ameに格納し、中身を確認する。

Code 10

fit1_ame <- marginaleffects(fit1, variables = "Satisfaction", 
                            newdata = datagrid(Female = c(0, 1)))

fit1_ame

  rowid     type         term      dydx std.error statistic      p.value
1     1 response Satisfaction -4.731682 0.5613401 -8.429261 3.478550e-17
2     2 response Satisfaction -3.830102 0.6110043 -6.268535 3.644599e-10
   conf.low conf.high predicted predicted_hi predicted_lo TempKyosan
1 -5.831888 -3.631476  24.41332     24.41143     24.41332     26.882
2 -5.027648 -2.632555  29.37297     29.37144     29.37297     26.882
  Satisfaction Interest Ideology   Age Female   eps
1        2.449 2.736333 6.340667 47.34      0 4e-04
2        2.449 2.736333 6.340667 47.34      1 4e-04

　非常に多くの情報が出力されるが、ここで興味のある列はFemaleの値（Female列）と、Satisfactionの限界効果（dydx列）、\(p\)値（p.value列）、95%信頼区間（conf.low列とconf.high列）である。それでは、select()関数でfit1_ameからこれらの列のみ出力してみよう。

　一つ注意してほしいのはselect()関数の呼び出し方である。Rでは同じ名前を持つ関数が複数のパッケージで提供される場合がある。これまで使ってきたselect()関数は{dplyr}パッケージのselect()関数だ。「え、私は{dplyr}パッケージなんか読み込んだ覚えないぞ」という人もいるかと思うが、{tidyverse}を読み込むと{dplyr}も読み込まれる（この場合、{tidyverse}は{dplyr}に依存すると言う）。しかし、select()関数は{MASS}というパッケージでも提供されており、これは{dplyr}のselect()関数とは別物である。実は、{marginaleffects}パッケージは{MASS}パッケージに依存している。したがって、「どのパッケージのselect()か」を具体的に指定しないと、Rは混乱してしまう。特定のパッケージの関数を指定するためにはパッケージ名::関数()と表記する。{MASS}も幅広くに使われ、多くのパッケージが{MASS}に依存する。したがって、select()関数は衝突しやすい関数であるため、普段からdplyr::select()という書き方をした方が良いかも知れない。

Code 11

fit1_ame %>%
  # p.value, conf.low, conf.high は隣接した列なので p.value:conf.high でOK
  dplyr::select(Female, dydx, p.value, conf.low, conf.high)

  Female      dydx      p.value  conf.low conf.high
1      0 -4.731682 3.478550e-17 -5.831888 -3.631476
2      1 -3.830102 3.644599e-10 -5.027648 -2.632555

　Female = 0の場合のSatisfactionの限界効果は約-4.732、Female = 1の場合のそれは約-3.830である。また、いずれも\(p\)値が0.05を下回ることから統計的に有意な関係であることがわかる。つまり、性別と関係なく、政治満足度は共産感情温度に負の影響を与えることが分かる。

　これらの結果を表でまとめると以下のようになる。表で\(p\)値を報告する場合、\(p\)値が非常に小さいケースがある（たとえば、\(p\) = 0.0000000235）。この場合、「\(p\) = 0.000」と表記せず、「\(p\) < 0.001」と表記すること。\(p\)値がぴったり0になることは実質的にあり得ない。

性別	平均限界効果	p値	95%信頼区間
性別	平均限界効果	p値	下限	上限
男性	−4.732	< 0.001	−5.832	−3.631
女性	−3.830	< 0.001	−5.028	−2.633

　続いて、fit2の場合のAgeの値が18、19、20、…、75の場合のSatisfactionの限界効果を求めるてみよう。回帰オブジェクト名はfit2、variablesの実引数は"Satisfaction"（"で囲む）、datagrid()の中身はAge = 18:75だ。結果を出力する際、今回はround()関数を使って\(p\)値を小数点3桁までにまとめてみよう。round()の第1引数は丸めたい変数名、第2引数は小数点の桁数である。

Code 12

fit2_ame <- marginaleffects(fit2, variables = "Satisfaction", 
                            newdata = datagrid(Age = 18:75))

fit2_ame %>%
  dplyr::select(Age, dydx, p.value, conf.low, conf.high) %>%
  mutate(p.value = round(p.value, 3)) # p値を小数点3桁で丸める

   Age        dydx p.value   conf.low   conf.high
1   18  0.48601719   0.594  -1.303101  2.27513577
2   19  0.32924919   0.711  -1.413868  2.07236610
3   20  0.17248118   0.842  -1.524974  1.86993684
4   21  0.01571317   0.985  -1.636448  1.66787385
5   22 -0.14105484   0.863  -1.748323  1.46621293
6   23 -0.29782284   0.709  -1.860632  1.26498639
7   24 -0.45459085   0.557  -1.973416  1.06423447
8   25 -0.61135886   0.417  -2.086713  0.86399499
9   26 -0.76812687   0.293  -2.200571  0.66431691
10  27 -0.92489487   0.192  -2.315040  0.46525006
11  28 -1.08166288   0.116  -2.430183  0.26685679
12  29 -1.23843089   0.063  -2.546063  0.06920099
13  30 -1.39519890   0.031  -2.662747 -0.12765044
14  31 -1.55196690   0.013  -2.780317 -0.32361652
15  32 -1.70873491   0.005  -2.898858 -0.51861172
16  33 -1.86550292   0.002  -3.018471 -0.71253512
17  34 -2.02227093   0.000  -3.139260 -0.90528203
18  35 -2.17903894   0.000  -3.261345 -1.09673307
19  36 -2.33580694   0.000  -3.384850 -1.28676404
20  37 -2.49257495   0.000  -3.509917 -1.47523333
21  38 -2.64934296   0.000  -3.636693 -1.66199246
22  39 -2.80611097   0.000  -3.765345 -1.84687694
23  40 -2.96287897   0.000  -3.896040 -2.02971777
24  41 -3.11964698   0.000  -4.028951 -2.21034272
25  42 -3.27641499   0.000  -4.164259 -2.38857123
26  43 -3.43318300   0.000  -4.302139 -2.56422714
27  44 -3.58995100   0.000  -4.442765 -2.73713655
28  45 -3.74671901   0.000  -4.586295 -2.90714257
29  46 -3.90348702   0.000  -4.732869 -3.07410516
30  47 -4.06025503   0.000  -4.882597 -3.23791326
31  48 -4.21702303   0.000  -5.035562 -3.39848441
32  49 -4.37379104   0.000  -5.191808 -3.55577396
33  50 -4.53055905   0.000  -5.351344 -3.70977442
34  51 -4.68732706   0.000  -5.514135 -3.86051920
35  52 -4.84409507   0.000  -5.680111 -4.00807938
36  53 -5.00086307   0.000  -5.849168 -4.15255820
37  54 -5.15763108   0.000  -6.021175 -4.29408711
38  55 -5.31439909   0.000  -6.195979 -4.43281922
39  56 -5.47116710   0.000  -6.373412 -4.56892224
40  57 -5.62793510   0.000  -6.553297 -4.70257280
41  58 -5.78470311   0.000  -6.735457 -4.83394883
42  59 -5.94147112   0.000  -6.919714 -4.96322790
43  60 -6.09823913   0.000  -7.105898 -5.09058066
44  61 -6.25500713   0.000  -7.293842 -5.21617238
45  62 -6.41177514   0.000  -7.483393 -5.34015699
46  63 -6.56854315   0.000  -7.674411 -5.46267557
47  64 -6.72531116   0.000  -7.866761 -5.58386096
48  65 -6.88207916   0.000  -8.060326 -5.70383261
49  66 -7.03884717   0.000  -8.254991 -5.82270311
50  67 -7.19561518   0.000  -8.450659 -5.94057103
51  68 -7.35238319   0.000  -8.647238 -6.05752826
52  69 -7.50915119   0.000  -8.844649 -6.17365348
53  70 -7.66591920   0.000  -9.042817 -6.28902160
54  71 -7.82268721   0.000  -9.241677 -6.40369745
55  72 -7.97945522   0.000  -9.441167 -6.51774362
56  73 -8.13622323   0.000  -9.641234 -6.63121250
57  74 -8.29299123   0.000  -9.841828 -6.74415400
58  75 -8.44975924   0.000 -10.042910 -6.85660863

　中級者以上の場合、mutate()内にacross()とwhere()を組み合わせて「数値型の列全てに対し、round(変数, 3)を適用する」といった処理もできる。覚える必要はないが、「Rを勉強していくと、コードがより効率的に書けるようになる」ということだけ覚えておこう。

Code 13

fit2_ame %>%
  dplyr::select(Age, dydx:conf.high) %>%
  mutate(across(where(is.numeric), round, 3))

   Age   dydx std.error statistic p.value conf.low conf.high
1   18  0.486     0.913     0.532   0.594   -1.303     2.275
2   19  0.329     0.889     0.370   0.711   -1.414     2.072
3   20  0.172     0.866     0.199   0.842   -1.525     1.870
4   21  0.016     0.843     0.019   0.985   -1.636     1.668
5   22 -0.141     0.820    -0.172   0.863   -1.748     1.466
6   23 -0.298     0.797    -0.374   0.709   -1.861     1.265
7   24 -0.455     0.775    -0.587   0.557   -1.973     1.064
8   25 -0.611     0.753    -0.812   0.417   -2.087     0.864
9   26 -0.768     0.731    -1.051   0.293   -2.201     0.664
10  27 -0.925     0.709    -1.304   0.192   -2.315     0.465
11  28 -1.082     0.688    -1.572   0.116   -2.430     0.267
12  29 -1.238     0.667    -1.856   0.063   -2.546     0.069
13  30 -1.395     0.647    -2.157   0.031   -2.663    -0.128
14  31 -1.552     0.627    -2.476   0.013   -2.780    -0.324
15  32 -1.709     0.607    -2.814   0.005   -2.899    -0.519
16  33 -1.866     0.588    -3.171   0.002   -3.018    -0.713
17  34 -2.022     0.570    -3.548   0.000   -3.139    -0.905
18  35 -2.179     0.552    -3.946   0.000   -3.261    -1.097
19  36 -2.336     0.535    -4.364   0.000   -3.385    -1.287
20  37 -2.493     0.519    -4.802   0.000   -3.510    -1.475
21  38 -2.649     0.504    -5.259   0.000   -3.637    -1.662
22  39 -2.806     0.489    -5.734   0.000   -3.765    -1.847
23  40 -2.963     0.476    -6.223   0.000   -3.896    -2.030
24  41 -3.120     0.464    -6.724   0.000   -4.029    -2.210
25  42 -3.276     0.453    -7.233   0.000   -4.164    -2.389
26  43 -3.433     0.443    -7.744   0.000   -4.302    -2.564
27  44 -3.590     0.435    -8.251   0.000   -4.443    -2.737
28  45 -3.747     0.428    -8.747   0.000   -4.586    -2.907
29  46 -3.903     0.423    -9.225   0.000   -4.733    -3.074
30  47 -4.060     0.420    -9.677   0.000   -4.883    -3.238
31  48 -4.217     0.418   -10.098   0.000   -5.036    -3.398
32  49 -4.374     0.417   -10.480   0.000   -5.192    -3.556
33  50 -4.531     0.419   -10.819   0.000   -5.351    -3.710
34  51 -4.687     0.422   -11.111   0.000   -5.514    -3.861
35  52 -4.844     0.427   -11.357   0.000   -5.680    -4.008
36  53 -5.001     0.433   -11.554   0.000   -5.849    -4.153
37  54 -5.158     0.441   -11.706   0.000   -6.021    -4.294
38  55 -5.314     0.450   -11.815   0.000   -6.196    -4.433
39  56 -5.471     0.460   -11.885   0.000   -6.373    -4.569
40  57 -5.628     0.472   -11.920   0.000   -6.553    -4.703
41  58 -5.785     0.485   -11.925   0.000   -6.735    -4.834
42  59 -5.941     0.499   -11.904   0.000   -6.920    -4.963
43  60 -6.098     0.514   -11.861   0.000   -7.106    -5.091
44  61 -6.255     0.530   -11.801   0.000   -7.294    -5.216
45  62 -6.412     0.547   -11.727   0.000   -7.483    -5.340
46  63 -6.569     0.564   -11.642   0.000   -7.674    -5.463
47  64 -6.725     0.582   -11.548   0.000   -7.867    -5.584
48  65 -6.882     0.601   -11.448   0.000   -8.060    -5.704
49  66 -7.039     0.620   -11.344   0.000   -8.255    -5.823
50  67 -7.196     0.640   -11.237   0.000   -8.451    -5.941
51  68 -7.352     0.661   -11.129   0.000   -8.647    -6.058
52  69 -7.509     0.681   -11.020   0.000   -8.845    -6.174
53  70 -7.666     0.703   -10.912   0.000   -9.043    -6.289
54  71 -7.823     0.724   -10.805   0.000   -9.242    -6.404
55  72 -7.979     0.746   -10.699   0.000   -9.441    -6.518
56  73 -8.136     0.768   -10.596   0.000   -9.641    -6.631
57  74 -8.293     0.790   -10.494   0.000   -9.842    -6.744
58  75 -8.450     0.813   -10.395   0.000  -10.043    -6.857

　18〜29歳の場合、政治満足度が共産感情温度に与える影響は確認できない（つまり、\(p \geq 0.05\)）。一方、回答者の年齢が30歳以上の場合、政治満足度は共産感情温度に負の影響を与えるが確認できる（つまり、\(p < 0.05\)）。

　交互作用を含む回帰モデルの場合、調整変数ごとの予測値の図と限界効果のグラフを示すのが一般的である。回帰分析の結果の可視化は第14回で解説するが、ここでは限界効果のグラフの読み方を紹介しよう。

　以下の図はfit1における政治満足度の（平均）限界効果である。男性でも女性でも95%信頼区間内に0が含まれておらず、どちらも負の限界効果である。これは「性別と関係なく政治満足度が高まると共産党に感情温度は下がる」ことを意味する。

コード

fit1_ame %>%
  mutate(Female = factor(Female, levels = 0:1, 
                         labels = c("男性", "女性"))) %>%
  ggplot() +
  geom_hline(yintercept = 0) +
  geom_pointrange(aes(x = Female, y = dydx, 
                      ymin = conf.low, ymax = conf.high)) +
  labs(x = "性別", y = "政治満足度の平均限界効果と95%信頼区間") +
  theme_bw()

　続いて、fit2における政治満足度の（平均）限界効果である。ここでは95%信頼区間が垂直線でなく、領域（ribbon）で示される。29歳までは95%信頼区間に0が含まれている。30歳からは限界効果が負となり、統計的に有意な関係を示している。。これは「30歳以上の有権者のみにおいて、政治満足度が高まると共産党に感情温度は下がる」ことを意味する。

コード

fit2_ame %>%
  ggplot(aes(x = Age)) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill = "gray80") +
  geom_hline(yintercept = 0) +
  geom_line(aes(y = dydx)) +
  labs(x = "年齢", y = "政治満足度の平均限界効果と95%信頼区間") +
  theme_bw()