回帰分析の解釈と限界効果

特殊講義「インターネットを活用した経済データの分析」講義資料

2/16/2017

1 回帰分析,回帰係数,限界効果

講義スライドで説明した通り,回帰分析の結果を解釈する上では単に回帰係数やその統計的有意性をみるだけではなく,実質的有意性や効果の大きさを計算する必要がある.この際,独立変数1 (単位) あたりの従属変数の変化量を示す限界効果 (marginal effect)を計算することが多い.また,単に限界効果の値を表や文章で提示するだけでは分かりにくいので,グラフの積極的な使用を推奨する.

Rを用いてある回帰モデルについて限界効果を計算しその結果を図示するには,R言語や回帰分析,行列計算に習熟している必要がある.とはいえ,限界効果の計算や図示を簡単に行なえる優れたパッケージが利用可能なので,講義時間の都合もあり,ここではパッケージを用いた最低限のコードを提示する.しかし,計量分析の手法やR言語に習熟するためには,なるべくこうしたパッケージに頼らず (計算過程やコードをブラックボックス化せず),自分で1から計算した方が理解が深まり,望ましい.

講義で説明した通り,線形回帰においても交互作用項を含む回帰式について限界効果を計算する際には,特に注意が必要になる. 以下では,(1) 交互作用項を含まない回帰式と (2) 交互作用項を含む回帰式のそれぞれについて,限界効果の計算・図示する.

実際の作業に入る前に,使用するパッケージとサンプルデータを読み込んでおく.

library(interplot)
library(sjPlot)
library(sjmisc)
library(stargazer)
library(stringr)
library(tidyverse)
# load saple datasets
data(efc)
data(mtcars)
mtcars = tbl_df(mtcars) %>%
    mutate_at(vars(am, gear), funs(factor)) # convert into factor
# set sjPlot options
set_theme(
    "forest",
    axis.title.size = .85,
    axis.textsize = .85,
    geom.errorbar.size = .5,
    legend.size = .85,
    geom.label.size = 3.5
    )

1.1 交互作用を含まない回帰式

サンプルデータmtcarsを使って,交互作用を含まない回帰分析を回し,mod_0オブジェクトに保存する. このサンプルデータの詳細は,Rのコンドールに?mtcarsと入力すれば説明が表示される.

stargazer(as.data.frame(mtcars), type="text") # summary statistics
## 
## ============================================
## Statistic N   Mean   St. Dev.  Min     Max  
## --------------------------------------------
## mpg       32 20.091   6.027   10.400 33.900 
## cyl       32  6.188   1.786     4       8   
## disp      32 230.722 123.939  71.100 472.000
## hp        32 146.688  68.563    52     335  
## drat      32  3.597   0.535   2.760   4.930 
## wt        32  3.217   0.978   1.513   5.424 
## qsec      32 17.849   1.787   14.500 22.900 
## vs        32  0.438   0.504     0       1   
## carb      32  2.812   1.615     1       8   
## --------------------------------------------
fm_0 = as.formula(mpg ~ qsec + wt + am) # set regression model
mod_0 = lm(fm_0, data = mtcars) # run regression
stargazer(mod_0, type = "text", single.row=TRUE,
          ci = TRUE # print 95% CIs
          )
## 
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                                 mpg            
## -----------------------------------------------
## qsec                  1.226*** (0.660, 1.792)  
## wt                  -3.917*** (-5.310, -2.523) 
## am1                   2.936** (0.171, 5.701)   
## Constant              9.618 (-4.023, 23.258)   
## -----------------------------------------------
## Observations                    32             
## R2                             0.850           
## Adjusted R2                    0.834           
## Residual Std. Error       2.459 (df = 28)      
## F Statistic           52.750*** (df = 3; 28)   
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

stargazer()関数はci=TRUEという引数を与えると回帰係数の95%信頼区間を出力してくれる.下の図では95%信頼区間をプロットするので,ここでは比べやすいように表でも95%信頼区間を出力した.

1.1.1 統計的有意性と回帰係数のプロット

さて,モデルに意味があるかはともかくとして,ひとまずサンプルデータを用いた回帰分析の結果が保存できた.また,mod_0オブジェクトの出力結果から,いくつか統計的に有意な係数があることもわかる.

しかし,表だけではどうしても分かりにくい.まず,sjPlot::sjp.lm()関数を用いて,回帰係数を図示しよう.

sjp.lm(mod_0)

これだけで,上の図が出力される (フォレスト・プロットと呼ぶ).点は回帰係数の推定値を,実線は回帰係数の95%信頼区間を示す. この図では,上のstargazer(mod_0, type = "text", single.row=TRUE)等で出力される統計量のいくつかは省略されてしまうが,回帰係数の大小を視覚的に提示するために表とフォレスト・プロットの両方を提示することが望ましい (あるいは,図を工夫して,余白に省略された統計量を含めてしまえばよい).

