1 Probability / 確率

1.1 Expectation and Variance of random variable

離散確率変数 \(X\) の期待値と分散は確率関数 \(f\) を用いて以下のように定義される．

\[ \mu = E(X) = \sum_j^k x_j f(x_j) \]

\[ \sigma^2 = V(X) = E[(X-E(X))^2] = \sum_j^k (x_j - E(X))^2 f(x_j) = E(X^2) - (E(X))^2 \]

連続確率変数の場合（ここでは \(f\) は確率密度関数）

\[ \mu = E(X) = \int_{-\infty}^{\infty} x f(x) dx \]

\[ \sigma^2 = V(X) = \int_{-\infty}^{\infty} (x - \mu)^2 f(x) dx\]

1.1.1 Example: Dice roll

例：サイコロの出目 \(X\) の期待値と分散 / Example of dice roll

\[ E(X) = \sum_{j=1}^6 x_j \frac{1}{6} = 1 \times \frac{1}{6} + 2 \times \frac{1}{6} + 3 \times \frac{1}{6} + 4 \times \frac{1}{6} + 5 \times \frac{1}{6} + 6 \times \frac{1}{6} = \frac{7}{2} \]

dice <- 1:6  # c(1, 2, 3, 4, 5, 6) と同じ
sum(dice * (1/6))  # definition of expectation
## [1] 3.5
mean(dice)  # mean
## [1] 3.5

\[ V(X) = \sum_{j=1}^6 (x_j - E(X))^2 \frac{1}{6} = \left( 1-\frac{7}{2} \right) ^2 \frac{1}{6} + \left( 2-\frac{7}{2} \right) ^2 \frac{1}{6} + \cdots + \left( 6-\frac{7}{2} \right) ^2 \frac{1}{6} \approx 2.917 \]

sum((dice - mean(dice))^2 * (1/6))  # variance: sum[(x-E(X))*f(x)]
## [1] 2.916667

なお，var(dice) は 3.5 となり，上記と一致しないが，これは var 関数が標本を引数として標本分散（不偏分散）を計算するためである． dice は標本ではなく母集団なので，この場合は母分散を計算しなければならない．

cf. 標本分散 \(s^2\) と母分散 \(\sigma^2\)

\[ s^2 = \frac{1}{n-1} \sum_i (x_i-\bar{x})^2, \quad \sigma^2 = \frac{1}{n} \sum_i (x_i-\mu)^2 \]

標本分散に \((n-1)/n\) を乗じて母分散に変換できる．

3.5 * (6-1) / 6  # approx. 2.917
## [1] 2.916667

\(V(X) = E(X^2) - (E(X))^2\) でもある．

sum(dice^2 * (1/6)) - sum(dice * (1/6))^2  # variance: E(X^2)-E(X)^2
## [1] 2.916667

1.2 Z score

標準正規分布の \(z = 1.96\) より左側の面積は確率 \(Pr(Z \le 1.96)\) である．

pnorm は与えられた \(z\) スコアに対応する確率を返す関数で，関数名は p (probability) + norm (normal distribution) から来ている．

pnorm(q = 1.96, mean = 0, sd = 1)
## [1] 0.9750021

その逆に，qnorm は与えられた確率に対応する \(z\) スコアを返す関数． 

qnorm(p = 0.975, mean = 0, sd = 1)
## [1] 1.959964
qnorm(p = 0.975)  # for N(0,1), the mean and sd arguments may be omitted 
## [1] 1.959964

確率密度． dnorm = d (density) + norm (normal distribution).

dnorm(x = 1.96, mean = 0, sd = 1)
## [1] 0.05844094

標準正規分布の確率密度 \(\phi\) を積分する（確率密度関数とX軸の間の面積を計算する）と，積分区間の間の値を取る確率が得られる．

\[ \Pr(-1.96 \le Z \le 1.96) = \int_{-1.96}^{1.96} \phi(z) dz \approx 0.95 \]

integrate(dnorm, -1.96, 1.96)  # 積分: 約 95%
## 0.9500042 with absolute error < 1e-11

この dnorm 関数は密度関数の描画にも使える．

