回帰分析とシミュレーション

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

2/16/2017

1 サンプルデータを用いた回帰分析

ここでは,carパッケージで利用できるサンプルデータPrestigeを用いて,R上で簡単な線形回帰モデルを実行する. Prestigeは1971年のカナダのセンサスデータに基づくデータセットである.次のコードを実行し,必要なパッケージとサンプルデータを呼び出す. これらのパッケージをインストールしていない場合は,各自インストールする.

library(arm)
library(car)
library(stargazer)
library(stringr)
library(tidyverse)
select = dplyr::select
data(Prestige)
(pres_tbl = tbl_df(Prestige))
## # A tibble: 102 × 6
##    education income women prestige census   type
## *      <dbl>  <int> <dbl>    <dbl>  <int> <fctr>
## 1      13.11  12351 11.16     68.8   1113   prof
## 2      12.26  25879  4.02     69.1   1130   prof
## 3      12.77   9271 15.70     63.4   1171   prof
## 4      11.42   8865  9.11     56.8   1175   prof
## 5      14.62   8403 11.68     73.5   2111   prof
## 6      15.64  11030  5.13     77.6   2113   prof
## 7      15.09   8258 25.65     72.6   2133   prof
## 8      15.44  14163  2.69     78.1   2141   prof
## 9      14.52  11377  1.03     73.1   2143   prof
## 10     14.64  11023  0.94     68.8   2153   prof
## # ... with 92 more rows

102の職業における教育水準や収入,女性の割合等が記録されたデータが読み込めた.なお,help(Prestige)とRに入力すると,データの詳細が表示される.

1.1 データの確認

このデータは,どのような姿をしているのだろうか.回帰分析に進む前に,まずはデータの構造をざっくりと把握するため,要約統計量をチェックし,を描こう. 要約統計量の計算と出力にはデフォルトのsummary()関数を用いてもよいが,回帰分析の出力にも使えるstargazer::stargazer()関数を使う.stargazer::stargazer()関数は,テキスト,html, LaTeXといった多様な出力フォーマットが指定できることに加え,さらにAER (American Economic Review) のような学術誌のフォーマットに沿った表を出力してくれる.この関数にPrestigeのようなデータを引数に与え,特にオプションを指定しなければ要約統計を出力する. また,図 (散布図行列) についてもデフォルトのpairs()関数があるが,car::scatterplotMatrix()関数を使うとより情報量の多い散布図行列 (近似線付き) を簡単に描ける (色やドット,文字が美しくないが,今は気にしない).

なお,census変数は職業の種類を示す変数 (解析には使わない) なので,ついでにデータから削除してしまおう.また,収入を示すincome変数は分布がやや右に歪んでいるので,自然対数をとっておく.

pres_tbl = pres_tbl %>%
    select(-census) %>%     # census変数を落とす
    mutate(income = log(income)) # 自然対数をとる
stargazer(as.data.frame(pres_tbl), type="text")   # 要約統計量
## 
## ===========================================
## Statistic  N   Mean  St. Dev.  Min    Max  
## -------------------------------------------
## education 102 10.738  2.728   6.380  15.970
## income    102 8.661   0.591   6.415  10.161
## women     102 28.979  31.725  0.000  97.510
## prestige  102 46.833  17.204  14.800 87.200
## -------------------------------------------
scatterplotMatrix(~ prestige + income + education + women + type, span = 0.5, data = pres_tbl, cex.labels = 1.1, cex.axis=.85)

この2つの関数の詳細については,それぞれ?stargazer?scatterplotMatrixとコンソールに入力して確認しておこう.特にscatterplotMatrix()関数については,なぜこれらの変数がプロットされているのか,また新たな引数を足すことができるのかといった点を確認するとよい. なお,実際のデータ解析でも回帰分析や高度な解析に進む前に,散布図や要約統計量,ヒストグラム等を用いてデータの姿を確認することがまず必要になる.

1.2 回帰分析

さて,このデータを用いて線形回帰モデルを考えよう.ここではまずprestigeを従属変数に,educationを独立変数を投入した (単) 回帰モデルを考える. Rで線形回帰 (OLS) を行なうには,lm()関数を用いる.

mod_0 <- lm(prestige ~ education, data= pres_tbl)

分析結果はmod_0に格納されている.結果を確認するには,デフォルトのsummary()関数かarm::display()関数を用いる.summary()関数を使うと有意確率を示す星にばかり目がいき,また出力もゴチャゴチャしていて見にくいので,ミニマルな出力のarm::display()関数の使用を勧める.

