1. モンテ・カルロシミュレーション
2. 偏回帰係数（傾き）の理解
3. 標準誤差の理解
4. 信頼区間の理解
References

1. モンテ・カルロシミュレーション

・回帰分析のしくみ（特に信頼区間の意味）を理解するために、簡単なモンテ・カルロシミュレーションを行う。
・シミュレーションでは、自分で母数を設定し、データを人工的に作り出す。
・特徴を知りたい分析手法（ここではOLS）を人工的に作り出したデータに当てはめ、母数を期待通りに推測できるかどうか確認する。
・ここでは、単回帰を例にシミュレーションを行う。 このシミュレーションを行う主な目的は以下の 3 つ。

  1. 最小二乗法が想定するデータ生成過程 (data generating process) を理解する
  2. 最小二乗推定量（OLSE: Ordinary Least Squares Estimator）の基本的な性質を理解する
  3. 信頼区間 (Confidence Interval) の意味を理解する

単回帰モデルは

\[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)

最初に母数として次の 3 つの値を設定した:　

切片：\(\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\) を含んでいる。

→ここで得られた信頼区間は 100% の確率で母数を含む

・以上の過程を、母数とサンプルサイズは変えずに何度も繰り返す。 ・同じコードを何度も繰り返し使うのは面倒なので、関数にする。

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に保存された。


2. 偏回帰係数（傾き）の理解

・説明変数 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つしかないのが普通。
・→ 自分のデータ分析が係数を「正しく」推定しているとは限らならい。
・そのために、推定の不確実性を明示することが求められる。

3. 標準誤差の理解

・\(\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)