curve(dnorm(x, mean = 0, sd = 1), xlim = c(-5, 5))

1.3 Generating random numbers / 乱数の生成

一様分布から乱数をドローするには，runif関数を使う．この関数名は r (random) + unif (uniform distribution) から来ている．

runif(n = サンプルサイズ, min = 下限, max = 上限) で乱数を生成．

runif(n = 5, min = 0, max = 1)
## [1] 0.29715894 0.27265571 0.46258601 0.67950573 0.09897916
runif(5, 0, 1)  # 引数の順番を守れば「n=」などを省略可
## [1] 0.04956765 0.46852701 0.44102886 0.23825383 0.42990961
runif(5)  # min = 0 かつ max = 1 の場合は省略可（デフォルトの設定）
## [1] 0.8443760 0.5234901 0.4353203 0.3651927 0.2525658

乱数をドローするたびに異なる値が出力される． 再現性のために乱数のシードを指定することができる．

set.seed(0)
runif(5)
## [1] 0.8966972 0.2655087 0.3721239 0.5728534 0.9082078
set.seed(0)
runif(5)
## [1] 0.8966972 0.2655087 0.3721239 0.5728534 0.9082078

正規乱数（正規分布からドローする乱数）は rnorm で生成できる．

rnorm(n = 5, mean = 0, sd = 10)
## [1] -8.356286 15.952808  3.295078 -8.204684  4.874291

1.3.1 Expectation and variance: Random variables from the uniform distribution

一様分布 \(U(a, b)\) から乱数を生成して標本平均と標本分散を計算し，それが理論値と一致するかを確認する．

\[ E(X) = \int_a^b x \frac{1}{b-a} dx = \frac{1}{b-a} \left[ \frac{x^2}{2} \right]_a^b = \frac{1}{b-a} \times \frac{b^2 - a^2}{2} = \frac{a+b}{2} \]

\[ V(X) = \int_a^b (x - \mu)^2 \frac{1}{b-a} dx = \frac{1}{b-a} \left[ \frac{ \left( x - \frac{a+b}{2} \right)^3}{3} \right]_a^b = \frac{(b - a)^2}{12} \]

\(X \sim U(-1, 1)\) とすると，期待値と分散は次の値となる．

\[ E(X) = \frac{a+b}{2} = \frac{-1+1}{2} = 0 \]

\[ V(X) = \frac{(b - a)^2}{12} = \frac{(1 - (-1))^2}{12} = \frac{1}{3} \]

「サンプルサイズ 100 で一様乱数を生成して平均と分散を計算する」を 1,000 回繰り返し，1,000個の平均と分散の分布を観察する． このようなシミュレーションを「モンテカルロ・シミュレーション」と呼ぶ．

  • for 文を使って，乱数生成，平均と分散の計算，記録を行う
    • 計算結果を格納するオブジェクトを runif_mean と runif_var という名前で作成．NULL を代入することで，空のオブジェクトを作成
    • c(直前までの計算結果のベクトル, 今回の計算結果) で append
runif_mean <- NULL  # 標本平均を格納しておくオブジェクト
runif_var <- NULL  # 標本分散を 〃
for (i in 1:1000) {
  set.seed(i)  # 再現性のため
  # 乱数を生成
  runif_i <- runif(n = 100, min = -1, max = 1)
  # 平均・分散を計算して上で作成したオブジェクトに格納 (append) 
  runif_mean <- c(runif_mean, mean(runif_i))
  runif_var <- c(runif_var, var(runif_i))
}
summary(runif_mean)
##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
## -1.842e-01 -4.050e-02 -4.024e-05 -1.030e-03  3.643e-02  2.225e-01
summary(runif_var)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.2413  0.3127  0.3327  0.3333  0.3545  0.4225

2 Central limit theorem / 中心極限定理

標本平均 \(\bar{x}\) は母平均 \(\mu\), 母分散 \(\sigma^2\), サンプルサイズ \(n\) で決まる正規分布に従う．

\[ \bar{x} \sim N \left( \mu, \frac{\sigma^2}{n} \right) \]

一様確率変数の標本平均の分布を求める．標本サイズが大きくなるほど正規分布に近づき，分布の分散が小さくなることが分かる．