summary(mod_0)
## 
## Call:
## lm(formula = prestige ~ education, data = pres_tbl)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -26.0397  -6.5228   0.6611   6.7430  18.1636 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -10.732      3.677  -2.919  0.00434 ** 
## education      5.361      0.332  16.148  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.103 on 100 degrees of freedom
## Multiple R-squared:  0.7228, Adjusted R-squared:   0.72 
## F-statistic: 260.8 on 1 and 100 DF,  p-value: < 2.2e-16
display(mod_0)
## lm(formula = prestige ~ education, data = pres_tbl)
##             coef.est coef.se
## (Intercept) -10.73     3.68 
## education     5.36     0.33 
## ---
## n = 102, k = 2
## residual sd = 9.10, R-Squared = 0.72

上のモデルではeducationのみを独立変数とした単回帰モデルを推定したが,2以上の独立変数を含む重回帰モデルもlm()関数で実行できる. また,stargazer()関数を用いると複数の回帰モデルの出力を見やすく並べられるので,試しておこう.

mod_1 <- lm(prestige ~ education + income, data= pres_tbl)
mod_2 <- lm(prestige ~ education + income + women, data= pres_tbl)
stargazer(mod_0, mod_1, mod_2, type = "text", single.row = TRUE)
## 
## ============================================================================================
##                                               Dependent variable:                           
##                     ------------------------------------------------------------------------
##                                                     prestige                                
##                               (1)                      (2)                     (3)          
## --------------------------------------------------------------------------------------------
## education               5.361*** (0.332)        4.002*** (0.312)        3.731*** (0.354)    
## income                                          11.437*** (1.437)       13.438*** (1.914)   
## women                                                                     0.047 (0.030)     
## Constant               -10.732*** (3.677)      -95.194*** (10.998)    -110.966*** (14.843)  
## --------------------------------------------------------------------------------------------
## Observations                  102                      102                     102          
## R2                           0.723                    0.831                   0.835         
## Adjusted R2                  0.720                    0.828                   0.830         
## Residual Std. Error     9.103 (df = 100)         7.145 (df = 99)         7.093 (df = 98)    
## F Statistic         260.751*** (df = 1; 100) 243.323*** (df = 2; 99) 165.428*** (df = 3; 98)
## ============================================================================================
## Note:                                                            *p<0.1; **p<0.05; ***p<0.01

stargazer()関数を用いる利点の一つは,こうしたモデルの比較を簡単に出力することができることにある.

2 シミュレーションによる回帰分析の理解

回帰分析は,モデルが正しく特定されていれば (かつ種々の仮定が満たされていれば),各要因 (独立変数) の影響の大きさ (回帰係数) やその不確実性 (信頼区間) を正しく推定してくれるのだろうか.また,現実のデータの解析では,当然ながら「正しい」回帰式は事前に分からない (分かるなら解析しようとも思わない).「正しくない回帰式」を用いた推定は,どのような問題を引き起こすのだろうか.

ここでは,以下の手順で簡単なシミュレーションを行ない,回帰分析 (線形回帰) を直感的に理解し,上記の疑問に回答しよう.

  1. OLSが想定するデータ生成プロセス (data generation process) に従った仮想のデータを多数 (ここでは3,000) 生成する
  2. 作成した仮想データのそれぞれについて,「正しい」回帰式を用いた推定を行なう (拡張として,「正しくない」回帰式を用いた推定も行なう)
  3. 推定された回帰係数が,「真の」影響 (データを作成したときのパラメータの値) を正しく捉えているかをみる
  4. 3,000の仮想データを用いた解析のうち,何%の解析がデータ生成メカニズムを正しく捉えているかをまとめる

なお,以下のシミュレーションは,A. Gelman & J. Hill. 2007. Data analysis using regression and multilevel/hierarchical models. Cambridge: Cambridge Univ Press, Chap.8と国際大学の矢内勇生氏の講義資料を参照している.以下で行なう,仮想データを用いて「統計手法により推定したパラメータ」を「真のパラメータ」(母数) を比較することで,当該手法の妥当性や性質を検討・理解するシミュレーションは一般的に “fake-data simulation” と呼ばれる (Gelman & Hill 2007: 155).

2.1 仮想データの作成

上記のように,ここでのシミュレーションの目的は仮想データを用いて回帰分析によって推定したパラメータと真のパラメータを比較することにある. 講義スライドの通り,線形回帰モデルは次のように表現できる.

\[\boldsymbol{y} \sim N(\boldsymbol{X}\boldsymbol{\beta}, \sigma^{2} \boldsymbol{I})\]

ここでは,まずこうした回帰モデルが想定するデータを人工的に生成する.回帰モデルが真のパラメータを「正しく」推定するなら,回帰分析によって推定されたパラメータの値は仮想データの生成に使用した真のパラメータと多くの場合一致すると考えられる.