また,ここでは問題になっていないが,単位が異なる独立変数を投入した回帰モデルでは,フォレスト・プロットが見にくくなることがしばしばある (e.g., cm単位の独立変数とkm単位の独立変数).その場合は,標準化偏回帰係数をプロットすればよい.上のコードに1つだけオプションを足すと,標準化偏回帰係数をプロットできる.

sjp.lm(mod_0, type="std")

x軸のラベルが自動的に変更されていることに注意しよう. なお,標準化偏回帰係数は「ある独立変数が 1標準偏差分変動したとき,標準化された従属変数が何単位変動するか」を示す.標準化偏回帰係数を用いることで,異なる単位の変数の影響の大きさを,同じ尺度で示すことができる.

1.1.2 実質的有意性と限界効果のプロット

さて,講義スライドでも説明した通り,データ解析においては統計的有意性 (statistical significance) と同程度以上に,実質的有意性 (substantial significance) が重要になる.上のフォレストプロットは統計的有意性を見分ける上で便利だが,実質的有意性を検討するためには限界効果を計算して図示することが望ましい.

実は,限界効果の計算と図示も,上記のコードを一箇所変えてやるだけで出力できる.たとえば,上記の回帰式に含まれるwt変数の限界効果は,次のコードで計算・プロットできる.

sjp.lm(mod_0, type="eff", facet.grid = FALSE, vars = "wt")

この図は,横軸にwt変数の値を,縦軸にはwt変数の値に対応する従属変数mpgの予測値をプロットしている (他の変数は平均値に固定).また,グレーのシェードは95%信頼区間を示す. wt変数の最小値と最大値は1.513, 5.424なので,横軸の範囲もこれに対応している (回帰分析の前に確認している要約統計量を確認のこと).

上のコードのvars = "wt"という部分を他の独立変数に対応するように変えてやれば,異なる独立変数の限界効果を計算・プロットできる.

sjp.lm(mod_0, type="eff", facet.grid = FALSE, vars = "qsec")

こうした限界効果を示す図を用いることで,「ある独立変数を任意の区間 (ここでは最小値から最大値) で動かしたとき,従属変数がどの程度変化するか (影響を受けるか)」を視覚的に示すことができる.

1.2 交互作用を含む回帰式

1.2.1 交互作用がある場合の限界効果

さて,上の例では交互作用項を含まない回帰モデルを想定していた.では,交互作用項を含む回帰モデルでは,どのような作業が必要になるのだろうか. 講義スライドの通り,交互作用項を含まない線形回帰モデル,たとえば

\[\hat{y} = \beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2} +\epsilon\]

では,ある独立変数\(x_1\)の限界効果は,

\[\frac{\partial \hat{y}}{\partial x_{1}} = \beta_{1}\]

となる.他方,たとえば

\[\hat{y} = \beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2}+\beta_{3}x_{1}x_{2} +\epsilon\]

のような交互作用項を含む線形回帰モデルでは,ある独立変数\(x_1\)の限界効果は,

\[\frac{\partial \hat{y}}{\partial x_{1}} = \beta_{1} + \beta_{3}x_{2}\]

となる.つまり,交互作用を構成する独立変数\(x_1\)の限界効果は,交互作用項の回帰係数\(\beta_3\)ともう1つの独立変数\(x_2\)の値に依存する (交互作用がある).

交互作用を含む回帰式では,限界効果をどのように示せばよいのだろうか.直感的にいえば,\(x_2\)を任意の (いくつかの) 値に固定したときの\(x_1\)の限界効果を示すプロットを多数作成することが必要になる.

1.2.2 交互作用を踏まえたプロット

mtcarsデータセットを用いて,交互作用を含む回帰モデルを考えよう.

set_label(mtcars$am) <- "Automatic (0) or Manual (1)"
m_am <- lm(mpg ~ am + wt + am * wt, data = mtcars)
stargazer(m_am, type = "text", single.row=TRUE)
## 
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                                 mpg            
## -----------------------------------------------
## am1                      14.878*** (4.264)     
## wt                       -3.786*** (0.786)     
## am1:wt                   -5.298*** (1.445)     
## Constant                 31.416*** (3.020)     
## -----------------------------------------------
## Observations                    32             
## R2                             0.833           
## Adjusted R2                    0.815           
## Residual Std. Error       2.591 (df = 28)      
## F Statistic           46.567*** (df = 3; 28)   
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

Rでは,交互作用項はwt*amのように表現する.なお,ここでの書き方であればas.formula(mpg ~ geat + wt + gear:wt)と,:を用いてもよい.また,as.formula(mpg ~ gear*wt) としても同様の結果が得られる.しかし,as.formula(mpg ~ gear:wt)としてしまうとwt変数とam変数が回帰式から除かれてしまうため,注意しよう.