2.1 \(n=1\)

これは乱数そのものの分布と等しい．

\[ \bar{x} = x_i, \quad x_i \sim U(0, 1) \]

x_bar <- NULL  # create empty (null) object
for (i in 1:10000) {
  x_bar_i <- mean(runif(n = 1, min = 0, max = 1))
  x_bar <- c(x_bar, x_bar_i)
}
hist(x_bar, breaks = 50, main = "Sample mean of uniform random variables with sample size n=1")

2.2 \(n=2\)

\[ \bar{x} = \frac{1}{2} \sum_{i=1}^{2} x_i, \quad x_i \sim U(0, 1) \]

x_bar <- NULL
for (i in 1:10000) {
  x_bar_i <- mean(runif(n = 2, min = 0, max = 1))
  x_bar <- c(x_bar, x_bar_i)
}
hist(x_bar, breaks = 50, main = "Sample mean of uniform random variables with sample size n=2")

2.3 \(n=10\)

\[ \bar{x} = \frac{1}{10} \sum_{i=1}^{10} x_i, \quad x_i \sim U(0, 1) \]

x_bar <- NULL
for (i in 1:10000) {
  x_bar_i <- mean(runif(n = 10, min = 0, max = 1))
  x_bar <- c(x_bar, x_bar_i)
}
hist(x_bar, breaks = 50, main = "Sample mean of uniform random variables with sample size n=10",
     xlim = c(0, 1))

2.4 For advanced students / 中級者向け

for 文を使わずに sapply 関数を使って次のように書くこともできる： 

hist(sapply(X = 1:10000, FUN = function (x) mean(runif(10)) ), breaks = 50)

（出力は省略）

3 Confidence interval / 信頼区間

3.1 Population mean with known \(\sigma\) / 母集団の標準偏差が既知の場合の母平均

標本平均が100，「母集団」の標準偏差が2，サンプルサイズが100の場合，信頼水準95%の信頼区間は次のように計算される．

\[ \bar{x} \pm z \frac{\sigma}{\sqrt{n}} = 100 \pm 1.96 \frac{2}{\sqrt{100}} = 100 \pm 0.392, \quad \mbox{95% CI}: [99.61, 100.39] \]

qnorm(p = 0.025, mean = 0, sd = 1)  # P(Z<z) = 0.025
## [1] -1.959964
qnorm(p = 0.975, mean = 0, sd = 1)  # P(Z<z) = 0.975
## [1] 1.959964
100 + qnorm(p = 0.025) * 2 / sqrt(100)  # lower limit
## [1] 99.60801
100 + qnorm(p = 0.975) * 2 / sqrt(100)  # upper limit
## [1] 100.392

3.2 Population mean with unknown \(\sigma\) / 母集団の標準偏差が未知の場合の母平均

標本平均が100，「標本」標準偏差が2，サンプルサイズが100の場合，信頼水準95%の信頼区間は次のようになる．

\[ \bar{x} \pm t \frac{s}{\sqrt{n}} = 100 \pm 1.98 \frac{2}{\sqrt{100}} = 100 \pm 0.396, \quad \mbox{95% CI}: [99.60, 100.40] \]

qt(p = 0.025, df = 100-1)  # P(T<t) = 0.025
## [1] -1.984217
100 + qt(p = 0.025, df = 100-1) * 2 / sqrt(100)  # lower limit
## [1] 99.60316
100 + qt(p = 0.975, df = 100-1) * 2 / sqrt(100)  # upper limit
## [1] 100.3968

3.2.1 Using t.test function

\(t\)-検定を行うt.testという関数は信頼区間をついでに出力してくれる．

まずは t.test 関数を使わない場合． 信頼水準95%とする．

wage <- c(1000, 1200, 1300, 1200, 1150, 1000, 1450, 1500, 1150, 1350)
mean(wage)  # mean 
## [1] 1230
sd(wage)  # standard deviation 
## [1] 170.2939
length(wage)  # sample size (length of the vector)
## [1] 10
mean(wage) + qt(p = 0.025, df = length(wage)-1) * sd(wage) / sqrt(length(wage))  # lower bound 
## [1] 1108.179
mean(wage) + qt(p = 0.975, df = length(wage)-1) * sd(wage) / sqrt(length(wage))  # upper bound 
## [1] 1351.821