仮想データを生成するには,まず上記の式のパラメータを決めてやる必要がある.具体的には,回帰式のパラメータをまず設定する.

beta_0 = 0.2    # intercept
beta_1 = 0.3    # x_1の回帰係数
beta_2 = -.2    # x_2の回帰係数
beta_3 = 1      # x_3の回帰係数
beta_4 = 0      # x_4の回帰係数.x_4はyの原因ではないので0
gamma = 0.1     # 後述 (交絡)
sigma = 1       # 誤差項の標準偏差
parameters = c(beta_0, beta_1, beta_2, beta_3)

もし仮想データを生み出したデータ生成プロセスを回帰式でうまく推定できれば,回帰分析で推定した回帰係数と上で設定したbeta_* (真のパラメータ) は近い値になる (gammasigmaも設定しているが,ここでは省略する.また,beta_4は後ほど検討する).

次いで,上のパラメータを用いて,仮想データを生成しよう.ここでは,x_* (実際の解析では「原因の候補」) については次のものがあると考える.

これらの変数の種類は現実のデータ解析でも観察され得る. 特に問題になるのは,ある従属変数yの原因と考えられる (あるいはその影響を検証したい) 変数x_2が,第3の変数x_3からも影響を受ける場合である.言い換えれば,「x_2は確かにyの原因ではあるが,yの原因でもありx_2の原因でもある第3の変数x_3」が存在する場合 (交絡変数がある場合),解析には特に注意しなければならない.

この点については後ほど検討することにして,まずyと4つのxからなる仮想データを生成してしまおう.4つのxは,それぞれ上記のような特徴をもつことに注意すること.

set.seed(123) # 再現可能なように乱数のシードを固定する
n_sample <- 1000    # 標本数
x_1 = runif(n_sample, -10, 10)  # yの原因の1つの変数
x_3 = runif(n_sample, -20, 20)  # yとx_2の共通の原因である変数
x_2 <- gamma*x_3 + rnorm(n_sample)  # yの原因だが,x_3から影響を受ける変数
x_4 = runif(n_sample, -10, 10)  # yの原因「ではない」変数
mu <- beta_0 + beta_1*x_1 + beta_2*x_2 + beta_3*x_3 + beta_4*x_4 # 平均値:「yの原因」になっている変数はここで決める
y <- rnorm(n_sample, mu, sd = sigma)    # 正規分布
(sim_dat <- data_frame(y, x_1, x_2, x_3, x_4))
## # A tibble: 1,000 × 5
##              y        x_1        x_2         x_3       x_4
##          <dbl>      <dbl>      <dbl>       <dbl>     <dbl>
## 1  -10.4281716 -4.2484496 -1.9013078  -9.0550906 -3.910716
## 2    6.3917183  5.7661027 -0.6644873   3.7546773  6.656376
## 3  -14.3678123 -1.8204616 -1.3772410 -13.5926075  1.872950
## 4   15.8450208  7.6603481  1.2815458  14.1372096  6.143933
## 5   17.7584319  8.8093457 -1.1583861  13.9095663 -4.118984
## 6   -4.0772340 -9.0888700  0.9521207  -0.8845275 -7.178296
## 7   11.0229916  0.5621098  1.3444942  10.9476849  7.772422
## 8   -4.9300360  7.8483809  1.5978077  -8.1839968 -9.834300
## 9  -17.8613880  1.0287003 -1.0522893 -17.3748755  1.382413
## 10  -0.7067274 -0.8677053 -0.6848268  -2.3786750  9.351004
## # ... with 990 more rows

これで,sim_datという仮想データを生成できた.

2.2 仮想データの回帰分析

さて,上で生成した仮想データsim_datは我々が指定したデータ生成プロセス (data generation process) にしたがって生成されているのだから,回帰分析がうまく働けば,推定された回帰係数は「真のパラメータ」の値 (0.2, 0.3, -0.2, 1) に近くなるだろう. これを確かめるため,生成した仮想データを用いて回帰分析を行なってみる.なお,x_4はyの原因ではないことが分かっている (そうなるように仮想データを作っている) ので,ひとまず回帰式から除いておく (後ほど検討する).

fit <- lm(y ~ x_1 + x_2 + x_3, data = sim_dat)
display(fit)
## lm(formula = y ~ x_1 + x_2 + x_3, data = sim_dat)
##             coef.est coef.se
## (Intercept)  0.20     0.03  
## x_1          0.30     0.01  
## x_2         -0.16     0.03  
## x_3          0.99     0.00  
## ---
## n = 1000, k = 4
## residual sd = 1.00, R-Squared = 0.99

2.2.1 「真のパラメータ」と推定結果の比較

一見うまく推定できているように見える.確認のため,仮想データを生成した「真のパラメータ」と回帰分析で推定したパラメータ (回帰係数) を比較してみよう.

