10 重回帰分析

Author
Affiliation

佐久間 智広

神戸大学大学院経営学研究科

Published

2024/06/17

Modified

2024/07/17

0.1 準備

Code 
if (!require("pacman")) install.packages("pacman")
pacman::p_load(tidyverse, magrittr, estimatr,
               car, modelsummary, ggrepel, patchwork,
               psych, ppcor, qgraph, AER, datarium, conflicted)
conflict_prefer("select", "dplyr")
conflict_prefer("filter", "dplyr")
1
複数のパッケージを一度に読み込めるpacmanパッケージが入っていない場合は，インストールする
2
パッケージの読み込み

主要なパッケージ

pacman

複数のパッケージを読み込めるパッケージ。そのパッケージの中のp_loadコマンドを使うと，（）内に指示したパッケージを読み込める。さらにそのパッケージがそもそもインストールされていない場合にはインストールした後に読み込んでくれる

tidyverse

Rのコードを楽にするたくさんのパッケージをパッケージにしたもの

magrittr

tidyverseに含まれる%>%などのコードをさらに拡張するもの。これがあると%$%とか%<>%が使える

1 前回の課題

kadai_8.csvには，ある大学の授業における学生19名の入学試験結果 (exam,700点満点)，入学から現在までのGPA (gpa)，高校の時の評定平均 (hschool, 10段階)が含まれています。ここから，GPAが入試成績とどのような関係にあるかを検証したいと考えています。

  1. describe関数を使って，gpaとexamの平均・標準偏差・中央値・最大値・最小値等を確認してください1。（コードのみでいいです）
  2. gpaとexamの相関係数，およびその検定をしてください。(コードと結果)
  3. gpaを従属変数，examを独立変数として回帰分析をしてください。その上で，係数とその統計的仮説検定の結果から，何が言えるかを説明してください（コードと説明のみ，結果は書く必要ないです）
  4. 回帰分析の結果から，predict関数を使って入試の点数がそれぞれ400,500,600,700点の時のgpaの予測値を計算したいです。まず，3で行った回帰分析の結果をkekka <- で保存しておいてください。次に以下のコードを実行して，まず入試の点数が400,500,600,700の時のデータを作り，newdataという名前で保存します。その後，predict(kekka, new = newdata)を使ってそれぞれの予測値を算出してください。（コードと結果）
newdata <- tibble(exam = c(400, 500, 600, 700))
  1. describe関数を使って，gpaとexamの平均・標準偏差・中央値・最大値・最小値等を確認してください2。（コードのみでいいです）
Code 
kadai <- read_csv("data/9_kadai.csv")
Code 
kadai |>  
  describe(skew = FALSE) |>  
  print(digits = 2)
        vars  n   mean    sd median    min    max  range    se
student    1 19  10.00  5.63  10.00   1.00  19.00  18.00  1.29
exam       2 19 534.05 71.91 518.00 440.00 690.00 250.00 16.50
gpa        3 19   2.51  0.76   2.63   1.14   3.87   2.73  0.17
hschool    4 19   7.04  0.92   7.30   5.50   8.80   3.30  0.21
  1. gpaとexamの相関係数，およびその検定をしてください。(コードと結果)
Code 
kadai %$% 
  cor.test(gpa,exam)

    Pearson's product-moment correlation

data:  gpa and exam
t = 4.8486, df = 17, p-value = 0.0001506
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.4703298 0.9034122
sample estimates:
      cor 
0.7617977
  1. gpaを従属変数，examを独立変数として回帰分析をしてください。その上で，係数とその統計的仮説検定の結果から，何が言えるかを説明してください（コードと説明のみ，結果は書く必要ないです）
Code 
lm(gpa ~ exam, data = kadai) -> kekka
msummary(kekka,
         gof_omit = "Log.Lik.|AIC|BIC|RMSE",
         title = "",          # タイトル
         stars = TRUE)