このモデルは交互作用項を含むので,上述の通りgear変数にせよwt変数にせよ,一方の変数の限界効果を計算・プロットするには他方の値を所与としなければならない. 手計算でこれを行なうとそれなりに手間がかかるが,sjPlot::sjp.int()関数を用いると簡単に行なえる.

sjp.int(m_am, type = "eff", facet.grid = FALSE, show.ci = TRUE) # two lines within single pane

sjp.int(m_am, type = "eff", facet.grid = TRUE, show.ci = TRUE) # two separated panes

sjPlot::sjp.int()関数は様々なオプションをもつが,type = "eff"という引数を与えると限界効果を計算してくれる.また,facet.gridFALSEにするかTRUEにするかで,上の図と下の図を描き分けることができる 上のプロットには2つのパネルがあるが,それぞれ次のようにam変数の値 (0, 1いずれかの値をとる) を所与とし,wt変数を最小値から最大値まで変化させたとき (横軸) の従属変数mpgの予測値 (縦軸) にプロットしている:

  • 左:am変数が0のとき
  • 右:am変数が1のとき

このプロットは少々見にくいので,0, 1という頭のラベルを変える.また,色も変えてみよう.なお,なぜ“Automatic”と“Manual”が0と1に対応するのかは,Rのコンソールに?mtcarsと入力し,このデータセットに含まれる変数の説明を読めば分かる.

sjp.int(m_am, type = "eff",
        facet.grid = TRUE, show.ci = TRUE,
        geom.colors = rep("black",2),
        legend.labels = c("Automatic", "Manual")
        )

この図から,am変数が“Manual”のとき,“Automatic”のときに比べて予測値を示す実線の傾きが大きくなっている.この図によって,独立変数wtの変化が従属変数mpgに与える影響はam変数が“Manual”のときの方が強くなることを,視覚的に示すことができる.

さて,上の分析結果は,am変数の効果もwt変数の値によって変化することを示唆する.さらに,回帰分析の出力を見るとwt変数と交互作用項の回帰係数は負だが,am変数の回帰係数は正になっている.ということは,wt変数の値によって,am変数が従属変数mpgに与える影響が正になる場合も負になる場合もあるということだろうか.

デフォルトでは,sjPlot::sjp.int()関数は「交互作用項を構成する変数 (ここではam*wtという交互作用項を投入したのでamwt) のうち,最初の変数 (ここではam変数)」 を,「プロットを分ける変数」として使用する (なので,上のプロットでは左右2枚がamによって分けられ,出力される).交互作用を含めたam変数の影響をプロットするには,swap.pred = TRUEという引数を足してやればよい.

sjp.int(m_am, type = "eff",
        facet.grid = TRUE, show.ci = TRUE,
        swap.pred = TRUE,
        legend.labels = c("wt=min(wt)", "wt=max(wt)")
        )

左のパネルはwt変数が最小値のときの,右のパネルはwt変数が最大値のときの結果をプロットしている. また,am変数は0 (“automatic”), 1 (“transmission”) という2つの値のいずれかをとる (ダミー変数という) ので,それぞれのパネルでは2つの予測値のみが描かれている. 上の予想通り,amが0から1に変化したとき,従属変数mpgに与える影響はwtの値によって正負の符号さえ異なる様子がみてとれる (この点を確認するには計算結果と図を追加した方がよいが,ここでは省略する).

なお,このプロットはwt変数が最小値のときと最大値のときを想定して描かれているが,mdrt.valuesという引数を変更すればこれを変更することもできる.Rのコンソールに?sjp.intと入力し,確認しておこう.

最後に,同じくmtcarsデータセットを用いるが別の交互作用を含む回帰分析の結果とプロットを,いくつか示しておく. 各自,出力結果が何を示すか解釈してみよう (二乗項もある種の交互作用項).

# An alternative interaction
m_gear <- lm(mpg ~ gear + wt + gear*wt, data = mtcars)
# stargazer(m_gear, type = "text", single.row=TRUE)
sjp.lm(m_gear) # forest plot

sjp.int(m_gear, type = "eff", # plot interaction
        facet.grid = TRUE, show.ci = TRUE,
        geom.colors = rep("black",3),
        legend.labels = str_c(rep("Gear ",3), seq(3,5))
        )

# Quadratic term
m_quad <- lm(mpg ~ gear + wt + I(wt^2), data = mtcars)
# sjp.lm(m_quad) ## forest plot
# plot interaction
interplot(m_quad, var1 = "wt", var2 = "wt") + 
    geom_hline(yintercept = 0) + 
    xlab("Automobile Weight (thousands lbs)") +
    ylab("Estimated Coefficient for\n Automobile Weight")