est_mat = cbind(parameters, summary(fit)$coef[,1:2],confint(fit))
colnames(est_mat) = c("true_beta", "est_beta", "se", "ci_lower", "ci_upper")
rownames(est_mat) = c("beta_0", "beta_1", "beta_2", "beta_3")
est_mat
##        true_beta   est_beta          se   ci_lower    ci_upper
## beta_0       0.2  0.1984271 0.031725450  0.1361707  0.26068349
## beta_1       0.3  0.2952745 0.005525376  0.2844318  0.30611722
## beta_2      -0.2 -0.1564711 0.031409737 -0.2181080 -0.09483425
## beta_3       1.0  0.9946395 0.004186310  0.9864245  1.00285446

ここでまとめた行列では,仮想データの生成に用いた真のパラメータを1列目に,回帰分析で推定したパラメータ (回帰係数) を2列目に,回帰係数の95%信頼区間を3列目と4列目に記録している.なお,回帰係数の信頼区間を求めるためにはstates::confint()関数を用いている.この関数は回帰分析の結果を引数とするので,ここではconfint(fit)と,線形回帰の結果を格納したfitオブジェクトを引数に指定している.confint()関数はデフォルトでは95%信頼区間を返すが,引数level0.950.90.99にしてやれば90%信頼区間や99%信頼区間が得られる.

なお,ここではコードを簡略化するためにconfint()関数を用いたが,回帰係数と標準誤差が分かっていれば信頼区間を自分で計算することもできる. 回帰係数を\(\beta_i\), その標準誤差を\(s_{i}\)とすると,\(\beta_{i}\)\(100\times(1-\alpha)\)%信頼区間は\(\beta_{i} \pm t_{n-k-1, 1-\alpha/2} \times s_{i}\)となる (\(n\)はサンプルサイズ,\(k\)は独立変数の数, \(t_{n-k-1, 1-\alpha/2}\)は自由度\(n-k-1\)\(t\)分布の\(1-\alpha/2\)%点).

自由度を引数として与えれば\(t\)分布 (に関連する値) を出力する関数はRに多数用意されているので,回帰係数・標準誤差・自由度が分かれば,次のコードでも95%信頼区間を求められる.

alpha = 0.05    # significance level
ci_mat = matrix(NA, nrow = 4, ncol = 2)
colnames(ci_mat) = c("ci_lower", "ci_upper")
deg_free = n_sample-4   # n-4 because the model includes intercept + 3 independent variables
t_thresh = qt(c(alpha/2, 1-alpha/2), df=deg_free)      # t-distribution
ci_mat[,1] = est_mat[,2] + t_thresh[1]*est_mat[,3] # 95% ci lower
ci_mat[,2] = est_mat[,2] + t_thresh[2]*est_mat[,3] # 95% ci upper

90%信頼区間や99%信頼区間を計算したければ,有意水準\(\alpha\)を指定しているalpha = 0.05を変更してやればよい (どんな数字に変更すればよいかは各自考える).以上のコードで出力された結果 (ci_matに格納されている) が,confint()関数で得た結果と等しいことを確認しておこう.

ci_mat
##        ci_lower    ci_upper
## [1,]  0.1361707  0.26068349
## [2,]  0.2844318  0.30611722
## [3,] -0.2181080 -0.09483425
## [4,]  0.9864245  1.00285446
confint(fit)
##                  2.5 %      97.5 %
## (Intercept)  0.1361707  0.26068349
## x_1          0.2844318  0.30611722
## x_2         -0.2181080 -0.09483425
## x_3          0.9864245  1.00285446

2.2.2 残る疑問

さて,回帰係数の推定値をみると,いずれのパラメータについても,回帰分析から得た推定値は真の値に近く,また信頼区間の中に真の値が入っているので,この回帰分析は各変数がyに与える影響を正しく捉えているといえそうだ.

しかし,いくつか疑問が残る.

  1. このシミュレーションで,たまたま「うまい結果」が出ただけかもしれない.この結果はどの程度信用できるのだろうか.
  2. x_2は確かにyの原因だが,x_3からも影響を受けている.x_3を回帰式から除くと,何か問題が出るのだろうか.
  3. この回帰式ではx_4を無視しているが,問題ないのだろうか.
  4. x_1は確かにyの原因だが,これを回帰式から除いたら何か問題が出るのだろうか.

以下では,シミュレーションを繰り返すことでこれら4つの疑問を検討しよう

2.3 シミュレーション (1):係数推定値と信頼区間

上の回帰分析は「うまくいった」ように見えるが,「たまたま」かもしれない.そこで,上と同じ要領で (1) n=1000のデータを生成し (2) 回帰分析をする,という作業を大量に繰り返そう.大量にこのシミュレーションを繰り返すことで,上で得た結果がどの程度信用できるのかを明らかにできる.

