1 Simple regression analysis / 単回帰分析

1つの説明変数をもつ次のような回帰モデル（単回帰モデル）を考える．

\[ Fertility = \alpha + \beta Exam + u \]

lm (linear model) 関数で推定する．

# str(swiss)  # built-in dataset
lm(formula = Fertility ~ Examination, data = swiss)

## 
## Call:
## lm(formula = Fertility ~ Examination, data = swiss)
## 
## Coefficients:
## (Intercept)  Examination  
##      86.819       -1.011

lm 関数で計算した回帰オブジェクトを summary 関数に渡すと決定係数なども表示される．

summary(lm(formula = Fertility ~ Examination, data = swiss))

## 
## Call:
## lm(formula = Fertility ~ Examination, data = swiss)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -25.9375  -6.0044  -0.3393   7.9239  19.7399 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  86.8185     3.2576  26.651  < 2e-16 ***
## Examination  -1.0113     0.1782  -5.675 9.45e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.642 on 45 degrees of freedom
## Multiple R-squared:  0.4172, Adjusted R-squared:  0.4042 
## F-statistic: 32.21 on 1 and 45 DF,  p-value: 9.45e-07

1.1 Estimation by hand

回帰係数の最小二乗（OLS）推定量は次のように与えられる．

\[ \hat{\alpha} = \bar{Y} - \hat{\beta} \bar{X}, \quad \hat{\beta} = \frac{S_{XY}}{S_{XX}} \] \[ S_{XX} = \sum_i (X_i - \bar{X})^2, \quad S_{XY} = \sum_i (X_i - \bar{X}) (Y_i - \bar{Y}) \]

まず $\bar{Y}, \bar{X}, S_{XX}, S_{XY}$ という4つの統計量を計算する．

y_bar <- mean(swiss$Fertility)  # Y bar
x_bar <- mean(swiss$Examination)  # X bar
s_xx <- sum((swiss$Examination - x_bar)^2)  # S_XX
s_xy <- sum((swiss$Examination - x_bar) * (swiss$Fertility - y_bar))  # S_XY
c(y_bar, x_bar, s_xx, s_xy)

## [1]    70.14255    16.48936  2927.74468 -2960.87872

計算した統計量を使って回帰係数を計算する．

beta_hat <- s_xy / s_xx
alpha_hat <- y_bar - beta_hat * x_bar
c(alpha_hat, beta_hat)

## [1] 86.818529 -1.011317

1.2 Estimation by numerical calculation / 数値計算による推定

OLSは残差二乗和の最小化によって回帰係数パラメタを推定する方法（モデルをデータに当てはめる原理）であるが，これを解析的に解くのではなく数値計算によって解いてみる．

残差二乗和の関数（目的関数） $Q$ の最小化は R に built-in されている汎用の最適化関数 optim によって行う．

1.2.1 Example of `optim`

例：次の二次関数は $x=2$ で最小値 3 をとる．

\[ f(x) = x^2 -4x + 7 = (x-2)^2 + 3 \]

最小値を取るような $x$ の値を optim 関数を使って求めてみる．

optim(par = パラメタの初期値ベクトル, fn = 目的関数) の形で指定
- $par ：最適化されたパラメタ（＝目的関数の引数；この例では $x$）
- $value ：最小化された目的関数の値
- $convergence ：最適化計算が収束していれば 0

obj_function <- function (x) x^2 - 4*x + 7
optim(par = 0, fn = obj_function, method = "BFGS")  # ニュートン法の亜種を利用

## $par
## [1] 2
## 
## $value
## [1] 3
## 
## $counts
## function gradient 
##        6        3 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

1.2.2 In the case of regression

回帰分析の場合，目的関数は残差二乗和 $Q$ であり，その引数は推定したい回帰係数パラメタ $(\alpha, \beta)$ となる．

\[ (\hat{\alpha}, \hat{\beta}) = \arg \min_{\alpha, \beta} Q(\alpha, \beta) \] \[ Q (\alpha, \beta) := \sum_i \hat{u}_i^2 = \sum_i (Y_i - \alpha - \beta X_i)^2 \]

sum_of_squared_residuals <- function (coef_vector) {
  # coef_vector[1]: alpha (intercept)
  # coef_vector[2]: beta (slope of X)
  y <- swiss$Fertility
  x <- swiss$Examination
  sum((y - coef_vector[1] - coef_vector[2] * x)^2)
}
sum_of_squared_residuals(coef_vector = c(1, 2))  # alpha=1, beta=2

## [1] 92200.11

sum_of_squared_residuals(coef_vector = c(1, 3))  # alpha=1, beta=3

## [1] 69485.91

optim(par = c(0, 0), fn = sum_of_squared_residuals)

## $par
## [1] 86.822437 -1.011885
## 
## $value
## [1] 4183.569
## 
## $counts
## function gradient 
##       93       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

$\hat{\alpha} = 86.8$ に固定して $\hat{\beta}$ の値によって残差二乗和がどのように変化するかをプロットすると，$\hat{\beta} = -1.01$ で最小化されていることが分かる．

（複雑なスクリプトなので，どのように動いているかは理解しなくてよい．）

sum_of_squared_residuals_beta <- function (beta1) {
  y <- swiss$Fertility
  x <- swiss$Examination
  sum((y - 86.8 - beta1 * x)^2)
}
beta_value <- seq(-3, 2, by = 0.01)
obj_fun <- sapply(beta_value, sum_of_squared_residuals_beta)
plot(x = beta_value, y = obj_fun, type = "l", xlim = c(-3, 2))
abline(v = -1.01, col = 2)

1.3 Coefficient of determination / 決定係数

\[ R^2 = \frac{\mbox{Variation in } Y \mbox{ explained by } X}{\mbox{Total variation in } Y} = \frac{S_{\hat{Y}\hat{Y}}}{S_{YY}} = 1 - \frac{\sum \hat{u}_i^2}{S_{YY}} \]

