Chapter 12 重回帰分析
12.1 重回帰分析とは
重回帰分析は、単純に言えば説明変数が２つ以上ある回帰式を用いる回帰分析である。
しかし、重回帰分析は、単純に多くの変数を入れることで、目的変数を説明できる部分が増える、というような目的で用いられるのではない。
特に、実証分析における因果推定では、重回帰分析における係数が、他の説明変数の影響を除いた影響として推定されることが重要である。
12.2 回帰分析における真のモデル
回帰分析では、そのデータを生み出している元となる存在があると考える(専門用語ではData Generating Processと呼ぶ)。
これは、現象を起こすメカニズムが存在し、そのメカニズムは神のみぞ知ると考えていると言ってもよい。
私たち分析者は、神のみぞ知るメカニズムに迫ろうとしているのだ。
実際の社会で起こっているメカニズムの本当の姿は（真のモデル）は、神のみぞ知るが、シミュレーションによる仮想世界を作ってその世界の真のメカニズムを知りながら、分析をすることはできる。いわば、仮想世界では神になることができるのだ。
12.2.1 シミュレーション
12.2.1.2 データの生成
今から行うのは、以下のような真のモデルがあると仮定して、データを生成することである。
ここでは、あるスーパーマーケットチェーンのお店ごとのデータを生成するとしよう。このスーパーマーケットチェーンでは全国で300店舗展開しているとする。 300店舗のうち、一部の店舗では広告を出稿することにした。ここで分析しようとしているのは、広告を出すことによって平均的に売上が上がるのかどうかである。
当然ながら売上は広告だけでは決まらない。ここでは、「神」となってデータを生成してみよう。 各店舗\(i\)の売上は、以下のモデルで決まるとしよう。これが「真のモデル」である。
\[ Sales_i = \alpha + \beta_1 Ads_i + \beta_2 Age_i + \beta_3 Income_i + \varepsilon_i \]
\(Sales_i\)は店舗\(i\)の売上、\(Ads_i\)は広告を出したかというダミー変数で出向していれば1, していなければ0を取る。\(Age_i\)は店長の年齢である。\(Income_i\)はその店舗が出店している地域の平均所得だとしよう。
単純なモデルであるが売上は、広告、店長の年齢、地域の平均所得で決まるとしよう。
しかし、広告を出すかどうか、という意思決定は外生的に決まらない。外生的とは、モデルの外で天から降ってきたように決まることをいう。しかし、この広告を出すかどうかは内生的に決まるものだとしよう。モデルの中の変数によって決まるとする。
\[ Ads_i = I\{100 \leq \gamma_1 Age_i + \gamma_2 Income_i + u_i\} \]
すこしややこしい数式だが、この数式が言っているのは次のとおりだ。店長の年齢と地域の平均所得で決まる数値が100より大きい店舗では広告を出稿する。
そして、仮想世界の神であるので、このモデルのパラメーターの数値を決めることができる。 以下のようにパラメーターの数値を決めた。
そして、いよいよデータの生成である。 データの生成には確率的な要素も入れ込むが、再現性を担保するためにシード値を決めておく。
そして、以下のようにデータを生成する。
少し説明が難しい要素もあるので、ここはコピー&ペーストでよいので自分のスクリプトに貼り付けて実行してみよう。
実行する際は、かならず上のset.seed(3)を実行してから実行すること。
# シミュレーションデータの作成
dat_sim = tibble(
shop_ID = 1:N2,
age = round(runif(N2,25,60)),
income = round(runif(N2,15,100))*10,
pop = round(runif(N2,10,1000))*100,
dogs = round(runif(N2,0,0.5),digits=2),
u_ads = rnorm(N2,0,sd=10),
z_ads = round(g1*age+g2*income,digits=1),
y_ads = z_ads + u_ads,
ads = ifelse(100 <= y_ads,1,0),
epsilon = rnorm(N2,sd=20),
sales = round(a1 + b1*ads + b2*age + b3*income + epsilon,digits=1)
) %>%
select(!contains("_ads") & !starts_with("u")) %>%
relocate(ads, .before=age)ここでは、上のモデルで使った変数に加えて、popという地域の人口に当たる変数と、dogsという地域の犬を飼っている割合という変数も生成している。
データの記述統計を見てみよう。
| Name | dat_sim |
| Number of rows | 300 |
| Number of columns | 8 |
| _______________________ | |
| Column type frequency: | |
| numeric | 8 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| shop_ID | 0 | 1 | 150.50 | 86.75 | 1.00 | 75.75 | 150.50 | 225.25 | 300.00 | ▇▇▇▇▇ |
| ads | 0 | 1 | 0.53 | 0.50 | 0.00 | 0.00 | 1.00 | 1.00 | 1.00 | ▇▁▁▁▇ |
| age | 0 | 1 | 42.11 | 9.81 | 25.00 | 34.00 | 42.00 | 51.00 | 60.00 | ▇▇▇▇▆ |
| income | 0 | 1 | 603.23 | 251.09 | 150.00 | 397.50 | 625.00 | 832.50 | 1000.00 | ▆▇▆▆▇ |
| pop | 0 | 1 | 49726.00 | 28658.81 | 1200.00 | 26000.00 | 48950.00 | 74075.00 | 100000.00 | ▇▇▇▇▇ |
| dogs | 0 | 1 | 0.24 | 0.14 | 0.00 | 0.13 | 0.25 | 0.35 | 0.50 | ▇▇▇▆▆ |
| epsilon | 0 | 1 | -1.28 | 20.71 | -64.81 | -16.14 | -0.31 | 12.57 | 51.91 | ▁▅▇▇▁ |
| sales | 0 | 1 | 789.90 | 73.83 | 624.80 | 735.03 | 789.45 | 846.12 | 961.00 | ▃▆▇▆▂ |
12.2.2 回帰モデルの推定
12.2.2.1 新たなパッケージの導入
ここで新しいパッケージを導入する。
fixestというパッケージをインストールしよう。
インストールするには、RStudioの右下ペーンのPackagesタブからInstallをクリックして名前を指定するか、以下のコードを実行する。
そして、忘れずにパッケージをロードする。
このパッケージは、回帰分析を行うパッケージである。
前章までは回帰分析をRにデフォルトで入っているlm()関数で行ってきた。
このパッケージに入っているfeols()関数は、基本的に同じように回帰分析を行うことができるが、より柔軟にいろんなことができるためこれを使っていく。
また、このパッケージに入っているetable()という関数を使うと、回帰分析の結果をよりわかりやすく表示することができる。
12.2.2.2 モデルの推定
ここでの疑問は「広告の出稿は売上増加に貢献するか」であった。
普通のデータなら、真のモデルはわからないが、ここでは我々は神なので真のモデルがわかっている。
まず、このデータから真のモデルと同じモデルを推定してみよう。
## model_true
## Dependent Var.: 支店の売上
##
## Constant 499.6*** (6.860)
## 広告ダミー 49.66*** (4.028)
## 店長の年齢(歳) 4.985*** (0.1299)
## 地域の平均所得(万円) 0.0899*** (0.0078)
## _______________ __________________
## S.E. type IID
## Observations 300
## R2 0.92132
## Adj. R2 0.92052
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1設定したパラメーターに近い数値が推定されている。
真のモデルを知っていれば広告の効果は50であるとわかる。 また、真のモデルに基づいたデータ分析でも推定値は49.66と非常に近い値が推定されている。
しかし、本当は分析者は真のモデルは知らない。 一般的には、単回帰分析から始めるが、単回帰分析をしてみるとどうなるだろうか。
\[ Sales_i = \alpha + \beta_1 Ads_i + \varepsilon_i \]
練習問題: 説明変数を広告とした単回帰モデルを推定し、結果を表示させなさい
## model_1
## Dependent Var.: 支店の売上
##
## Constant 733.4*** (4.261)
## 広告ダミー 107.3*** (5.872)
## _______________ ________________
## S.E. type IID
## Observations 300
## R2 0.52853
## Adj. R2 0.52694
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1推定値は107と出た。これは真の値である50に比べると２倍以上である。 つまり、このモデルでは広告の効果を２倍以上見積もっていることになる。
何度も言うように、本来ならば分析者は真のモデルや真の値を知らない。 そのため、単純に単回帰モデルを分析して、結果を解釈すると「広告の出稿は107万円の売上を増加させる」という間違った結果を得てしまう。
次に、分析者は他の変数を含めてみようと思うかもしれない。 地域の人口が売上に影響を与えていそうだ、という思って次のような重回帰モデルを推定するとしよう。
\[ Sales_i = \alpha + \beta_1 Ads_i + \beta_4 Pop_i + \varepsilon_i \]
練習問題: 説明変数を広告と人口(pop)とした重回帰モデルを推定し、結果を表示させなさい
## model_2
## Dependent Var.: 支店の売上
##
## Constant 737.1*** (6.782)
## 広告ダミー 107.1*** (5.885)
## pop -7.32e-5 (0.0001)
## _______________ _________________
## S.E. type IID
## Observations 300
## R2 0.52933
## Adj. R2 0.52616
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1結果は、popの係数の数値は非常に小さく、また統計的にゼロである可能性を棄却できなかった。 また、広告の効果も単回帰とは変わらない。
次に、分析者は店長の年齢は経験値などに比例しているので売上に関係があると思ったとする。
そこで、次に店長の年齢ageを含めたモデルを推定する。
練習問題: 説明変数を広告、人口(pop)、店長の年齢(age)とした重回帰モデルを推定し、結果を表示させなさい
## model_3
## Dependent Var.: 支店の売上
##
## Constant 552.0*** (6.929)
## 広告ダミー 86.14*** (2.977)
## pop -4.28e-5 (5.06e-5)
## 店長の年齢(歳) 4.623*** (0.1516)
## _______________ __________________
## S.E. type IID
## Observations 300
## R2 0.88635
## Adj. R2 0.88519
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1models_ads = list()
models_ads[["Model 1"]] = feols(sales ~ ads, se="iid", data=dat_sim)
models_ads[["Model 2"]] = feols(sales ~ ads + age, se="iid", data=dat_sim)
models_ads[["Model 3"]] = feols(sales ~ ads + age + income, se="iid", data=dat_sim)
models_ads[["Model 4"]] = feols(sales ~ ads + age + income + pop, se="iid", data=dat_sim)
models_ads[["Model 5"]] = feols(sales ~ ads + age + income + pop + dogs, se="iid", data=dat_sim)
etable(models_ads)## Model 1 Model 2 Model 3
## Dependent Var.: 支店の売上 支店の売上 支店の売上
##
## Constant 733.4*** (4.261) 549.7*** (6.372) 499.6*** (6.860)
## 広告ダミー 107.3*** (5.872) 86.26*** (2.973) 49.66*** (4.028)
## 店長の年齢(歳) 4.625*** (0.1515) 4.985*** (0.1299)
## 地域の平均所得(万円) 0.0899*** (0.0078)
## pop
## dogs
## _______________ ________________ _________________ __________________
## S.E. type IID IID IID
## Observations 300 300 300
## R2 0.52853 0.88607 0.92132
## Adj. R2 0.52694 0.88530 0.92052
##
## Model 4 Model 5
## Dependent Var.: 支店の売上 支店の売上
##
## Constant 501.6*** (7.241) 503.3*** (7.492)
## 広告ダミー 49.59*** (4.030) 49.58*** (4.031)
## 店長の年齢(歳) 4.982*** (0.1300) 4.982*** (0.1301)
## 地域の平均所得(万円) 0.0898*** (0.0078) 0.0901*** (0.0078)
## pop -3.7e-5 (4.21e-5) -3.67e-5 (4.21e-5)
## dogs -7.638 (8.689)
## _______________ __________________ __________________
## S.E. type IID IID
## Observations 300 300
## R2 0.92152 0.92173
## Adj. R2 0.92046 0.92040
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 112.3 モデルの当てはまり
12.3.1 決定係数
最小二乗法では残差二乗和が最小になる。
しかし、ゼロではない。これは、回帰モデルですべての観測値を完全に予測できているわけではないことを示している。
一般的には完全に予測することは不可能である。
この最小二乗法によって求めたモデルがどの程度データの変動を説明してくれているかを測る指標が決定係数。
決定係数は残差変動と総変動の比を１から引いたものとして定義される。
残差変動とは、モデルで説明できていないデータの変動である。つまりこれは残差二乗和である。
\[ \sum_{i=1}^{n}\hat{u}^2_{i} \]
- 総変動とはデータの目的変数のすべてのばらつきである。つまり、目的変数の平均とすべての観測値の差の二乗和である。
\[ \sum_{i=1}^{n} (y_{i} - \bar{y}_{i})^2 \] - 決定係数はすべての変動のうち、残差変動ではない部分の割合はどれぐらいか？ すなわち、モデルで説明される変動はどれぐらいか？を測っている。
\[ R^2 = 1 - \frac{\sum_{i=1}^{n}\hat{u}^2_{i}}{\sum_{i=1}^{n} (y_{i} - \bar{y}_{i})^2} \]
この数値は、高ければ高いほど、データを説明できている、ということができる。
じゃあ、この数値が高いモデルが「良いモデル」なのかというと、そうとは限らない。 その理由は２つある。 一つは、我々の目的はあくまで因果関係を実証することになる。決定係数は予測が当てはまっているか、というものを測るに過ぎない。 データ上、決定係数が高いが、不要な変数などを含めてしまっている可能性も存在し、それが因果関係を示すパラメータに影響を与えていたら、決定係数が高くても意味がない。
もう一つは、決定係数は説明変数が多ければ多いほど、高くなる傾向にある。 つまり、追加した説明変数に対して説明力（予測を改善する力）がなくても、見かけ上決定係数が上昇するのである。 その欠点を改善したものが次の自由度調整済み決定係数である。