1 Simple regression analysis / 単回帰分析

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

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

データは，R の built-in データセットの一つである swiss． 1988年のスイスにおける47の province レベルの社会経済変数を含む．

Fertility: 標準化された出生率
Agriculture: 男性の農業従事者割合
Examination: 徴兵検査における最高得点取得者割合
Education: 徴兵対象者のうち初等教育より高い教育を受けた人の割合

1.1 Estimation with `lm` function

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

これは以下のように回帰式が推定されたことを意味している．

\[ \hat{Fertility} = 86.8 - 1.01 Exam \]

Exam の回帰係数は「Exam が1単位高いと Fertility は 1.01 単位低い傾向にある」と解釈される．

「Exam を1単位引き上げると Fertility は 1.01 単位減少する」と解釈することはできない．この回帰係数はあくまで両変数の相関関係を記述しているだけであり因果効果を示すものではない．

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.2 Estimation by hand

回帰モデル

\[ Y_i = \alpha + \beta X_i + u_i \]

の係数パラメタ（$\alpha, \beta$）を最小二乗法（OLS）によって求めよう．

最小二乗法とは，誤差 $u$ （正確には残差）の二乗和を最小化するように係数パラメタを求める方法である．誤差の二乗和 $Q$ は以下で与えられる．

\[ Q = \sum_{i} \hat{u}_i^2 = \sum_{i} (Y_i - \alpha - \beta X_i)^2 \]

これは $\hat \alpha, \hat \beta$ のいずれについても下に凸の関数であるので，$Q$ をそれぞれで微分して 0 とおく．

\[ \frac{\partial Q}{\partial \alpha} = -2 \sum_i (Y_i - \alpha - \beta X_i) = 0 \]

\[ \frac{\partial Q}{\partial \beta} = -2 \sum_i x_i (Y_i - \alpha - \beta X_i) = 0 \]

上の連立方程式を解くことによって，最小二乗推定量は次のように与えられる．

\[ \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.3 Estimation by numerical calculation / 数値計算による推定

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

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

1.3.1 Example of `optim`

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

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

plot(function(x) x^2-4*x+7, xlim = c(-2,6))

最小値を取るような $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.3.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.4 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% を説明することができる，という意味．

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

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

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

1.5 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.5.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

Examination の係数は「Education の水準を固定したもとで Examination が1単位高い値を取ると Fertility は 0.557 単位低い値を取る傾向にある」と解釈される．

「単回帰分析では因果関係が推定できないが，重回帰分析では因果関係が推定できる」という主張は一般には誤り．

重回帰分析によって因果効果を推定するためには，処置変数の条件付き独立（かなり大雑把に言うと，交絡変数がすべてコントロール変数として正確な関数形でモデルに含められていること）が成立していなければならず，この仮定を社会経済データを用いた実証分析において満たすことは容易ではない．

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)$coef

## (Intercept) Examination        Edu2 
## 83.22071855 -0.66068646 -0.01034904

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

lm(Fertility ~ ExamEdu, data = swiss2)$coef

## (Intercept)     ExamEdu 
## 75.73211009 -0.02394142

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

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

##    (Intercept)    Examination I(Education^2) 
##    83.22071855    -0.66068646    -0.01034904

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

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

##           (Intercept) Examination:Education 
##           75.73211009           -0.02394142

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

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

##           (Intercept)           Examination             Education 
##          87.178103515          -0.625731412          -0.807551979 
## Examination:Education 
##           0.009201476

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

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

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

## (Intercept) Examination 
##  4.49251492 -0.01573623

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.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 と非常に相関の高い変数（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標準偏差当たりの変化量の比較

c(-1.195 * sd(swiss$Examination), -0.094 * sd(swiss$Agriculture))

## [1] -9.533571 -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 の分散をより説明する」と言える．

なお，t 値と p 値は標準化の前後で不変．

スケーリング（分散を 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 Take home messages

回帰分析：summary(lm(formula = Y ~ X1 + X2, data = dataset))
説明変数の二乗は I(X1^2)
2つの説明変数の掛け算（交差項）は X1:X2

Data Analysis Using Statistical Packages: Regression Analysis

Sho Kuroda / 黒田翔

2024年3月 (Last update: 2025年4月)

1 Simple regression analysis / 単回帰分析

1.1 Estimation with `lm` function

1.2 Estimation by hand

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

1.3.1 Example of `optim`

1.3.2 In the case of regression

1.4 Coefficient of determination / 決定係数

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

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

2 Multiple regression analysis / 重回帰分析

3 Advanced method / 応用

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

3.2 Dummy variables ダミー変数

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 Take home messages

Data Analysis Using Statistical Packages: Regression Analysis

Sho Kuroda / 黒田翔

2024年3月 (Last update: 2025年4月)

1 Simple regression analysis / 単回帰分析

1.1 Estimation with lm function

1.2 Estimation by hand

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

1.3.1 Example of optim

1.3.2 In the case of regression

1.4 Coefficient of determination / 決定係数

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

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

2 Multiple regression analysis / 重回帰分析

3 Advanced method / 応用

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

3.2 Dummy variables ダミー変数

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 Take home messages

1.1 Estimation with `lm` function

1.3.1 Example of `optim`