\[ S_{YY} = \sum_i (Y_i - \bar{Y})^2, \quad S_{\hat{Y} \hat{Y}} = \sum_i (\hat{Y}_i - \bar{\hat{Y}})^2 \]

s_yy <- sum((swiss$Fertility - y_bar)^2)
y_pred <- alpha_hat + beta_hat * swiss$Examination  # y hat
s_yhyh <- sum((y_pred - y_bar)^2)
s_yhyh / s_yy  # R^2

## [1] 0.4171645

resid_swiss <- swiss$Fertility - y_pred  # residuals
1 - sum(resid_swiss^2) / s_yy  # R^2 ... should be the same as above

## [1] 0.4171645

決定係数が 0.417 ということは，説明変数 (Examination) によって被説明変数 (Fertility) の変動のうち 41.7% を説明することができる，という意味．

余裕のある履修者向けの宿題：自由度調整済み決定係数を計算．

\[ R^2 = 1 - \frac{n - 1}{n - K - 1} \frac{\sum \hat{u}_i^2}{S_{YY}} \]

$K$: 説明変数の数（定数項は含めない）

1.4 Standard error and test / 標準誤差と検定

帰無仮説「$H_0: \beta = 0$」の検定を行うために，$\hat{\beta}$ の標準誤差（s.e.）を求める．

\[ V(\hat{\beta}) = \frac{\sigma^2}{S_{XX}}, \quad s.e.(\hat{\beta}) = \sqrt{V(\hat{\beta})} \]

誤差分散 $\sigma^2$ は残差 $\hat{u}$ を利用して次のように推定する．

\[ \hat{\sigma}^2 = s^2 = \frac{1}{n - 2} \sum \hat{u}_i^2 \]

sum(resid_swiss^2) / (nrow(swiss) - 2)  # estimated sigma^2

## [1] 92.96816

beta_se <- sqrt(sum(resid_swiss^2)/(nrow(swiss) - 2) / s_xx)  # standard error
beta_se

## [1] 0.1781971

\[ t = \frac{\hat{\beta}}{\sqrt{V(\hat{\beta})}} \]

beta_hat / beta_se

## [1] -5.675275

$t < 0$ なので次のように p 値を計算できる．

\[ p = 2 \times \Pr(T \le t) \]

2 * pt(q = beta_hat / beta_se, df = nrow(swiss) - 2)

## [1] 9.450437e-07

補足：$t>0$ の場合も考慮してより一般的に書けば $p = 2 \times [1-\Pr(|t| \le T)]$．

2 * (1 - pt(q = abs(beta_hat / beta_se), df = nrow(swiss) - 2))

## [1] 9.450437e-07

1.4.1 cf. Critical value approach / 臨界値を用いて検定を行う場合

t分布の2.5%臨界値（両側5%に対応）は？

qt(p = 0.025, df = nrow(swiss) - 2)  # df: degree of freedom (自由度)

## [1] -2.014103

qt(p = c(0.005, 0.025, 0.975, 0.995), df = nrow(swiss) - 2)  # 1% and 5%

## [1] -2.689585 -2.014103  2.014103  2.689585

よって，$t = -5.7$ は有意水準1%に対応する臨界値 (-2.7) よりも小さいので，1%水準で有意である．

自由度が大きければ標準正規分布の場合とほとんど同じ臨界値になる．

qt(p = 0.975, df = 10000)

## [1] 1.960201

2 Multiple regression analysis / 重回帰分析

2つの説明変数をもつ次のような回帰モデル（重回帰モデル）を考える．

\[ Fertility = \alpha + \beta_1 Exam + \beta_2 Educ + u \]

summary(lm(Fertility ~ Examination + Education, data = swiss))

## 
## Call:
## lm(formula = Fertility ~ Examination + Education, data = swiss)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.9935  -6.8894  -0.3621   7.1640  19.2634 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  85.2533     3.0855  27.630   <2e-16 ***
## Examination  -0.5572     0.2319  -2.402   0.0206 *  
## Education    -0.5395     0.1924  -2.803   0.0075 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.982 on 44 degrees of freedom
## Multiple R-squared:  0.5055, Adjusted R-squared:  0.483 
## F-statistic: 22.49 on 2 and 44 DF,  p-value: 1.87e-07

2.1 Estimation by hand

重回帰モデルの回帰係数のOLS推定量は行列形式で表現した方が分かりやすい．

式中の $\mathbf{X}, \mathbf{Y}$ は（$X, Y$ とは異なり）太字かつ立体（斜体ではない）で，これは行列またはベクトルを表す．

つまり，以下のように定義されている． $X_{nm}$ は $n$ 番目の個体の $m$ 番目の説明変数を表す． $Y_{n}$ は $n$ 番目の個体の被説明変数（アウトカム）を表す．

\[ \mathbf{X} = \left[ \matrix{X_{11} & X_{12} & X_{13} & \cdots \\ X_{21} & X_{22} & X_{23} & \cdots \\ X_{31} & X_{32} & X_{33} & \cdots \\ \vdots & \vdots & \vdots & \ddots} \right], \quad \mathbf{Y} = \left[ \matrix{Y_1 \\ Y_2 \\ Y_3 \\ \vdots} \right] , \quad \boldsymbol{\beta} = \left[ \matrix{\beta_1 \\ \beta_2 \\ \beta_3 \\ \vdots} \right] \]

行列表記の回帰モデル．

\[ \mathbf{Y} = \mathbf{X} \boldsymbol{\beta} + \mathbf{u} \]

残差二乗和 $Q$ を $\boldsymbol{\beta}$ で微分して 0 とおくことでOLS推定量が得られる．

\[ Q = \mathbf{u}'\mathbf{u} = (\mathbf{Y} - \mathbf{X} \boldsymbol{\beta})' (\mathbf{Y} - \mathbf{X} \boldsymbol{\beta}), \quad \frac{\partial Q}{\partial \boldsymbol{\beta}} =0 \]