同じ結果が t.test 関数で確認できる．

信頼水準 (\(1 - \alpha\)) はデフォルトで95%に設定されている． conf.level 引数を指定すればそれ以外の信頼水準に変更できる．

t.test(wage)  # default: 95% CI
## 
##  One Sample t-test
## 
## data:  wage
## t = 22.841, df = 9, p-value = 2.806e-09
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  1108.179 1351.821
## sample estimates:
## mean of x 
##      1230
t.test(wage, conf.level = 0.99)  # 99% CI
## 
##  One Sample t-test
## 
## data:  wage
## t = 22.841, df = 9, p-value = 2.806e-09
## alternative hypothesis: true mean is not equal to 0
## 99 percent confidence interval:
##  1054.991 1405.009
## sample estimates:
## mean of x 
##      1230

3.3 Population proportion / 母比率

標本の比率 \(p\)，サンプルサイズ \(n\)，z 値 (such that \(\Pr(Z \le z) = 1 - \alpha\)) より次のように母比率の信頼区間が計算される．

\[ p \pm z \sqrt{\frac{p (1-p)}{n}} \]

Example of Banerjee, Duflo, and Glennerster (BMJ 2010)

382人の子供のうち148人がワクチンの予防接種を受けた． 

n1 <- 148; n0 <- 382  # intervention B
p <- n1 / n0
p
## [1] 0.3874346
p + qnorm(p = .025) * sqrt(p * (1 - p) / n0)  # lower limit of 95% CI
## [1] 0.3385815
p + qnorm(p = .975) * sqrt(p * (1 - p) / n0)  # upper limit
## [1] 0.4362876

注：この数値は論文中で報告されている数値とは異なる． 地理的なブロック（村）を用いてクラスター無作為化対照試験が実行されたため，著者らはこのデータの階層性を考慮したより複雑な計算をしていると思われる．

4 Statistical hypothesis testing / 統計的仮説検定

4.1 One-sample \(t\) test for population mean / 1標本の母平均の \(t\) 検定

4.1.1 Calculation from the sample mean and sd of the data / データの標本平均・標準偏差から計算

賃金の母集団が正規分布に従っており，\(\sigma\) は未知，\(\bar{x} = 1230, s=170.3, n = 10, \alpha=0.05\) とする．

\[ H_0: \mu = 1100, \quad H_1: \mu \ne 1100 \]

\[ t = \frac{\bar{x} - \mu_0}{s/\sqrt{n}} = \frac{1230 - 1110}{170.3 / \sqrt{10}} \approx 2.41 \]

t_value <- (1230 - 1100) / (170.3 / sqrt(10))  # t value
t_value
## [1] 2.413952

\[ \Pr (T \le 2.41) = 0.98 , \quad p = 2 \times (1 - \Pr (T \le 2.41)) \approx 0.039 \]

pt(q = t_value, df = 10 - 1)  # Pr(T <= t)
## [1] 0.9805023
1 - pt(q = t_value, df = 10 - 1)  # one sided p value ... p = Pr(T > t)
## [1] 0.01949773
2 * (1 - pt(q = t_value, df = 10 - 1))  # two sided p value
## [1] 0.03899546

よって，\(H_0\) は棄却され \(H_1\) が採択される．

4.1.2 Calculation from the data itself / データそのものから計算

t.test(ベクトル, mu = 帰無仮説の仮説値) と指定．

wage <- c(1000, 1200, 1300, 1200, 1150, 1000, 1450, 1500, 1150, 1350)
mean(wage); sd(wage)
## [1] 1230
## [1] 170.2939
t.test(wage, mu = 1100)  # default: two sided ... H1 is [mu != 1100]
## 
##  One Sample t-test
## 
## data:  wage
## t = 2.414, df = 9, p-value = 0.03899
## alternative hypothesis: true mean is not equal to 1100
## 95 percent confidence interval:
##  1108.179 1351.821
## sample estimates:
## mean of x 
##      1230

