Chapter 12 回帰分析
【Section 12.9に関連動画を紹介しています。】
回帰分析とは、従属変数と独立変数の関係を数式（モデル）で表し、そのパラメータを推定する分析方法です。 ここでは、もっとも基本的な回帰分析である線形回帰（Linear Regression）を扱います。
12.1 線形回帰とは
線形回帰では、従属変数 \(y\) と説明変数 \(x\) があるとき、\(y\) と \(x\) の関係の以下の式（回帰式）で表します。
別の言い方をすると、 \(y\) を \(x\) の１次関数で表すと言うことです。 傾きを示す \(\beta\) は \(x\) と \(y\) の関係を示す重要な数値で、係数（Coefficient）と呼ばれます。 これがプラスであれば \(x\) と \(y\) に正の関係があることになりますし、マイナスであれば負の関係、0であれば関係がないことになります。 \(\beta\) は \(x\) が1単位高いとき、\(y\) が \(\beta\) 単位高い、という意味になります。 回帰分析では変数の単位を常に把握して分析するよう心がけましょう。
\(\alpha\) や \(\beta\) は最小二乗法（OLS）という方法で推定します。 この文書ではRでの操作方法に特化して説明するので、詳しくは各自調べてください。
12.2 Rでの線形回帰
Rでは回帰式をy ~ xと表記します。（ここはChapter 11で学んだとt検定同じ表記です）
線形回帰を行う関数がlm()で（lmはlinear modelの略）、以下のように表記します。
lm(回帰式, data = データフレーム)例として、従属変数を年収inc、独立変数を教育年数educとした線形回帰分析を行ってみましょう。
回帰式は \(inc = \alpha + \beta educ\) です。
lm(inc ~ educ, data = saving)##
## Call:
## lm(formula = inc ~ educ, data = saving)
##
## Coefficients:
## (Intercept) educ
## 1342.7 742.5Coefficientsのinterceptが\(\alpha\)、educが\(\beta\)を表します。 \(\beta\)を見ると、教育年数が1年高いと年収が742.5ドル高くなっていることがわかります。
lm()を実施すると、以上のように係数の推定値を得ることができますが、統計的推測のための数値（t値、p値）がわかりません。
lm()の結果をsummary()に渡すと、詳細を見ることができます。
lm(inc ~ educ, data = saving) %>%
summary()##
## Call:
## lm(formula = inc ~ educ, data = saving)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7570 -3297 -1288 1617 20743
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1342.7 1763.5 0.761 0.448
## educ 742.5 146.1 5.084 1.78e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4993 on 98 degrees of freedom
## Multiple R-squared: 0.2087, Adjusted R-squared: 0.2006
## F-statistic: 25.84 on 1 and 98 DF, p-value: 1.777e-06Coefficientsで以下の数値が表示されるようになりました。 それぞれの意味は教科書を確認してください。
Estimate: 係数Std. Error: 標準誤差t value: t値Pr(>|t|): p値
educのp値は1.78e-06と表記されています。
これは、\(1.78 \times 10^{-6} = 1.78 \times 0.000001 = 0.00000178\)を表します。
とても小さい値なので、educの係数 \(\beta\) は有意に負と言えます。
他にも、決定係数 \(R^2\) などの数値も記載されています。
summary()で詳細を見ることができますが、やや見づらい印象があります。
表形式にするため、Chapter 11と同じくbroomパッケージの関数tidy()を使ってみましょう。
lm(inc ~ educ, data = saving) %>%
tidy()## # A tibble: 2 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 1343. 1764. 0.761 0.448
## 2 educ 743. 146. 5.08 0.00000178表形式で見やすくなりました。
★練習問題
- 収入
incを世帯人数sizeに回帰し、係数の意味を解釈してください。
12.3 重回帰分析
回帰分析には複数の説明変数を含め、それぞれの説明変数と被説明変数との変数を検証することができます。 複数の説明変数がある回帰分析を重回帰分析、一方上で説明した１つの説明変数がある回帰分析を単回帰分析と呼びます。 回帰式は説明変数を \(x_1, x_2\) としたとき、以下のようになります。
Rでの回帰式の表記は、y ~ x1 + x2と、説明変数を+でつなげます。
ここでは、収入incを教育年数educと世帯人数sizeに回帰してみましょう。
lm(inc ~ educ + size, data = saving) %>%
tidy()## # A tibble: 3 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 3027. 2283. 1.33 0.188
## 2 educ 743. 146. 5.10 0.00000171
## 3 size -389. 335. -1.16 0.249回帰式の解釈は単回帰分析と同様です。 ここでは、教育年数の結果は先ほどと大きく変わりませんでした。
説明変数はさらに3つ、4つ…と増やすことができ、Rではy ~ x1 + x2 + x3 + x4 + ...と+でつなげていきます。
説明変数が多く回帰式が長くなる場合は、回帰式を一旦オブジェクトとして保存しておくとコードが読みやすくなります。
例えば、上記の回帰分析にさらに年齢age、黒人ダミーblackを加えてみましょう。
equation <- inc ~ educ + size + age + black
lm(equation, data = saving) %>%
tidy()## # A tibble: 5 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -10005. 3934. -2.54 0.0126
## 2 educ 857. 144. 5.97 0.0000000408
## 3 size -101. 320. -0.317 0.752
## 4 age 271. 66.3 4.09 0.0000917
## 5 black -553. 1878. -0.294 0.7694つの説明変数による回帰分析の結果が出力できました。黒人ダミーの解釈については、次項をご覧ください。
★練習問題
- 貯蓄額
savを被説明変数、教育年数educ、世帯人数size、年齢ageを説明変数とする回帰分析を行い、それぞれの係数を解釈しなさい。
12.4 ダミー変数の利用
条件を満たすグループに1、条件を満たさないグループに0をとるダミー変数は、そのまま回帰分析に組み込むことができます。 ただし、解釈には注意が必要です。 ダミー変数の係数の解釈は、「条件を満たすグループは満たさないグループと比べて、被説明変数が（係数）大きい」となります。
ここでは、収入incを黒人ダミーblackに回帰してみましょう。
lm(inc ~ black, data = saving)##
## Call:
## lm(formula = inc ~ black, data = saving)
##
## Coefficients:
## (Intercept) black
## 10181 -3418黒人ダミーの係数は-3418です。これは、「黒人はそれ以外と比べて、年収が3418ドル低い」ことを意味します。
12.5 カテゴリ変数の利用
カテゴリ変数はそのまま回帰分析に組み込むことができません。 代わりにカテゴリの数-1個のダミー変数を作成し、回帰分析に組み込みます。 例えば、A, B, Cという3つのカテゴリがある場合、
- Bであれば1、そうでなければ0をとる「Bダミー」
- Cであれば1、そうでなければ0をとる「Cダミー」
を作成し、回帰分析の独立変数として両方入れます。 係数の解釈には注意が必要です。それぞれの係数の解釈は、ダミー変数として入っていない「参照レベル」との比較となります。
- Bダミーの係数は「カテゴリBはカテゴリAと比べて被説明変数が（係数）大きい」
- Cダミーの係数は「カテゴリCはカテゴリAと比べて被説明変数が（係数）大きい」
カテゴリBとカテゴリCを比べたい場合は、双方の係数を比較してください。
Rではカテゴリ変数がfactor型になっていれば、自動的にダミー変数を作成して回帰分析に入れてくれます。（もちろん、自分で作成しても構いませんし、そのほうが良い場合もあります）
ここでは、年齢のカテゴリ変数age_categoryを作成し、年収incに対する回帰分析に組み込んでみましょう。
saving_with_age_category <-
saving %>%
mutate(age_category = case_when(age < 30 ~ "20s",
age >= 30 & age < 40 ~ "30s",
age >= 40 & age < 50 ~ "40s",
age >= 50 ~ "50s"
)
)
lm(inc ~ age_category, data = saving_with_age_category)##
## Call:
## lm(formula = inc ~ age_category, data = saving_with_age_category)
##
## Coefficients:
## (Intercept) age_category30s age_category40s age_category50s
## 7685 1330 3761 3885ここでは自動的に“20s”が参照レベルとなり、回帰分析には含まれていません。 各係数は20代と比べて収入がどれくらい高いかを表します。
- 30代は20代より年収が1330ドル高い
- 40代は20代より年収が3761ドル高い
- 50代は20代より年収が3885ドル高い
それぞれ比較すると、年収は20代<30代<40代<50代と年齢が上がるにつれ高くなることがわかります。 しかし、40代と50代では差はほぼありません。
12.6 交互作用の導入
ある説明変数と被説明変数のとの関係が別の説明変数によって変化するような状況のことを交互作用がある、と言います。
例えば、収入と教育年数の関係が黒人ダミーによって異なる、といった場合です。
Rの回帰式では、交互作用は2つの変数を:でつなげることで表現できます。
上記の例をRで実施してみましょう。
lm(inc ~ educ + black + educ:black, data = saving) %>%
tidy()## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 1595. 1926. 0.828 0.410
## 2 educ 727. 157. 4.63 0.0000115
## 3 black -525. 5773. -0.0909 0.928
## 4 educ:black -63.1 615. -0.103 0.918交互作用の解釈は少し複雑になります。
まずここでは教育年数educの係数が727なので、教育年数1年高いと収入が727ドル高いことになります。
これは参照レベルである「黒人以外」での教育年数と収入の関係になります。
交互作用educ:blackは-63ですが、これは黒人での教育年数と収入の係数が727-63=664であることを意味します。
なので、黒人では教育年数が1年高いと収入が664ドル高いということになります。
ただし、ここでは交互作用のp値が大きく有意ではありません。
すなわち、「黒人とそれ以外では教育年数と収入の関係が同じ」という仮説は棄却できていません。
交互作用の分析の際は、交互作用に使用した変数単体も必ず回帰式に含めるようにしてください。
*を使用すると、単体・相互作用をすべて自動的に作成・分析してくれます。
lm(inc ~ educ*black, data = saving) %>%
tidy()## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 1595. 1926. 0.828 0.410
## 2 educ 727. 157. 4.63 0.0000115
## 3 black -525. 5773. -0.0909 0.928
## 4 educ:black -63.1 615. -0.103 0.918前の分析と全く同じ結果となっています。
★練習問題
- 収入
incと年齢ageの関係が黒人ダミーによって異なるかどうか、交互作用を含めた回帰分析を行い、その結果を解釈しなさい。
12.7 回帰分析の表をまとめる（modelsummary）
ここまでは、回帰分析の結果はtidy()を用いて整理し、卒業論文に貼り付ける想定で話を進めてきました。
回帰分析の数が少ない場合は、この方法でも十分でしょう。
しかし、実施した回帰分析が多い場合、同じような表が続き冗長になります。
ここでは、modelsummaryパッケージを用いて、複数の回帰分析の結果をまとめた表を作成してみる方法を説明します。
modelsummaryの説明には、以下のページを参考にさせていただきました。詳しい情報を知りたい方は、こちらもご覧ください。
https://keita43a.hatenablog.com/entry/2020/05/29/210250
まずはパッケージをインストールし、呼び出しておきましょう。以降の操作でkableExtraやgtといった他のパッケージのインストールを要求されたら同様にインストールしておきましょう。
install.packages("modelsummary")
library(modelsummary)１つの回帰分析を表にする方法から説明します。
まず被説明変数を収入inc、説明変数を教育年数educとする単回帰分析の結果をオブジェクト（ここではreg1）として保存します。
そのオブジェクトを引数として、msummary()を実行すると、表が出力されます。
reg1 <- lm(inc ~ educ, data = saving)
msummary(reg1)| Model 1 | |
|---|---|
| (Intercept) | 1342.745 |
| (1763.546) | |
| educ | 742.530 |
| (146.062) | |
| Num.Obs. | 100 |
| R2 | 0.209 |
| R2 Adj. | 0.201 |
| AIC | 1990.9 |
| BIC | 1998.7 |
| Log.Lik. | −992.455 |
| F | 25.844 |
デフォルトで表示されるのは係数と標準誤差（かっこ内）です。
統計量としては、Num.Obs.が観測数、R2が決定係数、R2 Adj.が調整済み決定係数を表します。（その他の説明は省略）
表はさまざまな形式で出力できます。例えばWord形式で保存したい場合は、以下のようにオプションを加えます。
msummary(reg1, 'result.docx')デフォルトではさまざまな統計量が表示されますが、論文に書くにはやや冗長です。 gof_omitオプションを使ってAIC, BIC, Log Likelihood, F値を落としてみましょう。
msummary(reg1, gof_omit = 'AIC|BIC|Log.Lik.|F')| Model 1 | |
|---|---|
| (Intercept) | 1342.745 |
| (1763.546) | |
| educ | 742.530 |
| (146.062) | |
| Num.Obs. | 100 |
| R2 | 0.209 |
| R2 Adj. | 0.201 |
statisticオプションを使うと、統計的推測に使用する統計量（デフォルトは標準誤差）を指定できます。
たとえば95%信頼区間を使用したい場合は以下のように指定します。
msummary(reg1, statistic = 'conf.int', conf_level = .95, gof_omit = 'AIC|BIC|Log.Lik.|F')| Model 1 | |
|---|---|
| (Intercept) | 1342.745 |
| [−2156.954, 4842.445] | |
| educ | 742.530 |
| [452.674, 1032.385] | |
| Num.Obs. | 100 |
| R2 | 0.209 |
| R2 Adj. | 0.201 |
今はあまりおすすめしませんが、有意水準を示すアスタリスク（*）を示したい場合は、estimateオプションをestimate = "{estimate}{stars}"と書くと表示できます。
２つ以上の回帰分析の結果をまとめるには、結果をlist形式にまとめて保存し、そのオブジェクトを引数に指定します。
被説明変数を収入inc、説明変数を教育年数educと年齢ageとする重回帰分析を加えて実行してみましょう。
regs <-
list(
"model1" = lm(inc ~ educ, data = saving),
"model2" = lm(inc ~ educ + age, data = saving)
)
msummary(regs, gof_omit = 'AIC|BIC|Log.Lik.|F')| model1 | model2 | |
|---|---|---|
| (Intercept) | 1342.745 | −10858.405 |
| (1763.546) | (3250.631) | |
| educ | 742.530 | 869.852 |
| (146.062) | (137.566) | |
| age | 276.677 | |
| (63.872) | ||
| Num.Obs. | 100 | 100 |
| R2 | 0.209 | 0.337 |
| R2 Adj. | 0.201 | 0.323 |
２つの回帰分析の結果が並べて表示されました。
12.8 回帰分析の表をまとめる（stargazer）
【注意】以前のバージョンではstargazerパッケージを用いた方法を紹介していましたが、パッケージの更新が止まっているようです。前節で紹介したmodlsummaryパッケージを用いる方法を推奨しますが、stargazerの記述も記録のために残しておきます。（将来的に削除する可能性もあります。
ここでは、stargazerパッケージを用いて、複数の回帰分析の結果をまとめた表を作成してみる方法を説明します。
以下の2つの回帰分析をまとめて表にしてみましょう。
- 被説明変数を収入
inc、説明変数を教育年数educとする単回帰分析 - 被説明変数を収入
inc、説明変数を教育年数educと年齢ageとする重回帰分析
stargazer()では、回帰分析結果を保存したオブジェクトを引数とします。
regression1 <- lm(inc ~ educ, data = saving) #単回帰分析の結果をオブジェクトに保存
regression2 <- lm(inc ~ educ + age, data = saving) #重回帰分析の結果をオブジェクトに保存
stargazer(regression1, regression2, type = "html", out = "test.doc")| Dependent variable: | ||
| inc | ||
| (1) | (2) | |
| educ | 742.530*** | 869.852*** |
| (146.062) | (137.566) | |
| age | 276.677*** | |
| (63.872) | ||
| Constant | 1,342.745 | -10,858.410*** |
| (1,763.546) | (3,250.631) | |
| Observations | 100 | 100 |
| R2 | 0.209 | 0.337 |
| Adjusted R2 | 0.201 | 0.323 |
| Residual Std. Error | 4,992.593 (df = 98) | 4,593.594 (df = 97) |
| F Statistic | 25.844*** (df = 1; 98) | 24.646*** (df = 2; 97) |
| Note: | p<0.1; p<0.05; p<0.01 | |
1・2行目では、lm()で回帰分析を行い、結果をオブジェクトとしてregression1,regression2にそれぞれ入力しています。
stargazer()では、表に含める回帰分析の結果を引数とし、さらにオプションを指定します。
typeはデフォルトではLaTeXになっており多くの方は使えないと思いますので、ここではHTMLを指定して表形式に整えます。
outを指定すると、結果をファイルとして保存できます。
このファイルの表を卒業論文に貼り付け、必要に応じて加工すると良いでしょう。
表は列ごとに回帰分析の結果を示しています。 各行は説明変数を表し、かっこの無い上側が係数、かっこのある下側が標準誤差を表しています。 また、アスタリスク*で有意水準を表しています。 ***が有意水準1%で有意、**が有意水準5%で有意、*が有意水準10%で有意であることを示します。
線の下は以下のものを表しています。重要な情報なので、そのまま残しておきましょう。
Observation: 観測数- \(R^2\): 決定係数
Adjusted\(R^2\): 修正済み決定係数Residual Std. Error: 残差の標準誤差F Statistic: F値（F検定の結果も合わせて表示）