\[ \hat{\boldsymbol{\beta}} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y} \]

まずは $\mathbf{X}, \mathbf{Y}$ を作成する．説明変数には1（定数項）が含まれている点に注意しよう．

as.matrix(データフレーム) ：行列オブジェクトに変換

Y <- as.matrix(swiss[, c("Fertility")])
X <- as.matrix(cbind(1, swiss[, c("Examination", "Education")]))

転置は t，逆行列は solve，行列の積は %*% で計算できるので，回帰係数のOLS推定量は以下のように計算される．

beta_multi <- solve(t(X) %*% X) %*% t(X) %*% Y
beta_multi

##                   [,1]
## 1           85.2532753
## Examination -0.5572183
## Education   -0.5394570

2.2 Standard error / 標準誤差

\[ V(\hat{\boldsymbol{\beta}}) = \sigma^2 (\mathbf{X}'\mathbf{X})^{-1} \]

回帰係数の検定を行うためには誤差分散 $\sigma^2$ の推定値を求める必要がある． $k$ は（定数項を除いた）説明変数の数なので2である．

\[ \hat{\sigma}^2 = \frac{1}{n-(1+k)} \sum_i (Y_i - \hat{Y}_i)^2 \]

resid_multi <- Y - X %*% beta_multi
s2 <- sum(resid_multi^2) / (nrow(swiss)-3)  # estimates of sigma^2
s2

## [1] 80.67295

回帰係数の推定量ベクトルの分散は $\sigma^2 (\mathbf{X}'\mathbf{X})^{-1}$ （分散共分散行列）なので，個々の回帰係数の推定量の分散はこの分散共分散行列の対角成分である．

s2 * solve(t(X) %*% X)  # variance-covariance matrix of beta

##                      1 Examination   Education
## 1            9.5202985 -0.54480732  0.10745077
## Examination -0.5448073  0.05379495 -0.03117276
## Education    0.1074508 -0.03117276  0.03703237

diag(s2 * solve(t(X) %*% X))  # extract diagonals

##           1 Examination   Education 
##  9.52029846  0.05379495  0.03703237

標準誤差はその平方根となる．すなわち，$j$ 番目の説明変数の標準誤差は次のように計算される．

\[ s.e.(\hat{\beta}_j) = \sqrt{\{\sigma^2 (\mathbf{X}'\mathbf{X})^{-1} \}_{jj}} \]

sqrt(diag(s2 * solve(t(X) %*% X)))  # standard errors

##           1 Examination   Education 
##   3.0854981   0.2319374   0.1924380

回帰係数の推定値を標準誤差で割れば t 値が求められる．

t_value <- beta_multi / sqrt(diag(s2 * solve(t(X) %*% X)))
t_value

##                  [,1]
## 1           27.630312
## Examination -2.402451
## Education   -2.803277

t 値が分かれば p 値も計算できる．

\[ p = 2 \times (1 - \Pr(|t| \le T)) \]

2 * (1 - pt(q = abs(t_value), df = nrow(swiss) - 3))

##                    [,1]
## 1           0.000000000
## Examination 0.020571604
## Education   0.007497224

3 Advanced method / 応用

3.1 Adding quardatic or interaction terms / 二次・交差項

新しい変数を作成するなどのデータフレームの操作をするためには tidyverse に含まれる関数 (mutate, rename, group_by, etc.) を用いると簡単．

library(tidyverse)
swiss2 <- swiss %>%
  mutate(Edu2 = Education^2,
         ExamEdu = Examination * Education)

# tidyverse を使わない場合
# swiss2 <- swiss
# swiss2$Edu2 <- swiss$Education^2

\[ Fertility = \alpha + \beta_1 Exam + \beta_2 Educ^2 + u \]

lm(Fertility ~ Examination + Edu2, data = swiss2)

## 
## Call:
## lm(formula = Fertility ~ Examination + Edu2, data = swiss2)
## 
## Coefficients:
## (Intercept)  Examination         Edu2  
##    83.22072     -0.66069     -0.01035

\[ Fertility = \alpha + \beta_1 Exam \times Educ + u \]

lm(Fertility ~ ExamEdu, data = swiss2)

## 
## Call:
## lm(formula = Fertility ~ ExamEdu, data = swiss2)
## 
## Coefficients:
## (Intercept)      ExamEdu  
##    75.73211     -0.02394

データセットを直接操作せずに，回帰モデルを指定する際に I 関数を使って二乗項を表現することもできる．

lm(Fertility ~ Examination + I(Education^2), data = swiss)

## 
## Call:
## lm(formula = Fertility ~ Examination + I(Education^2), data = swiss)
## 
## Coefficients:
##    (Intercept)     Examination  I(Education^2)  
##       83.22072        -0.66069        -0.01035

「var1 : var2」は「var1 × var2」を表し，「var1 * var2」は「var1 + var2 + (var1 × var2)」を表す．

lm(Fertility ~ Examination : Education, data = swiss)

## 
## Call:
## lm(formula = Fertility ~ Examination:Education, data = swiss)
## 
## Coefficients:
##           (Intercept)  Examination:Education  
##              75.73211               -0.02394

\[ Fertility = \alpha + \beta_1 Exam + \beta_2 Educ + \beta_3 Exam \times Educ + u \]

lm(Fertility ~ Examination * Education, data = swiss)

## 
## Call:
## lm(formula = Fertility ~ Examination * Education, data = swiss)
## 
## Coefficients:
##           (Intercept)            Examination              Education  
##             87.178104              -0.625731              -0.807552  
## Examination:Education  
##              0.009201

説明変数や被説明変数に対数を取る場合は I なしで log を使ってよい．

\[ \log(Fertility) = \alpha + \beta \ Exam + u \]

lm(log(Fertility) ~ Examination, data = swiss)

## 
## Call:
## lm(formula = log(Fertility) ~ Examination, data = swiss)
## 
## Coefficients:
## (Intercept)  Examination  
##     4.49251     -0.01574

3.2 Dummy variables ダミー変数

Examination 変数が平均よりも大きい場合に 1 をとるダミー変数を定義してみる．

swiss2 <- swiss %>%
  mutate(exam_dummy = ifelse(Examination > mean(Examination), 1, 0))

\[ Fertility = \alpha + \beta \ 1(Exam > \bar{Exam}) + u \]

summary(lm(Fertility ~ exam_dummy, data = swiss2))

## 
## Call:
## lm(formula = Fertility ~ exam_dummy, data = swiss2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -31.359  -5.557   1.945   6.793  16.441 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   76.059      2.020  37.657  < 2e-16 ***
## exam_dummy   -13.904      3.096  -4.491 4.91e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.5 on 45 degrees of freedom
## Multiple R-squared:  0.3095, Adjusted R-squared:  0.2941 
## F-statistic: 20.17 on 1 and 45 DF,  p-value: 4.907e-05

以下のように回帰式を指定することでダミー変数を表現することもできる． Examination > mean(Examination) という変数は FALSE と TRUE のいずれかをとり，そのうちの一方（TRUE）に該当する場合に 1 をとるダミー変数として推定が行われる．すなわち，この変数の回帰係数は，もう一方（FALSE）の場合をベースラインとして，TRUE の場合に切片が相対的にどれだけ異なるかを表す．

lm(Fertility ~ I(Examination > mean(Examination)), data = swiss2)

（出力は省略）

3.2.1 Relationship to the test of difference of means / 平均値の差の検定との関係

exam_dummy の回帰係数 (-13.9) は，「Examination が平均を上回るサブサンプルにおける Fertility の平均」と「Examination が平均以下のサブサンプルにおける Fertility の平均」の差に等しい．

swiss2 %>%
  group_by(exam_dummy) %>%
  summarise(mean_fertility = mean(Fertility))

## # A tibble: 2 × 2
##   exam_dummy mean_fertility
##        <dbl>          <dbl>
## 1          0           76.1
## 2          1           62.2

また，回帰係数の t 検定は2標本の母平均の差の検定（等分散を仮定）に対応する．

t.test(Fertility ~ exam_dummy, data = swiss2, var.equal = TRUE)

## 
##  Two Sample t-test
## 
## data:  Fertility by exam_dummy
## t = 4.4906, df = 45, p-value = 4.907e-05
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
##   7.667982 20.140536
## sample estimates:
## mean in group 0 mean in group 1 
##        76.05926        62.15500

このように，線形回帰モデルにおける回帰係数に対する検定は幾つかの検定と本質的に同じことを行っているものとして解釈できる．詳細は Lindeløv (2019) Common statistical tests are linear models を参照．

3.3 Categorical variables

複数の値を持つカテゴリー変数（質的変数）はそのまま回帰式に含めることで自動的にダミー変数として処理される．

いずれか一つのカテゴリーはベースラインとされる（すべてのカテゴリーをダミー変数として含めると多重共線性により推定できなくなるため，一つは除外される）．

swiss データで，以下の6地域をカテゴリー変数として利用する（ChatGPT が提示してくれたもので，正確性・妥当性は未検証）．

A: ジュラカントン, B: ヴォー州, C: フリブール州, D: ネウシャテル州, E: ジュネーブ州, F: ヴァレー州

swiss2 <- swiss
swiss2$region <- 
  c("A", "A", "A", "A", "A", "A", "C", "C", "C", "C",
    "C", "B", "B", "B", "B", "B", "B", "B", "B", "B",
    "B", "B", "B", "B", "B", "B", "B", "B", "B", "B",
    "F", "F", "F", "F", "F", "F", "F", "F", "D", "D",
    "D", "D", "D", "D", "E", "E", "E")
table(swiss2$region)

## 
##  A  B  C  D  E  F 
##  6 19  5  6  3  8

lm(Fertility ~ Examination + region, swiss2)$coef

## (Intercept) Examination     regionB     regionC     regionD     regionE 
##   90.339456   -0.741199  -12.142443    6.053810   -1.829081  -30.976148 
##     regionF 
##   -8.810613

regionB の回帰係数 -12.1 は，region = A をベースラインとしたときの B の切片の違い（平均的に Fertility がどれだけ異なるか）を表す．

3.3.1 Change the reference group

ベースラインのカテゴリーを変更するには relevel 関数を使う．

swiss2$region <- relevel(as.factor(swiss2$region), ref = "B")
lm(Fertility ~ Examination + region, swiss2)$coef

## (Intercept) Examination     regionA     regionC     regionD     regionE 
##   78.197013   -0.741199   12.142443   18.196253   10.313362  -18.833705 
##     regionF 
##    3.331830

3.4 Multicollinearity / 多重共線性

3.4.1 Example of perfect multico. / 完全な多重共線性の例

lm(formula = Fertility ~ Examination + I(2 * Examination), data = swiss)$coef

##        (Intercept)        Examination I(2 * Examination) 
##          86.818529          -1.011317                 NA

cor(swiss$Examination, 2 * swiss$Examination)

## [1] 1

R では，該当する変数（のうち一つ）が自動的に drop される．

3.4.2 Example of imperfect multico. / 不完全な多重共線性の例

説明変数 Examination と非常に相関の高い変数が説明変数として追加される場合．

set.seed(2)
swiss_multico <- swiss %>% 
  mutate(Exam_with_noise = Examination + rnorm(nrow(swiss), sd = 1))
cor(swiss_multico$Examination, swiss_multico$Exam_with_noise)

## [1] 0.9900738

summary(lm(formula = Fertility ~ Examination + Exam_with_noise, data = swiss_multico))$coef

##                   Estimate Std. Error    t value     Pr(>|t|)
## (Intercept)     86.6475143   3.202510 27.0561290 4.660541e-29
## Examination      0.9822164   1.245743  0.7884584 4.346563e-01
## Exam_with_noise -1.9706250   1.219204 -1.6163209 1.131734e-01

（参考）Examination のみの場合．

summary(lm(formula = Fertility ~ Examination, data = swiss_multico))$coef

##              Estimate Std. Error   t value     Pr(>|t|)
## (Intercept) 86.818529  3.2576034 26.651043 3.353924e-29
## Examination -1.011317  0.1781971 -5.675275 9.450437e-07

3.4.3 Using VIF for diagnostics / VIF による診断

変数 $j$ の分散拡大係数（variance inflation factor, VIF）は，変数 $j$ をそれ以外の変数に回帰したときの決定係数 $R_j^2$ を使って以下のように計算される．

\[ VIF_j = \frac{1}{1 - R_j^2} \]

あくまで一つの経験則でしかないが，「$VIF > 10$ ならば多重共線性が生じている」と判断される．

lm_exam <- lm(Examination ~ Education + Agriculture + Catholic, data = swiss)
r2_exam <- summary(lm_exam)$r.squared
1 / (1 - r2_exam)  # VIF

## [1] 3.675372

car パッケージの vif 関数を使うと便利．

# install.packages("car")
library(car)
lm0 <- lm(formula = Fertility ~ Examination + Education + Agriculture + Catholic, data = swiss)
car::vif(lm0)

## Examination   Education Agriculture    Catholic 
##    3.675372    2.689400    2.147410    1.856475

3.5 Comparing the size of the coefficients / 回帰係数の大きさを比較する

例：Examination と Agriculture の回帰係数の（絶対値の）大きさを比較する．

lm(formula = Fertility ~ Examination + Agriculture, data = swiss)$coef

## (Intercept) Examination Agriculture 
## 94.60968827 -1.19502875 -0.09399751

上の結果から「$|-1.195|>|-0.094|$ なので Examination の方が Fertility の分散をより説明する」と結論付けることはできない．

なぜなら，係数の大きさは変数のスケール（ばらつき）に依存するから．

たとえば，線形モデル $y = \beta_0 + \beta_1 x + u$ において金額ベースの変数 $x$ の単位を「円」から「千円」に変えると，推定される値そのもの ($\hat{\beta}_1$) は1千倍になるが，推定された回帰係数の「解釈」は全く変わらない（この変数の重要性が1千倍になるわけではない）．

3.5.1 Comparing the change per one SD / 1標準偏差当たりの変化量の比較

-1.195 * sd(swiss$Examination)

## [1] -9.533571

-0.094 * sd(swiss$Agriculture)

## [1] -2.134854

「$|-9.5|>|-2.1|$ なので Examination の方が Fertility の分散をより説明する」と言える．

3.5.2 Standardized partial regression coefficient / 標準偏回帰係数

より一般的な方法は，標準化（平均 0，分散 1 になるように変換）した変数を使って回帰分析をするもの．

swiss_scaled <- swiss %>% 
  mutate(Fertility_scaled = (Fertility - mean(Fertility)) / sd(Fertility),
         Examination_scaled = (Examination - mean(Examination)) / sd(Examination),
         Agriculture_scaled = (Agriculture - mean(Agriculture)) / sd(Agriculture))
lm(formula = Fertility_scaled ~ Examination_scaled + Agriculture_scaled, data = swiss_scaled)$coef

##        (Intercept) Examination_scaled Agriculture_scaled 
##       3.846942e-16      -7.632109e-01      -1.708973e-01

「$|-0.76|>|-0.17|$ なので Examination の方が Fertility の分散をより説明する」と言える．

スケーリング（分散を 1 にする）とセンタリング（平均を 0 にする）を同時に行ってくれる scale という関数を使えばより簡単．

lm(formula = Fertility ~ Examination + Agriculture, data = as.data.frame(scale(swiss)))

（出力は省略）

3.6 Model selection / モデル選択

被説明変数の予測を目的として回帰分析を行う場合，複数のモデル候補があった場合に，予測力が最も優れたモデルを一つ選択したい場合があるだろう．しかし，説明変数を追加すれば決定係数は必ず増加するが，それは標本に特有のデータの変動（あるいはノイズ）に過適合しているだけかもしれない．これを防ぐために，モデルの説明力とモデルの簡潔さのトレードオフの間でバランスを取るための何らかの指標が必要となる．

自由度調整済みの決定係数はモデル選択の趣旨には当てはまっているが，実務的により広く利用されるのはAIC（赤池情報量規準）と呼ばれる規準（基準）．これはモデルを最尤推定することを前提として「最尤推定量の対数尤度」と「パラメタの数（＝定数項を含む説明変数の数＝K+1）」のトレードオフとして記述される．

\[ AIC = -2 \log L + 2(K+1) \]

対数尤度は大きいほど当てはまりが良く，パラメタ数は少ない方がよい．すなわち，複数のモデル間で AIC が最も小さいモデルが望ましい．

注：

因果推論を目的として回帰分析を行う場合には一般にモデル選択は行わない場合が多い．分析の目的に応じてモデル選択を行うかどうかを決めよう．
- 典型的な AIC の使いどころは，時系列分析におけるラグの次数決定．
AICは「真のモデル」を選択する基準ではない．あくまで予測力の観点で（比較対象の中で相対的に最も）望ましいモデルを選択する基準である．
モデル選択後に回帰係数の検定を行わない方が良い．少なくとも，オーソドックスな t 検定をナイーブに適用するのはNG．

3.6.1 Find AIC

lm1 <- lm(Fertility ~ Examination, data = swiss)
lm2 <- lm(Fertility ~ Examination + Education, data = swiss)
lm3 <- lm(Fertility ~ Examination + Education + Agriculture + Catholic, data = swiss)
AIC(lm1, lm2, lm3)

##     df      AIC
## lm1  3 350.3525
## lm2  4 344.6292
## lm3  6 332.4121

3.6.2 Stepwise model selection with AIC

考え得る説明変数をすべて含めたモデル (full model) からスタートして stepwise で最適なモデルを選択することもできる．

step 関数を使うと AIC を規準としてモデル選択が行われる．

lm_full <- lm(Fertility ~ Examination + Education + Agriculture + Catholic + Infant.Mortality, data = swiss)
step(lm_full, trace = 0)  # trace = 1 にすると計算過程が出力される

## 
## Call:
## lm(formula = Fertility ~ Education + Agriculture + Catholic + 
##     Infant.Mortality, data = swiss)
## 
## Coefficients:
##      (Intercept)         Education       Agriculture          Catholic  
##          62.1013           -0.9803           -0.1546            0.1247  
## Infant.Mortality  
##           1.0784

4 Topics in Econometrics

4.1 Properties of fitted values and residuals

再掲：一変数の回帰モデルの係数パラメタ推定

y_bar <- mean(swiss$Fertility)  # Y bar
x_bar <- mean(swiss$Examination)  # X bar
s_xx <- sum((swiss$Examination - x_bar)^2)  # S_XX
s_xy <- sum((swiss$Examination - x_bar) * (swiss$Fertility - y_bar))  # S_XX
beta_hat <- s_xy / s_xx
alpha_hat <- y_bar - beta_hat * x_bar

アウトカムの当てはめ値 $\hat{Y}_i = \hat{\alpha} + \hat{\beta} X_i$ の平均はアウトカムの平均と一致する．

\[ \bar{\hat{Y}} = \bar{Y} \]

y_pred <- alpha_hat + beta_hat * swiss$Examination  # y hat
c(mean(y_pred), y_bar)  # should be matched

## [1] 70.14255 70.14255

残差（誤差変数の当てはめ値, $\hat{u} = Y - \hat{Y}$）に関する幾つかの性質．

\[ \sum_i \hat{u}_i = 0, \quad \sum_i \hat{u}_i X_i = 0 \]

resid_swiss <- swiss$Fertility - y_pred  # residuals
sum(resid_swiss)  # should be zero

## [1] 1.065814e-13

sum(resid_swiss * swiss$Examination)  # eX ... should be zero

## [1] 1.58451e-12

\[ S_{X\hat{u}} = \sum_i (X - \bar{X}) (\hat{u} - \bar{\hat{u}}) = 0, \quad S_{Y\hat{u}} = \sum_i (Y - \bar{Y}) (\hat{u} - \bar{\hat{u}}) = 0 \]

sum((swiss$Examination - x_bar) * (resid_swiss - mean(resid_swiss)))  # S_Xe ... should be zero

## [1] -1.501022e-13

sum((y_pred - y_bar) * (resid_swiss - mean(resid_swiss)))  # S_Ye ... should be zero

## [1] -8.354428e-14

4.2 Joint hypothesis test / 結合有意性検定

「Education と Examination のどちらの回帰パラメタも 0 に等しい」という帰無仮説を検定するためには Joint test を行う．

\[ H_0: \beta_{Educ} = \beta_{Exam} = 0 \]

制約の本数 $G$ を用いて，検定統計量 $F$ は次のように計算される．

\[ F = \chi^2 \frac{1}{G} = \frac{(Q_R - Q) / G}{Q / (n-k-1)} \sim F(G, n-k-1) \]

$Q$ は制約なしのモデル（対立仮説に対応．つまり，Education と Examination の2つの変数を説明変数として含むモデル）の残差二乗和を表す．

\[ Fertility = \alpha + \beta_{Educ} Educ + \beta_{Exam} Exam+ u, \quad Q = \sum_i \hat{u}^2_i \]

$Q_R$ は制約あり (restricted) のモデル（帰無仮説に対応．つまり，Education と Examination の2つの変数を説明変数として含まないモデル）の残差二乗和を表す．

\[ Fertility = \alpha + v, \quad Q_R = \sum_i \hat{v}^2_i \]

lm_swiss_unrestricted <- lm(Fertility ~ Examination + Education, data = swiss)
lm_swiss_restricted <- lm(Fertility ~ 1, data = swiss)  # 定数項のみのモデル
resid_unrestricted <- sum(lm_swiss_unrestricted$resid^2)
resid_restricted <- sum(lm_swiss_restricted$resid^2)

計算した $Q, Q_R$ を用いて $F$ 値を計算する．

f_value <- ((resid_restricted - resid_unrestricted) / 2) / (resid_unrestricted / (nrow(swiss) - 2 - 1))
f_value

## [1] 22.48799

1 - pf(q = f_value, df1 = 2, df2 = nrow(swiss) - 2 - 1)  # p value

## [1] 1.8705e-07

臨界値と比較する場合．

qf(p = 0.05, df1 = 2, df2 = nrow(swiss) - 2 - 1, lower.tail = F)  # critical value

## [1] 3.209278

# qf(p = 0.95, df1 = 2, df2 = nrow(swiss) - 2 - 1, lower.tail = T)  # same as above

$F \approx 22.5 > \mbox{critical value}$ なので，帰無仮説は5%有意水準で棄却される．すなわち，Education と Examination という2つの説明変数のうち少なくともいずれか一方の回帰パラメタは0ではないと結論付けられる．

4.2.1 Using `car::linearHypothesis` function

car というパッケージの linearHypothesis という関数で検定できる．「car::linearHypothesis」のようにパッケージ名と関数名を :: で繋いで使用すれば，当該パッケージを読み込まなくても関数を使える．

lm_swiss <- lm(Fertility ~ Examination + Education, data = swiss)
summary(lm_swiss)

## 
## Call:
## lm(formula = Fertility ~ Examination + Education, data = swiss)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.9935  -6.8894  -0.3621   7.1640  19.2634 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  85.2533     3.0855  27.630   <2e-16 ***
## Examination  -0.5572     0.2319  -2.402   0.0206 *  
## Education    -0.5395     0.1924  -2.803   0.0075 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.982 on 44 degrees of freedom
## Multiple R-squared:  0.5055, Adjusted R-squared:  0.483 
## F-statistic: 22.49 on 2 and 44 DF,  p-value: 1.87e-07

car::linearHypothesis(lm_swiss, "Examination = 0")

## Linear hypothesis test
## 
## Hypothesis:
## Examination = 0
## 
## Model 1: restricted model
## Model 2: Fertility ~ Examination + Education
## 
##   Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
## 1     45 4015.2                              
## 2     44 3549.6  1    465.63 5.7718 0.02057 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

p値はどちらも 0.0206 となっていることが確認できる．

\[ H_0: \beta_{Educ} = \beta_{Exam} = 0 \]

car::linearHypothesis(lm_swiss, c("Examination = 0", "Education = 0"))

（出力は省略）

\[ H_0: \beta_{Educ} + \beta_{Exam} = 1 \]

car::linearHypothesis(lm_swiss, "Examination + Education = 1")

（出力は省略）

4.2.2 Advanced: ANOVA

分散分析について知っている履修者は次の2つの分析の関係を検討してみよう．

被説明変数をカテゴリー変数ダミーに回帰する線形回帰モデル (lm 関数で推定)
一元配置分散分析 (aov 関数で推定)

summary(lm(Fertility ~ region, data = swiss2))$fstatistic

##    value    numdf    dendf 
## 29.07634  5.00000 41.00000

anova(aov(Fertility ~ region, data = swiss2))$`F value`

## [1] 29.07634       NA

4.3 White standard error / ホワイトの標準誤差

残差 $\hat{u} = Y - \hat{Y}$ を用いた不均一分散に頑健な標準誤差（White によるオリジナルの方法；説明変数が1つの場合）．

\[ s.e.(\hat{\beta}) = \sqrt{ \frac{\sum_i (X_i - \bar{X})^2 \hat{u}_i^2}{\left[\sum_i (X_i - \bar{X})^2\right]^2} } \]

以下の回帰モデルの $\beta$ について考える．

\[ Fertility = \alpha + \beta \ Exam + u \]

まずは通常の標準誤差を計算．

summary(lm(Fertility ~ Examination, data = swiss))$coef

##              Estimate Std. Error   t value     Pr(>|t|)
## (Intercept) 86.818529  3.2576034 26.651043 3.353924e-29
## Examination -1.011317  0.1781971 -5.675275 9.450437e-07

White の標準誤差．

x <- swiss$Examination
lm_swiss <- lm(Fertility ~ Examination, data = swiss)
sqrt(sum((x - mean(x))^2 * lm_swiss$resid^2) / (sum((x - mean(x))^2))^2)  # White SE (original)

## [1] 0.1776713

誤差分散を不偏推定するために自由度を調整する場合（一般に HC1 として知られる）．

\[ s.e._{HC1}(\hat{\beta}) = \sqrt{ \frac{n}{n-(k+1)} \frac{\sum_i (X_i - \bar{X})^2 \hat{u}_i^2}{\left[\sum_i (X_i - \bar{X})^2\right]^2} } \]

sqrt(nrow(swiss) / (nrow(swiss) - 2) * sum((x - mean(x))^2 * lm_swiss$resid^2) / (sum((x - mean(x))^2))^2)

## [1] 0.1815766

fixest パッケージで推定する場合． feols 関数で vcov = "hetero" (誤差項の variance-covariance について heteroskedasticity を許容する，という意味) 引数を指定する．

# install.packages("fixest")
library(fixest)
feols(Fertility ~ Examination, data = swiss, vcov = "hetero")

## OLS estimation, Dep. Var.: Fertility
## Observations: 47 
## Standard-errors: Heteroskedasticity-robust 
##             Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept) 86.81853   3.175114 27.34344 < 2.2e-16 ***
## Examination -1.01132   0.181577 -5.56965 1.353e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 9.43462   Adj. R2: 0.404213

4.4 Cluster robust SE / クラスターロバスト標準誤差

swiss データにおいて，以下の6地域をクラスターとして利用する（再掲）．

A: ジュラカントン, B: ヴォー州, C: フリブール州, D: ネウシャテル州, E: ジュネーブ州, F: ヴァレー州

library(fixest)
table(swiss2$region)  # 再掲

## 
##  B  A  C  D  E  F 
## 19  6  5  6  3  8

feols(Fertility ~ Examination, data = swiss2, cluster = ~ region)

## OLS estimation, Dep. Var.: Fertility
## Observations: 47 
## Standard-errors: Clustered (region) 
##             Estimate Std. Error  t value   Pr(>|t|)    
## (Intercept) 86.81853   4.708134 18.44011 8.6278e-06 ***
## Examination -1.01132   0.281775 -3.58910 1.5723e-02 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 9.43462   Adj. R2: 0.404213

4.5 OLS vs. WLS Monte Carlo simulation / モンテカルロシミュレーション

誤差の不均一性をもたらすメカニズムが既知で当該要因を観測可能な場合は OLS よりも加重最小二乗法（WLS）の方が効率よく推定できる．

鹿野繁樹『新しい計量経済学』（日本評論社）第11章で行われたシミュレーションを再現．

OLS:

\[ Y_i = 1 + X_i + \sqrt{h} \tilde{u}_i \] \[ X_i \sim N(0, 1) , \quad \tilde{u}_i \sim \mbox{Uni}(0, 4 \sqrt{12}), \quad h_i \sim \mbox{Uni}(0, 1) \]

WLS:

\[ \frac{Y_i}{\sqrt{h_i}} = \frac{1}{\sqrt{h_i}} + \beta \frac{X_i}{\sqrt{h_i}} + u_i, \quad \tilde{Y}_i = I_i + \beta \tilde{X}_i + \tilde{u}_i \]

注：$E(u_i) \ne 0$ であるため，WLSで定数項を0にして推定してはいけない．

ols_matrix <- wls_matrix <- NULL
for (i in 1:1000) {
  set.seed(i)
  n <- 50
  x <- rnorm(n, 0, 1)
  u_tilde <- runif(n, 0, 4*sqrt(12))
  h <- runif(n, 0, 1)
  # OLS
  y <- 1 + x + u_tilde * h
  ols_matrix <- rbind(ols_matrix, lm(y ~ x)$coef)
  # WLS
  y_tilde <- y / sqrt(h)
  Ii <- 1 / sqrt(h)
  x_tilde <- x / sqrt(h)
  wls_matrix <- rbind(wls_matrix, lm(y_tilde ~ x_tilde)$coef)
}

$X$ の回帰係数の推定値の分布を観察する．

summary(wls_matrix[, 2]); summary(ols_matrix[, 2])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -0.7191  0.8522  1.0012  1.0113  1.1636  2.5938

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -0.3889  0.7342  0.9998  1.0167  1.3030  2.3599

hist(wls_matrix[, 2], breaks = 30)
hist(ols_matrix[, 2], breaks = 20, add = T, col = rgb(1,0,0, alpha = .2))
legend(x = "topright", col = c(8, 2), legend = c("WLS", "OLS"), lwd = 5)

WLS で推定すると推定値のばらつきが小さい．このことは，誤差の不均一性を正確にモデル化できる場合にはWLSを適用することで推定精度を改善できる可能性があることを示している．

5 Take home messages

回帰分析：summary(lm(formula = Y ~ X1 + X2, data = dataset))
説明変数の二乗は I(X1^2)
2つの説明変数の掛け算（交差項）は X1:X2
誤差項の不均一分散に頑健な標準誤差を計算するには feols(Y ~ X1 + X2, data = dataset, vcov = "hetero")

Data Analysis Using Statistical Packages: Regression Analysis

Sho Kuroda / 黒田翔

2024年3月 (Last update: 2024年7月)

1 Simple regression analysis / 単回帰分析

1.1 Estimation by hand

1.2 Estimation by numerical calculation / 数値計算による推定

1.2.1 Example of `optim`

1.2.2 In the case of regression

1.3 Coefficient of determination / 決定係数

1.4 Standard error and test / 標準誤差と検定

1.4.1 cf. Critical value approach / 臨界値を用いて検定を行う場合

2 Multiple regression analysis / 重回帰分析

2.1 Estimation by hand

2.2 Standard error / 標準誤差

3 Advanced method / 応用

3.1 Adding quardatic or interaction terms / 二次・交差項

3.2 Dummy variables ダミー変数

3.2.1 Relationship to the test of difference of means / 平均値の差の検定との関係

3.3 Categorical variables

3.3.1 Change the reference group

3.4 Multicollinearity / 多重共線性

3.4.1 Example of perfect multico. / 完全な多重共線性の例

3.4.2 Example of imperfect multico. / 不完全な多重共線性の例

3.4.3 Using VIF for diagnostics / VIF による診断

3.5 Comparing the size of the coefficients / 回帰係数の大きさを比較する

3.5.1 Comparing the change per one SD / 1標準偏差当たりの変化量の比較

3.5.2 Standardized partial regression coefficient / 標準偏回帰係数

3.6 Model selection / モデル選択

3.6.1 Find AIC

3.6.2 Stepwise model selection with AIC

4 Topics in Econometrics

4.1 Properties of fitted values and residuals

4.2 Joint hypothesis test / 結合有意性検定

4.2.1 Using `car::linearHypothesis` function

4.2.2 Advanced: ANOVA

4.3 White standard error / ホワイトの標準誤差

4.4 Cluster robust SE / クラスターロバスト標準誤差

4.5 OLS vs. WLS Monte Carlo simulation / モンテカルロシミュレーション

5 Take home messages

Data Analysis Using Statistical Packages: Regression Analysis

Sho Kuroda / 黒田翔

2024年3月 (Last update: 2024年7月)

1 Simple regression analysis / 単回帰分析

1.1 Estimation by hand

1.2 Estimation by numerical calculation / 数値計算による推定

1.2.1 Example of optim

1.2.2 In the case of regression

1.3 Coefficient of determination / 決定係数

1.4 Standard error and test / 標準誤差と検定

1.4.1 cf. Critical value approach / 臨界値を用いて検定を行う場合

2 Multiple regression analysis / 重回帰分析

2.1 Estimation by hand

2.2 Standard error / 標準誤差

3 Advanced method / 応用

3.1 Adding quardatic or interaction terms / 二次・交差項

3.2 Dummy variables ダミー変数

3.2.1 Relationship to the test of difference of means / 平均値の差の検定との関係

3.3 Categorical variables

3.3.1 Change the reference group

3.4 Multicollinearity / 多重共線性

3.4.1 Example of perfect multico. / 完全な多重共線性の例

3.4.2 Example of imperfect multico. / 不完全な多重共線性の例

3.4.3 Using VIF for diagnostics / VIF による診断

3.5 Comparing the size of the coefficients / 回帰係数の大きさを比較する

3.5.1 Comparing the change per one SD / 1標準偏差当たりの変化量の比較

3.5.2 Standardized partial regression coefficient / 標準偏回帰係数

3.6 Model selection / モデル選択

3.6.1 Find AIC

3.6.2 Stepwise model selection with AIC

4 Topics in Econometrics

4.1 Properties of fitted values and residuals

4.2 Joint hypothesis test / 結合有意性検定

4.2.1 Using car::linearHypothesis function

4.2.2 Advanced: ANOVA

4.3 White standard error / ホワイトの標準誤差

4.4 Cluster robust SE / クラスターロバスト標準誤差

4.5 OLS vs. WLS Monte Carlo simulation / モンテカルロシミュレーション

5 Take home messages

1.2.1 Example of `optim`

4.2.1 Using `car::linearHypothesis` function