片側検定をする場合は alternative = "greater" または = "less" 引数で指定．

t.test(wage, mu = 1100, alternative = "greater")  # one sided ... H1 is [mu > 1100]
## 
##  One Sample t-test
## 
## data:  wage
## t = 2.414, df = 9, p-value = 0.01949
## alternative hypothesis: true mean is greater than 1100
## 95 percent confidence interval:
##  1131.284      Inf
## sample estimates:
## mean of x 
##      1230

このように信頼区間も同時に計算してくれる．

4.2 Two-sample \(t\) test for population mean / 2標本の母平均の \(t\) 検定

4.2.1 Independent \(t\) test

4.2.1.1 Example of “wage” data

t.test(x = 1つ目のベクトル, y = 2つ目のベクトル) のように指定する．

wage_jp <- c(1000, 1200, 1300, 1200, 1150, 1000, 1450, 1500, 1150, 1350)  # Japan
wage_us <- c(900, 1300, 1200, 800, 1600, 850, 1000, 950)  # US
t.test(wage_jp, wage_us)  # default: Welch's t test (assuming unequal variance)
## 
##  Welch Two Sample t-test
## 
## data:  wage_jp and wage_us
## t = 1.4041, df = 11.205, p-value = 0.1874
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -87.42307 397.42307
## sample estimates:
## mean of x mean of y 
##      1230      1075
t.test(wage_jp, wage_us, var.equal = TRUE)  # t test (assuming equal variance)
## 
##  Two Sample t-test
## 
## data:  wage_jp and wage_us
## t = 1.479, df = 16, p-value = 0.1586
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -67.16377 377.16377
## sample estimates:
## mean of x mean of y 
##      1230      1075

4.2.1.2 Titanic data

データ読み込み．

titanic <- read.csv("https://raw.githubusercontent.com/kurodaecon/bs/main/data/titanic3_csv.csv")

性別によって年齢の平均値が異なるかどうかの検定．

t.test(x = titanic$age[titanic$sex == "female"], y = titanic$age[titanic$sex == "male"])
## 
##  Welch Two Sample t-test
## 
## data:  titanic$age[titanic$sex == "female"] and titanic$age[titanic$sex == "male"]
## t = -2.0497, df = 798.36, p-value = 0.04072
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -3.71593356 -0.08034711
## sample estimates:
## mean of x mean of y 
##  28.68709  30.58523

t.test(変数 ~ グループ, data = dataset) で，2つのグループ間の検定が可能．

t.test(age ~ sex, data = titanic)  # same as above
## 
##  Welch Two Sample t-test
## 
## data:  age by sex
## t = -2.0497, df = 798.36, p-value = 0.04072
## alternative hypothesis: true difference in means between group female and group male is not equal to 0
## 95 percent confidence interval:
##  -3.71593356 -0.08034711
## sample estimates:
## mean in group female   mean in group male 
##             28.68709             30.58523

\(p < 0.05\) なので，「年齢の平均値は男性と女性で等しい」という帰無仮説を有意水準5%で棄却できる． 男性の標本平均の方が高いという事実と合わせると，「男性の方が有意に平均年齢が高い」と言える．

サブサンプルに対して検定する場合やデータフレームを加工して作成した変数を用いて検定する場合はパイプ演算子を使うと便利（tidyverse が必要）．

  • t.test の直前のパイプ演算子より前の部分（データフレーム）が t.test 関数の第一引数にならないため，t.test 関数内で data = . を明示的に書く必要がある（. がパイプ演算子の前の部分を受け取る）
library(tidyverse)
titanic %>% 
  filter(pclass == 1) %>%
  t.test(age ~ sex, data = .)
## 
##  Welch Two Sample t-test
## 
## data:  age by sex
## t = -2.3283, df = 278.85, p-value = 0.02061
## alternative hypothesis: true difference in means between group female and group male is not equal to 0
## 95 percent confidence interval:
##  -7.3664537 -0.6169013
## sample estimates:
## mean in group female   mean in group male 
##             37.03759             41.02927

4.2.2 Paired (dependent) \(t\) test

対応のある t 検定は t.test の引数で paired = TRUE 引数を追加して計算する．