2.3.1 関数の作成とシミュレーションの設定

まずは,上で行なった作業をfor文で繰り返し,(1) と (2) の作業を3000回繰り返す.乱数を使ってデータを生成しているので,3000回のシミュレーションから3000個の異なる仮想データと解析結果が得られる.この作業から得た3000個の結果を使って,3000回中何回のシミュレーションにおいて回帰分析が真のパラメータを正しく推定できたかを検討しよう.

シミュレーションのため,まずデータを生成するための関数を定義する.

generate_sample <- function(
    n_sample = 1000,
    beta_0 = 0.2,
    beta_1 = 0.3,
    beta_2 = -0.2,
    beta_3 = 1,
    beta_4 = 0,
    gamma = 0.1,
    sigma = 1) {

    ## ----------------------------------------
    ## parameters
    betaz = c(beta_0, beta_1, beta_2, beta_3, beta_4)
    ## ----------------------------------------
    ## draw variables
    x_1 = runif(n_sample, -10, 10)
    x_3 = runif(n_sample, -20, 20)  ## confounder
    x_2 = gamma*x_3 + rnorm(n_sample, mean = 30, sd = 10)
    x_4 = runif(n_sample, -10, 10)
    ## ----------------------------------------
    ## set model and draw samples
    mu <- betaz[1] + betaz[2]*x_1 + betaz[3]*x_2 + betaz[4]*x_3 + betaz[5]*x_4
    y <- rnorm(n_sample, mu, sd = sigma)
    ## ----------------------------------------
    ## output
    sim_dat = data_frame(y, x_1, x_2, x_3, x_4)
    betaz = data_frame(parameter = c("intercept", str_c("x_", seq(1,4,by=1))), beta_true = betaz)
    return(list(betaz, sim_dat))
}

ここで定義したgenerate_sample()関数は,上の要領で生成した仮想データと使用したパラメータをリストオブジェクトとして返す.リストオブジェクトは,Rのオブジェクトを「なんでも」格納できる便利なオブジェクトで,Rで解析を行なう際に重宝する.リストのi番目の要素にアクセスするには,[[i]]とする.たとえば,sampleがリストオブジェクトなら,sample[[i]]とすればi番目の要素を取得できる.return(list(betaz, sim_dat))から分かる通り,この関数の出力は1番目の要素にパラメータのベクトルを,2番目の要素に生成した仮想データのtibbleオブジェクトを保持している.したがって,この関数の出力の2番目の要素にアクセスすれば,生成された仮想データを取得できる.

なお,この関数では新たなパラメータとしてbeta_4が追加されていることにも注意してほしい.ただし,beta_4=0とデフォルトの値を設定しているので,得られる結果は上のものと (乱数シードが同じなら) 同一になる.

このgenerate_sample()関数を定義してやることでコードを簡素化でき,またbeta_*のような (母数) パラメータも容易に変更できるようになった.もっとも,この関数は改善点が大量にあるが,ひとまず分かりやすさを優先して我慢する.

generate_sample()関数を利用して,3000回のシミュレーションを行なおう.

## ------------------------------------------------------------------------
set.seed(12345) ## set random seed
n_sample = 10^2*5   ## N samples
n_sim = 10^3*3  ## N simulations
## ------------------------------------------------------------------------
est_labz = c("beta", "se", "ci_lower", "ci_upper")
reg_mod = as.formula(y ~ x_1 + x_2 + x_3)
## ------------------------------------------------------------------------
## Loop
for (i in 1:n_sim) {

    ## ------------------------------------------------------------------------
    ## draw samples and store parameter values
    drawn_data = generate_sample(n_sample)
    para_dat = drawn_data[[1]]
    sim_dat = drawn_data[[2]]
    ## ------------------------------------------------------------------------
    ## linear regression
    fit <- lm(reg_mod, data = sim_dat)
    ## pull out estimates
    tmp_mat = cbind(summary(fit)$coef[,1:2],confint(fit))
    colnames(tmp_mat) = est_labz
    tmp_mat = tbl_df(as.data.frame(tmp_mat))
    tmp_mat$parameter = c("intercept", rownames(summary(fit)$coef)[2:nrow(summary(fit)$coef)])
    tmp_mat$round = i
    ## ------------------------------------------------------------------------
    ## store estimates
    if (i > 1) {
        est_mat = bind_rows(est_mat, tmp_mat)
    } else {
        est_mat = tmp_mat
    }
}   ## end for