(1)
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
(Intercept) -1.797+
(0.896)
exam 0.008***
(0.002)
Num.Obs. 19
R2 0.580
R2 Adj. 0.556
F 23.509

examの係数は0.008で，0.1%水準で統計的に有意。これは，入試の点数が1点高いと，gpaが0.008高くなると期待されることを意味する。例えば入試の点数が100点高いと，gpaが0.8高くなることが期待される。

  1. 回帰分析の結果から，predict関数を使って入試の点数がそれぞれ400, 500, 600, 700点の時のgpaの予測値を計算したいです。まず，3で行った回帰分析の結果をkekka <- で保存しておいてください。次に以下のコードを実行して，まず入試の点数が400,500,600,700の時のデータを作り，newdataという名前で保存します。その後，predict(kekka,new = newdata)を使ってそれぞれの予測値を算出してください。（コードと結果）
newdata <- tibble(exam = c(400, 500, 600, 700))


Code 
newdata <- tibble(exam = c(400, 500, 600, 700))
predict(kekka, new = newdata)
       1        2        3        4 
1.429565 2.236329 3.043092 3.849855

2 重回帰分析とは

2.1 前回やったこと：回帰分析

  • 回帰分析では，2つの変数間の関係を線形関係（\(y = \beta_0 + \beta_1x\)という形）にモデル化して推定しました。

  • しかし，ある現象（例えば前回の例だったアイスクリーム屋さんの客数）が，1つの要因（気温）によって決まるということは稀です。

    • 曜日，天気，価格…

  • どれが最も重要なのかを知るために，例えば

    • \(\text{客数} = \alpha_0 + \alpha_1 \text{気温} + u_1\)
    • \(\text{客数} = \gamma_0 + \gamma_1 \text{価格} + u_2\)

  • みたいな形で別々に客数との関係を推定するのはあまり効率的ではありません。そこで使われるのが重回帰分析です。

2.2 概要

  • 重回帰分析では，以下のように複数の独立変数を一つの式に含めます。

\[ y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + u \]

  • これは図のように2つの独立変数が従属変数に与える影響を表しています。

  • このように複数の要因を一度に分析するものを重回帰分析と言います（独立変数が一つの回帰分析を単回帰分析と言います）

2.3 重回帰分析の良いところ

他の要因を一定とした場合の関係を推定できる

  • 重回帰分析は，単回帰分析を複数行うこととは違い，「他の要因をコントロールした」係数を推定できる点が特徴です。
    • その他の要因を一定としたとき，ある独立変数が従属変数とどのように関係しているかがわかる
      • この「他の条件がすべて等しければ」はceteris paribusというラテン語で表記されたりします
      • 上記の例の場合，価格の係数は「気温を一定としたとき，価格は客数とどのように関係しているか」が分かる


欠落変数バイアスの除去

  • 独立変数，従属変数両方に関係する要因を含めずに回帰分析をすることで欠落変数バイアスが生じる
    • 回帰分析の仮定が満たされず，推定結果がゆがむ
  • その要因を独立変数として回帰モデルに含めると，バイアスが除去される

2.4 重回帰分析の意義他の要因を一定とした場合の関係を推定できる

2.4.1 例: ラーメンの評価

  • ラーメン好き500人に聞いた，ラーメンの全体評価と個別要素の評価
    • スープ
    • 麺
    • チャーシュー
    • その他トッピング
  • 仮説：ラーメンの評価を決めるもの…スープ！
  • 下記のような単回帰分析を推定

\[ ramen_i = \beta_0 + \beta_1 soup_i +\epsilon_i \]

  • ramen: ラーメンの総合評価
  • soup: スープの評価

moccos.jp

moccos.jp

\(\beta_1\)が正で，統計的に有意だった。ここからスープ評価が大事だとわかる。

ここで

スープだけでなく麺の食感も大事でしょう

チャーシューの味や分厚さが大事

味玉の美味しさでラーメン屋の仕事の丁寧さがわかる