wage_w <- c(1000, 1200, 1300, 1200, 1150, 1000, 1450, 1500, 1150, 1350)  # wife
wage_h <- c(900, 1300, 1200, 800, 1600, 850, 1000, 950, 1200, 1400)  # husband
# length(wage_w); length(wage_h)  # = 10
t.test(wage_w, wage_h, paired = TRUE)
## 
##  Paired t-test
## 
## data:  wage_w and wage_h
## t = 1.1638, df = 9, p-value = 0.2744
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -103.8108  323.8108
## sample estimates:
## mean difference 
##             110
t.test(wage_w - wage_h)  # same as an ordinary "one-sample t test"
## 
##  One Sample t-test
## 
## data:  wage_w - wage_h
## t = 1.1638, df = 9, p-value = 0.2744
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -103.8108  323.8108
## sample estimates:
## mean of x 
##       110

上記の通り「要素ごとの差」という一標本の検定になっている．

4.3 Multiple testing problem: Simulation / 多重検定の問題に関するシミュレーション

「標準正規分布 \(N(0, 1)\) から \(n=100\) の乱数をドローして，その平均が 0 に等しいかどうかを検定する（\(H_0: \mu = 0\) とする一標本の母平均の検定）」という手続きを1,000回繰り返し，検定で計算されたp値の分布を観察する．

設定上，帰無仮説は棄却されないはずだが，現実には偶然に大きな値（または小さな値）ばかりがドローされることで標本平均が 0 から大きく乖離し，帰無仮説が棄却される場合がありうる．

t.test 関数で計算される p 値は t.test(...)$p.value で抽出できる．

t.test(x = 1:10)$p.value
## [1] 0.000278196
p_value_list <- NULL
for (i in 1:1000) {
  p_value_i <- t.test(rnorm(100, mean = 0, sd = 1))$p.value
  p_value_list <- c(p_value_list, p_value_i)
}
hist(p_value_list, breaks = 20)

母集団分布の母平均は \(\mu = 0\) であるからこの帰無仮説は棄却されないことが望まれるが，上記の通り5%程度の確率で誤って棄却されてしまう（第一種の過誤）．

この簡単なシミュレーションは有意水準 \(\alpha\) が何を表しているかを理解するのに役立つと同時に，多重検定の問題を理解するのにも役立つ．

4.3.1 Multiple tests performed independently

複数の標本（または複数の標本組み合わせ）に関する帰無仮説を独立して検定すると，そのうち少なくとも一つで誤って帰無仮説を棄却する確率が事前に設定した有意水準（5%など）を上回ってしまう．

1,000 の標本のうち「1つ」について \(H_0: \mu = 0\) を検定した場合，帰無仮説が誤って棄却されてしまう割合は以下の通り（概ね5%）．

table(p_value_list <= 0.05)
## 
## FALSE  TRUE 
##   941    59

1,000 の標本のうち「2つ」について \(H_0: \mu = 0\) を検定した場合，少なくともどちらか一方で帰無仮説が誤って棄却されてしまう割合は以下の通り（概ね10%）．

\[ 1 - \mbox{(2回とも棄却されない確率)} = 1 - (1 - 0.05)^2 = 0.0975 \]

  • sapply(X = ベクトル, FUN = 関数) ：ベクトルの要素を1つずつ関数で評価しその結果をベクトルで返す
  • sample(ベクトル, 個数) ：指定した個数分だけベクトルからランダムにサンプリング
  • any(論理値ベクトル) ：論理値ベクトルの要素の1つ以上が TRUE の場合に TRUE を返す
table(sapply(X = 1:1000, FUN = function (x) any(sample(p_value_list, 2) <= 0.05)))
## 
## FALSE  TRUE 
##   902    98

1,000 の標本のうち「5つ」について \(H_0: \mu = 0\) を検定した場合，少なくとも1つで帰無仮説が誤って棄却されてしまう割合は以下の通り（概ね23%）．

\[ 1 - (1 - 0.05)^5 \approx 0.226 \]

table(sapply(X = 1:1000, FUN = function (x) any(sample(p_value_list, 5) <= 0.05)))
## 
## FALSE  TRUE 
##   754   246

.