1回のシミュレーションでは,2.2の例で得たものと同様に4行のアウトプットが得られる.これを3000回繰り返し,その結果を縦に結合している. したがって,このシミュレーションの繰り返しがうまくいっていれば,シミュレーション結果を保持するest_matは12,000行のオブジェクトになる.以下で確認するように,無事12,000行のオブジェクトが生成されている.

2.3.2 シミュレーション結果と信頼区間

意図した通り,12,000行のデータができているか確かめよう.ついでに,(1) 仮想データの生成に使った「真のパラメータ」と,(2) 回帰分析から得た係数推定値の95%信頼区間が「真のパラメータ」を含んでいるか否かを判別する変数を追加しておく (beta_truewithin95ci).

est_mat = est_mat %>%
    left_join(para_dat, by = "parameter")   %>%
    mutate(
        within95ci = ifelse( (beta_true >= ci_lower & beta_true <= ci_upper), TRUE, FALSE) ## 「真のパラメータ」が95%信頼区間に入っているか否か
    )   %>%
    select(
        round, parameter,
        contains("beta"),
        se:ci_upper, within95ci ## 列を並び替える
    )
est_mat
## # A tibble: 12,000 × 8
##    round parameter        beta beta_true          se    ci_lower
##    <int>     <chr>       <dbl>     <dbl>       <dbl>       <dbl>
## 1      1 intercept  0.30474997       0.2 0.143229491  0.02333864
## 2      1       x_1  0.28852734       0.3 0.007883426  0.27303831
## 3      1       x_2 -0.20472588      -0.2 0.004518395 -0.21360343
## 4      1       x_3  1.00784207       1.0 0.004118984  0.99974926
## 5      2 intercept  0.04319759       0.2 0.140616435 -0.23307972
## 6      2       x_1  0.28522879       0.3 0.007134705  0.27121082
## 7      2       x_2 -0.19582967      -0.2 0.004496793 -0.20466478
## 8      2       x_3  0.99460376       1.0 0.003641890  0.98744832
## 9      3 intercept  0.10285867       0.2 0.140386396 -0.17296667
## 10     3       x_1  0.29877057       0.3 0.007886414  0.28327567
## # ... with 11,990 more rows, and 2 more variables: ci_upper <dbl>,
## #   within95ci <lgl>
est_mat %>% filter(within95ci) %>% nrow
## [1] 11399

最後の行では,「回帰係数の95%信頼区間に,データを生成した『真のパラメータ』が入る」データのみを抽出している.12,000より少ない11399行が抽出されているので,「すべてのシミュレーションにおいて,すべての回帰係数が正しく推定されている」訳ではないようだ.

では,シミュレーションは全体として (平均として),どの程度真のパラメータをうまく推定しているのだろうか.図にして検討しよう. この回帰分析では,切片 (beta_0) を含めて4つのパラメータについての結果を得ている (x_4は今のところ省略している).上で指定している通り,「真の」パラメータの値はそれぞれ,0.2, 0.3, -0.2, 1だったことを思う出そう.まずは,beta_1について結果を確認する.

n_binz = 60
fig_height = 2
res_plot = ggplot(est_mat, aes(x = beta, y = ..density..)) +
    geom_histogram(color = "black", bins = n_binz)
# dev.new(width = fig_height*(1+sqrt(2)), height = fig_height)
res_plot %+% 
    filter(est_mat, parameter == est_mat$parameter[2]) +
    geom_vline(xintercept=est_mat$beta_true[2], color="green2") +
    xlab(str_c("Estimated coefficient of ", est_mat$parameter[2]))

このコードではまずggplot()関数をシミュレーション結果を格納したデータ全体を引数として実行し,%+% filter(est_mat, parameter == est_mat$parameter[2])で,beta_1変数の結果のみを抽出している.なお,ggplot2::geom_vline()関数は垂直の線を描く関数,ggplot2::xlab()は横軸のラベルを書く (変更する) 関数である.なお,est_mat$parameter[2]はx_1を返すので,その係数であるbeta_1がプロットされる.この説明を読んで「日本語でおk」と思った場合は,迷わずRのコンソールにest_matと打ち込み,parameter変数の2行目に入っている文字列を確認しよう

この図から分かるように,回帰係数の推定値は緑の線 (真の値) を中心に分布していることから,シミュレーション全体としてみたとき,線形回帰は平均して真のパラメータ (緑の線) をうまく推定しているといえる. この統計的性質を不偏性 (unbiasedness; 標本抽出を繰り返して得られる統計量の平均が真の値 [母数] に一致する性質) と呼ぶ.個々のシミュレーション単位で見ると,信頼区間に真のパラメータ値を含まない「ダメな」推定結果がある,という点に注意しよう. 言い換えると,線形回帰は常に正しい推定をする訳ではない.上のヒストグラムが示すように,あるデータセット (標本) が与えられ,それにを用いて回帰分析を行なったとき,真の値を過小評価することも,過大評価することもある