とか言ってくる人も？結局どれがどれぐらい大事？


じゃあ関係ありそうな要因を全て入れてみて

\[ ramen_i = \beta_0 + \beta_1 soup_i + \beta_2 men_i + \beta_3 niku_i + \beta_4 negi_i + \epsilon_i \tag{1}\]

とできる3。こうすることで

  • soupの係数\(\beta_1\)は，「麺，チャーシュー，トッピングの評価を一定とした時のスープの評価とラーメン総合評価との関係」を表します。
  • ここで，関係を検証したい変数以外は，条件を一定にする，つまりコントロールするために使われているのでコントロール変数と呼びます。

つまり，上記のような疑問に（ある意味）対処した上でのスープの評価が大事かどうかがわかる。

2.5 欠落変数バイアスの問題

  • 重回帰分析のもう1つの重要な役割が，欠落変数バイアスに対処できる点にあります。
    • コントロールするべき変数をコントロールしないと，分析結果が不適切になります。
    • 不偏性，一致性が失われ，歪んだ推定結果になる
    • これを内生性の問題と呼びます。詳しくは別の回で説明します

3 分析方法

3.1 データ

今回使うのは，以前課題で使ったアイスクリーム屋さんの例です (アイスクリーム統計学より)

データの読み込み

Code 
ice10 <- read_csv("data/10_ice.csv")
ice10

データには，その日の最高気温saikou，最低気温saitei，客数kyakuが含まれています。

3.2 単回帰分析

Code 
reg1 <- lm(kyaku ~ saikou, data = ice10)
reg2 <- lm(kyaku ~ saitei, data = ice10) 

regs <-list(reg1, reg2)
msummary(regs,
         stars = TRUE,
         gof_omit = "Log.Lik.|AIC|BIC|RMSE|F",
         )
(1) (2)
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
(Intercept) -229.982** 118.600
(73.787) (151.256)
saikou 17.248***
(2.300)
saitei 8.100
(6.029)
Num.Obs. 20 20
R2 0.758 0.091
R2 Adj. 0.744 0.041

これは，最高気温が客数と有意ににプラスの関係にある一方で，最低気温と客数に関係が見られない，ということを示唆しています。

関連して，普通の相関係数を求めてみる。データのうち，通し番号(Num)は必要ないので，dplyr::select()コマンドで，使う変数(saikou, saitei, kyaku)だけを選択し，そのあとcor()関数を使う。

Code 
cor5 <- ice10 |>  
  select(saikou, saitei, kyaku) |>
  cor() 

cor5
          saikou    saitei     kyaku
saikou 1.0000000 0.7058315 0.8703519
saitei 0.7058315 1.0000000 0.3019116
kyaku  0.8703519 0.3019116 1.0000000


最高気温と客数は強い正の相関（0.87），最低気温と客数はあまり強くない相関（0.30）最高気温と最低気温も強い正の相関（0.71）

3.3 重回帰分析

重回帰分析は，単回帰分析の独立変数部分を複数にして，「+」で繋ぐ。

Code 
reg3 <- lm(kyaku ~ saitei + saikou, data = ice10) 
summary(reg3)

Call:
lm(formula = kyaku ~ saitei + saikou, data = ice10)

Residuals:
    Min      1Q  Median      3Q     Max 
-17.274  -9.342  -2.920  10.810  26.392 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -90.640     37.718  -2.403   0.0279 *  
saitei       -16.703      2.012  -8.301 2.20e-07 ***
saikou        25.957      1.486  17.463 2.71e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 13.52 on 17 degrees of freedom
Multiple R-squared:  0.952, Adjusted R-squared:  0.9464 
F-statistic: 168.6 on 2 and 17 DF,  p-value: 6.161e-12

重回帰分析を先ほどの2本の単回帰分析と並べてみる

Code 
regs <-list(reg1, reg2, reg3) 
msummary(regs,
         stars = TRUE,
         gof_omit = "Log.Lik.|AIC|BIC|RMSE|F",
         )
