1. モンテ・カルロシミュレーション
2. 偏回帰係数(傾き)の理解
3. 標準誤差の理解
4. 信頼区間の理解
References
・回帰分析のしくみ(特に信頼区間の意味)を理解するために、簡単なモンテ・カルロシミュレーションを行う。
・シミュレーションでは、自分で母数を設定し、データを人工的に作り出す。
・特徴を知りたい分析手法(ここではOLS)を人工的に作り出したデータに当てはめ、母数を期待通りに推測できるかどうか確認する。
・ここでは、単回帰を例にシミュレーションを行う。 このシミュレーションを行う主な目的は以下の 3 つ。
単回帰モデルは
\[y_i \sim \mathrm{N}(\alpha_1 + \alpha_2 x_i, \sigma^2)・・・i = 1, 2, \dots, n\]
と表すことができる。
設定する母数は次の3つである。
alpha1 <- 10 # 切片の値は 10
alpha2 <- 3 # 傾き(独立変数の係数)は 3
sigma <- 6 # 誤差の標準偏差は 6
次に、単回帰モデルが想定するデータ生成過程に従って、データを生成する。
・サンプルサイズ = 100 と指定
n <- 100
・同一のサンプル結果を得られるよう特定のノイズを設定
set.seed(2016)
・最大値10の一様分布から長さ n の x ベクトルを無作為抽出
x <- runif(n, 0, 10)
y の値を指定
y <- rnorm(n, alpha1 + alpha2 * x, sd = sigma)
・設定に従って正規分布から抽出した x と y をデータフレームに入れる
df <- data.frame(y=y, x=x)
head(df) # データの文頭を表示
y x
1 14.452955 1.801636
2 8.950190 1.429437
3 20.157514 8.416465
4 7.245589 1.335746
5 25.262810 4.775025
6 3.883445 1.212584
tail(df) # データの終わりを表示
y x
95 38.43990 9.8589020
96 15.22791 0.9887315
97 38.14979 7.9868813
98 24.89843 5.2812564
99 16.12850 0.4392445
100 30.47304 3.5779110
・散布図にして直線を当てはめる
library("ggplot2")
scat <- ggplot(df, aes(x, y)) +
geom_point() +
geom_smooth(method = "lm")
print(scat)
切片:\(\alpha_1=10\)
傾き:\(\alpha_2=3\)
標準偏差:\(\sigma=6\)
このサンプルから推定値を計算する:
reg1 <- lm(y ~ x, data = df)
summary(reg1)
Call:
lm(formula = y ~ x, data = df)
Residuals:
Min 1Q Median 3Q Max
-15.0923 -3.3533 0.2759 3.4820 16.4403
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 11.4369 1.2201 9.374 2.77e-15 ***
x 2.8293 0.2301 12.297 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6.534 on 98 degrees of freedom
Multiple R-squared: 0.6068, Adjusted R-squared: 0.6027
F-statistic: 151.2 on 1 and 98 DF, p-value: < 2.2e-16
切片: Intercept = 11.44
傾き: coef.est = 2.83
標準誤差: residual stanard error = 6.53
それぞれの母数 (10, 3, 6) に近い値が得られた。
・confint(reg1) とコマンドを入力すると係数の95パーセント信頼区間が得られる。
confint(reg1)
2.5 % 97.5 %
(Intercept) 9.015618 13.858253
x 2.372720 3.285918
・この 95 パーセント信頼区間は最初に設定した母数 \(\alpha_1=10, \alpha_2=3\) を含んでいる。
・以上の過程を、母数とサンプルサイズは変えずに何度も繰り返す。 ・同じコードを何度も繰り返し使うのは面倒なので、関数にする。
sim_ols1 <- function(alpha, sigma, n = 100, trials = 10000, x.rng = c(0, 10)) {
## Function to simulate simple OLS regressions
## Arguments:
## alpha: vector of coefficient parameters
## sigma: standard deviation of the error
## n: sample size, default n=100
## trials: the number of simulation trials
## default trials=10000
## x.rng: range of the explanatory variable x
## default x.rng=c(0,10)
## Return:
## df: data frame containing
## (1) parameter estimates,
## (2) se of coefficient estimates, and
## (3) 95% CI of coefficient estimates
## for each trial of the simulation.
## 結果を保存するために df という名前を付けてデータフレームを作る
empty <- rep(NA, trials)
df <- data.frame( # データフレームには次の 9 つの変数を含める
a1 = empty,
a1.se = empty,
a1.lower = empty,
a1.upper = empty,
a2 = empty,
a2.se = empty,
a2.lower = empty,
a2.upper = empty,
sigma.hat = empty
)
for (i in 1:trials) { ## loop for trials
## generate data
x <- runif(n, x.rng[1], x.rng[2])
y <- rnorm(n, alpha[1] + alpha[2] * x, sigma)
## fit OLS
fit <- lm(y ~ x)
## sum of squared residuals
sigma.hat <- summary(fit)$sigma
## coefficient estimates
a1 <- coef(fit)[1]
a2 <- coef(fit)[2]
## standard errors of coefficient estimates
a1.se <- sqrt(summary(fit)$cov.unscaled[1,1]) * sigma.hat
a2.se <- sqrt(summary(fit)$cov.unscaled[2,2]) * sigma.hat
## 95% CI
a1.ci95 <- confint(fit)[1,]
a2.ci95 <- confint(fit)[2,]
## save the results
df[i,] <- c(a1, a1.se, a1.ci95,
a2, a2.se, a2.ci95,
sigma.hat)
}
## return the data frame
return(df)
}
・この関数を利用して、実際にシミュレーションを行う。
・母数とサンプルサイズ(ここでは規定値のn=100を利用)を指定し、データの生成と回帰式の推定を1,000回繰り返してみる。
alpha1 <- 10; alpha2 <- 3; sigma = 6
sim.1 <- sim_ols1(alpha = c(alpha1, alpha2), sigma = sigma, trials = 1000)
head(sim.1)
a1 a1.se a1.lower a1.upper a2 a2.se a2.lower
1 9.641121 1.315498 7.030558 12.25168 3.076103 0.2273122 2.625009
2 10.506001 1.149149 8.225553 12.78645 2.897796 0.1895133 2.521713
3 8.806736 1.391324 6.045699 11.56777 3.182437 0.2242534 2.737413
4 9.911978 1.091111 7.746705 12.07725 2.970146 0.2051283 2.563075
5 9.528927 1.290571 6.967831 12.09002 2.957848 0.2121729 2.536798
6 10.168580 1.246947 7.694055 12.64310 3.104525 0.2120863 2.683646
a2.upper sigma.hat
1 3.527197 6.634995
2 3.273879 5.383881
3 3.627460 6.163088
4 3.377216 5.709748
5 3.378898 6.329047
6 3.525403 6.410671
これで、シミュレーションで得られた数字はsim.1に保存された。
・説明変数 x の係数の推定値 \(\alpha_2\) の分布を確認してみよう。
・私たちは母数の傾き \(\alpha_2 =3\) であることを知っているが、OLSはこの数値をうまく推定してくれただろうか?
hist.a2 <- ggplot(sim.1, aes(x = a2, y=..density..)) + geom_histogram(binwidth = .1)
hist.a2 <- hist.a2 + labs(y = "Probability Density",
title = "Simulation Result: Distribution of alpha_2") +
geom_vline(xintercept = 3, color = "green")
print(hist.a2)
・上記のヒストグラムが示すように、OLSは平均すると、母数として設定した\(\alpha_2 =3\)をうまく推定している。
・これを不偏性 (unbiasedness)という。
・しかし、常に正しい推定をするわけではない。
・係数の値を過小推定することもあれば、過大推定することもある。
・実際の分析では、データセットが1つしかないのが普通。
・→ 自分のデータ分析が係数を「正しく」推定しているとは限らならい。
・そのために、推定の不確実性を明示することが求められる。
・\(\alpha_2\) の標準誤差をヒストグラムで表示する。
hist.se <- ggplot(sim.1, aes(x = a2.se, y =..density..)) +
geom_histogram(binwidth = .01) +
labs(x = "Standard Error: alpha_2", y = "Probability Density", title = "Simulation Result: Distribution of Standard Error alpha_2")
print(hist.se)
・標準誤差 (se) 自体が推定量なので、その値はサンプルごとに異なる(つまり、分布する)。
・標準誤差 (se) は次の式で求められる。
\[\frac{\hat{\alpha} - \alpha}{\mathrm{se}}\]
・これは自由度 \(n-k\)(ここでは、100 - 2 = 98)の \(t\) 分布に従うはず。
・ここで求めた標準誤差 (se) をヒストグラムにして、自由度98の \(t\) 分布の確率密度曲線を重ね書きすると。
library("dplyr")
sim.1 <- mutate(sim.1, a2.t = (a2 - alpha2) / a2.se)
## t分布の確率密度を求め、true.tという名前のデータフレームに保存する
true.t <- data.frame(x = seq(-4, 4, length = 100))
true.t <- mutate(true.t, density = dt(x, df = 98))
hist.t <- ggplot(sim.1, aes(x = a2.t, y =..density..)) +
geom_histogram(binwidth = 0.5) +
geom_line(data = true.t, aes(x = x, y = density), color = "violetred")
hist.t <- hist.t +
labs(x = expression(frac(hat(alpha) - alpha, "se")),
y="Probability Density", title="Standardied alpha_2 and t-distributio with df=98")
print(hist.t)