ここで,各々のシミュレーションでの仮想データは,完全に同一のデータ生成過程によって生成されていることに注意してほしい.つまり,背後に働く「原因」やメカニズムは全く同一にもかかわらず,あるときは真の値に近い推定値が得られ,別のときには真の値から遠い推定値が得られたということを意味する.信頼区間は,こうした推定に伴う不確実性を示す統計量である.すなわち,信頼区間は,同一のデータ生成過程によって生成されたデータセットに対し同一の解析を繰り返したとき,「得られた信頼区間のうちのどの程度が真のパラメータ値を区間内に含むか」という不確実性の程度を示す. 言い換えれば,95%信頼区間の意味は,「あるデータセットから推定した95%信頼区間の中に,真のパラメータの値が入る確率が95%」ということではないことに注意しよう (浅野・矢内 2013: 107-110).

ここでのシミュレーションでは3,000個の仮想データを生成しているが,現実のデータ解析では1個のデータしか得られない.したがって,自分の目の前にあるデータを解析したとき,正しい係数が推定されているとは限らない.このため,係数の推定値だけでなく,その不確実性を報告しなければならない.標準誤差や信頼区間は,こうした不確実性についての情報を教えてくれる.

2.3.3 課題

  1. est_matオブジェクトから,beta_0からbeta_3のいずれかについての結果を抽出する
  2. そのうち,「回帰係数の95%信頼区間に,データを生成した『真のパラメータ』が入る」データのみを抽出し,行数を数える
    1. の行数を3,000で割ると,おおよそ何%になるだろうか.また,信頼区間の意味を踏まえると,その数字は何を意味するのだろうか.

2.3.4 作図の自動化

同じ要領で,beta_0, beta_2, beta_3の推定値の分布も図示しよう.beta_1についてプロットしたときはres_plot %+% filter(est_mat, parameter == est_mat$parameter[2]) + geom_vline(xintercept=est_mat$beta_true[2], color="green2")とした. このコードのbeta_1を抽出しているparameter == est_mat$parameter[2]という部分と,緑の線の位置を指定しているxintercept=est_mat$beta_true[2]という部分を変更してやればよい. つまり,[2]の部分を[1], [3], [4]としてやれば,それぞれbeta_0, beta_2, beta_3の推定値がプロットできる.数字以外同じコードを何度も書くのは面倒なので,for文で処理する.

## Plot beta_0, beta_2, and beta_3
plot_loop = c(1, 3, 4)
fig_height = 2
for (i in plot_loop) {
    current_plot = res_plot %+% 
        filter(est_mat, parameter == est_mat$parameter[i]) +
        geom_vline(xintercept=est_mat$beta_true[i], color="green2") +
        xlab(str_c("Estimated coefficient of ", est_mat$parameter[i]))
    # dev.new(width = fig_height*(1+sqrt(2)), height = fig_height)
    print(current_plot)
}

beta_1と同様に,beta_0, beta_2, beta_3ともに真の値を中心に推定値が分布していることから,平均としてうまく推定できたことが分かる. なお,上のコードをRで実行する際にはdev.new()の先頭にある#を削除すること (削除しないと,beta_3のヒストグラムのみが残ってしまう).

2.4 シミュレーション (2):交絡変数と欠落変数バイアス

2.4.1 ここまでのシミュレーションでの想定とその変更

さて,上のシミュレーションではデータ生成過程を捉えた「正しいモデル」による分析を行なっていた.しかし,現実のデータを解析する際には,「正しいモデル」は分からない (だから解析する).では,「本当は正しいモデルに入っている変数を無視/追加してしまったとき」(「正しくない」モデルで推定してしまったとき) どのような問題が生じるのだろうか

この問題を考えるにあたって,ここまでのシミュレーションは,それぞれのx_*を次のように設定していたことをまず思い出そう.

  • yの原因であり,他のxとは無関係の変数:x_1
  • yの原因であり,他の変数x_3からも影響を受ける変数:x_2
  • yの原因であり,他の変数x_2にも影響を与える変数:x_3
  • yの原因ではない変数:x_4

ここで,x_3のように他の (yの原因である) xとyの両方に影響する変数を,交絡変数 (confounder) と呼ぶ. 「正しいモデル」は,x_1, x_2, x_3のすべてを含んでいた.