(1) (2) (3)
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
(Intercept) -229.982** 118.600 -90.640*
(73.787) (151.256) (37.718)
saikou 17.248*** 25.957***
(2.300) (1.486)
saitei 8.100 -16.703***
(6.029) (2.012)
Num.Obs. 20 20 20
R2 0.758 0.091 0.952
R2 Adj. 0.744 0.041 0.946
  • 最低気温を調整したら，最高気温と客数の関係はより強く（係数17 → 26, p < 0.001）
  • 最低気温は，有意に負の関係に（係数-17, p < 0.001）

3.4 決定係数と自由度調整済み決定係数

  • 分析結果の当てはまり度合いを表す指標（のひとつ）が決定係数 (\(R^2\))
    • データのばらつきのうち，推定された回帰式で説明できる割合を表す。
      • 回帰分析がどれぐらいデータを予測できているか？
    • 0から1まで，1に近いほど当てはまりが良い
  • 単回帰分析ではR-squaredとAdjusted R-squaredという2つの種類は大差ないという話だった
  • 決定係数には「独立変数を増やせば増やすほど改善する」という特徴がある。
    • いい変数でもそうでなくても改善する
    • 変数入れれば入れるほど良いのか？
  • この，変数入れれば入れるほど改善してしまうというのを調整したのが自由度調整済み決定係数Adjusted R-squared
    • 重回帰分析ではこっちを見る

3.5 多重共線性

  • 説明変数の間で強い相関があるものをモデルに取り入れることで引き起こされる問題
  • 3つまたはそれ以上の変数が共線性（予測変数間に強い相関）をもつと多重共線性を引き起こす
  • 多重共線性が起こると推定結果がぐちゃぐちゃになることも。

回帰から問題のある変数のうち1つを取り除く分散拡大要因 (Variance Inflation Factor: VIF)の値が5または10を上回る変数がないか確認する

Code 
lm(kyaku ~ saitei + saikou, data = ice10) |>  
  vif()
  saitei   saikou 
1.992818 1.992818

今回は問題ない

  • 問題があった場合，数値が高い変数のどちらか一方を取り除く

4 実習課題

前回と同じデータで分析を続けます。9_kadai.csvには，ある大学の授業における学生19名の入学試験結果 (exam,700点満点)，入学から現在までのGPA (gpa)，高校の時の評定平均 (hschool, 10段階)が含まれています。ここから，GPAが入学時点で得られる情報とどのような関係にあるかを検証したいと考えています。

  1. gpaとexam,hschool3つの変数の記述統計量をdescribe関数を使って算出してください。（コード）
  2. gpaとexam,hschool3つの変数の相関係数を計算してください。cor()関数で，3つ以上の変数が含まれるデータを()内に指定すると，相関表が出力されます。（コード）
  3. gpaを従属変数，examを独立変数として単回帰分析をしてください（前回課題と同じ）
  4. gpaを従属変数，hschoolを独立変数として単回帰分析をしてください。その上で，係数とその統計的仮説検定の結果から，何が言えるかを説明してください（コードと説明のみ，結果は書く必要ないです）
  5. gpaを従属変数，examとhschoolを独立変数として，重回帰分析をしてください。その上で，係数とその統計的仮説検定の結果から，何が言えるかを説明してください（コードと説明のみ，結果は書く必要ないです）

/

Footnotes

  1. describe()関数は，psychパッケージの関数です。関数が見つからないというエラーが出た場合は，パッケージのインストール，読み込みがされているか確認してください。この授業の最初の 準備 のところで読み込んでいます。↩︎

  2. describe()関数は，psychパッケージの関数です。関数が見つからないというエラーが出た場合は，パッケージのインストール，読み込みがされているか確認してください。この授業の最初の 準備 のところで読み込んでいます。↩︎

  3. men, niku, negiはそれぞれ麺，チャーシュー，トッピングの評価↩︎