では,交絡変数であるx_3を回帰モデルから外すと,その影響を受けているx_2の係数の推定結果は,どのような影響を受けるだろうか.逆に,x_4を回帰モデルに加えると,どのような影響が生じるだろうか.また,一見x_2と無関係そうに見えるx_1を外した場合,x_2の係数推定値は影響を受けるのだろうか. これらの疑問に回答するシミュレーションを行なえば,上で提示した

  1. x_2は確かにyの原因だが,x_3からも影響を受けている.x_3を回帰式から除くと,何か問題が出るのだろうか.
  2. この回帰式ではx_4を無視しているが,問題ないのだろうか.
  3. x_1は確かにyの原因だが,これを回帰式から除いたら何か問題が出るのだろうか.

という疑問に回答できる.

常に回帰式に投入しているx_2とその係数beta_2に着目して回答を先に言っておくと,

  1. x_3を回帰式で無視した場合:beta_2の推定値について,平均して誤った推定結果が出てしまう (欠落変数バイアス)
  2. x_4を回帰式に投入した場合:beta_2の推定値に強い影響はない
  3. x_1を回帰式で無視した場合:beta_2の推定値に強い影響はない

2.のケースのように,交絡変数を無視したことで推定結果に生じるバイアスを,欠落変数バイアス (omitted variable bias)と呼ぶ. 実際の解析でもこのバイアスが生じないよう特に注意しなければならず,交絡変数は必ず回帰式に投入するなどして,必ずその影響をコントロールしなければならない.

2.4.2 例:回帰式にx_4を加えた場合

以上の3つのケースをそれぞれシミュレーションで確かめてみよう.実は,上で行った3000回のシミュレーションと作図のコードはある程度一般的に書いてあるので,1行変えてやるだけで上記3つのケースについてのシミュレーションを実行できる.具体的には,回帰式を指定している

reg_mod = as.formula(y ~ x_1 + x_2 + x_3)

という行を変更してやればよい.

以下の例では,この回帰式を指定する行を書き換えて2のケース (yの原因ではないx_4を追加した回帰式) を実行している.beta_4 (x_4の回帰係数) の推定値は真の値である0を中心に分布し,他の回帰係数も真の値を中心に正しく推定されていることが分かる.ただし,以下のコードでは切片 (intercept) の推定値のヒストグラムは省略している (切片の値そのものに興味があることは希なため).

set.seed(12345) ## set random seed
n_sample = 10^2*5   ## N samples
n_sim = 10^3*3  ## N simulations
## ------------------------------------------------------------------------
est_labz = c("beta", "se", "ci_lower", "ci_upper")
reg_mod = as.formula(y ~ x_1 + x_2 + x_4)
## ------------------------------------------------------------------------
## Loop
for (i in 1:n_sim) {

    ## ------------------------------------------------------------------------
    ## draw samples and store parameter values
    drawn_data = generate_sample(n_sample)
    para_dat = drawn_data[[1]]
    sim_dat = drawn_data[[2]]
    ## ------------------------------------------------------------------------
    ## linear regression
    fit <- lm(reg_mod, data = sim_dat)
    ## ------------------------------------------------------------------------
    ## store estimates
    tmp_mat = cbind(summary(fit)$coef[,1:2],confint(fit))   ## pull out estimates
    colnames(tmp_mat) = est_labz
    tmp_mat = tbl_df(as.data.frame(tmp_mat))
    tmp_mat$parameter = c("intercept", rownames(summary(fit)$coef)[2:nrow(summary(fit)$coef)])
    tmp_mat$round = i
    if (i > 1) {
        est_mat = bind_rows(est_mat, tmp_mat)
    } else {
        est_mat = tmp_mat
    }
}   ## end for

## Modification
est_mat = est_mat %>%
    left_join(para_dat, by = "parameter")   %>%
    mutate(
        within95ci = ifelse( (beta_true >= ci_lower & beta_true <= ci_upper), TRUE, FALSE)
        )   %>%
    select(
        round, parameter,
        contains("beta"),
        se:ci_upper, within95ci 
        )

## Plot
n_binz = 60
fig_height = 2
res_plot = ggplot(est_mat, aes(x = beta, y = ..density..)) + geom_histogram(color = "black", bins = n_binz)
plot_loop = 2:length(unique(est_mat$parameter))
for (i in plot_loop) {
    current_plot = res_plot %+% 
        filter(est_mat, parameter == est_mat$parameter[i]) +
        geom_vline(xintercept=est_mat$beta_true[i], color="green2") +
        xlab(str_c("Estimated coefficient of ", est_mat$parameter[i]))
    # dev.new(width = fig_height*(1+sqrt(2)), height = fig_height)
    print(current_plot)
}

2.4.3 課題:回帰式からx_1を除いたシミュレーションとx_3を除いたシミュレーション

回帰式をどのように変更すれば1. と3. のケースを実行できるか,またx_2の回帰係数はどのような影響を受ける (受けない) か,各自考えて確かめてみよう. なお,上のコードをRで実行する際にはdev.new()の先頭にある#を削除すること.