3 第5章基礎統計学復習

経営学部の1回生がおおよそ学んでいるであろう統計学の内容を、Rコードとともに復習します。テキストの内容に加えて、コードの効率化や関数化についても解説します。

# コード5-1
pacman::p_load(tidyverse, readxl, pwr, knitr)

まずサイコロを作ります。 1から6の数値からなるベクトルdiceを作ります。つぎにsample()関数を用いて、diceの中から1つの数をランダムに取り出します。

sample()関数の引数は以下の通りです。

x: 抽出元のベクトル
size: 抽出する数
replace: 置換を許可するかどうか

# コード5-2
set.seed(442)
# サイコロ
dice <- 1:6
# サイコロを1回振る
d <- sample(dice, size = 1)
# 結果の表示
d

[1] 6

同じコードを3回以上書いているなと思ったら、ループ処理や関数化を検討しましょう。 for()関数を使ったり、サイコロを投げる回数を引数にした関数を作成したりする方法があります。

d1 <- sample(dice, size = 10, replace = TRUE)
d2 <- sample(dice, size = 10, replace = TRUE)
d3 <- sample(dice, size = 10, replace = TRUE)

をまとめて書くと、次のようになります。

d <- list() # 空のリストを作成
for (i in 1:3) { # ループ処理
  d[[i]] <- sample(dice, size = 10, replace = TRUE)
}
d # 結果の表示

[[1]]
 [1] 4 4 3 2 6 3 5 2 6 5

[[2]]
 [1] 3 2 6 6 5 4 6 4 3 1

[[3]]
 [1] 6 1 1 5 2 1 4 5 2 2

これをさらに発展させて、後で変更する可能性のある変数である試行回数やサイコロを振る回数を別に設定しておいて、コードを読みやすく、また変更しやすくします。

trial_count <- 3 # 試行回数
sample_sizes <- 10 # サイコロを振る回数

# 空の数列を作成
d_mean <- numeric(trial_count) 
for (i in 1:trial_count) {
  d_mean[i] <- mean(sample(dice, size = sample_sizes, replace = TRUE))
}
d_mean

[1] 3.2 3.3 4.0

試行回数trialsとサイコロを振る回数sample_sizesを引数にした関数を作成すると、さらに便利になります。ここでは、自作関数dice_mean()を定義してみましょう。デフォルト値として、試行回数を3回、サイコロを振る回数を10回に設定しています。

dice_mean <- function(trial_count = 3, sample_sizes = 10) {
  d_mean <- numeric(trial_count) # 空のリストを作成
  for (i in 1:trial_count) {
    d_mean[i] <- mean(sample(dice, size = sample_sizes, replace = TRUE))
  }
  return(d_mean)
}

dice_mean() # デフォルト値で実行

[1] 3.4 2.5 4.4

dice_mean(trial_count = 5, sample_sizes = 20) # 引数を指定して実行

[1] 3.15 3.55 3.65 3.50 3.25

dice_mean(5, 20) # 省略形でもOK

[1] 3.80 4.30 2.75 3.70 3.05

この関数を使って、サイコロを10回、100回、1,000回振ったときの平均をそれぞれ3回ずつ試行してみましょう。

set.seed(352)
size = c(10, 100, 1000)
d_mean <- list()

for (i in seq_along(size)) {
  d_mean[[i]] <- dice_mean(trial_count = 3, sample_sizes = size[i])
}

d_mean_mat <- do.call(rbind, d_mean)

# 列名を付与
colnames(d_mean_mat) <- paste0("試行", 1:ncol(d_mean_mat))
# 行名を付与
rownames(d_mean_mat) <- paste0("サイコロを", size, "回振る")
# 表として出力
knitr::kable(
  d_mean_mat,
  caption = "サイコロの標本平均（各サイズ3回ずつ）",
  align = "ccc"
)

サイコロの標本平均（各サイズ3回ずつ）
	試行1	試行2	試行3
サイコロを10回振る	3.200	3.000	2.400
サイコロを100回振る	3.440	3.340	3.560
サイコロを1000回振る	3.433	3.446	3.554

tidyverseのpurrrパッケージを使うと、さらに簡潔に書けます。

set.seed(352)
trial_count <- 3 # 試行回数
sizes <- c(10, 100, 1000) # サイコロを振る回数
# map_dfr()でデータフレームを直接作成
d_mean_df <- map_dfr(sizes, function(n) {
  # replicate()で3回繰り返し試行    
  results <- replicate(
    trial_count, 
    mean(sample(dice, size = n, replace = TRUE))
    )
  set_names(results, paste0("試行", 1:trial_count))
  })

result <- d_mean_df |> 
  mutate(条件 = paste0("サイコロを", sizes, "回振る")) |> 
  relocate(条件) # 「条件」列を一番左へ移動

# 出力
kable(
  result,
  caption = "サイコロの標本平均（各サイズ3回ずつ）",
  align = "cccc"
)

サイコロの標本平均（各サイズ3回ずつ）
条件	試行1	試行2	試行3
サイコロを10回振る	3.200	3.000	2.400
サイコロを100回振る	3.440	3.340	3.560
サイコロを1000回振る	3.433	3.446	3.554

3.1 区間推定

いま、標準正規分布に従う母集団から抽出した無作為標本Z_1, Z_2, ..., Z_nを考えます。標準正規分布の確率分布関数\phi(z)を用いると、次のように表せます。 P \left( -z_{\frac{\alpha}{2}} \leq \frac{\bar{Z} - \mu}{\frac{\sigma}{\sqrt{n}}} \leq z_{\frac{\alpha}{2}} \right) = 1 - \alpha

3.1.1 電球の寿命データで区間推定

例として、次のような白熱電球の寿命データが手元にあるものとしましょう。ここで、以下のことが分かっているものとします。

平均寿命が1800時間である。
寿命の標準偏差$$は180時間である。
平均寿命は正規分布に従う。

図にすると次のようになります。

averate <- 1800
standard_dev <- 180
rnorm(1000, mean = averate, sd = standard_dev) |>
  as_tibble() |>
  ggplot(aes(x = value)) +
  stat_function(
    fun = dnorm,
    args = list(mean = averate, sd = standard_dev),
    color = "red",
    size = 1
    ) +
  theme_minimal()

この分布から抽出したと想定される16個の電球の寿命データが以下の通りです。

# コード5-5
bulb <- c(# 白熱電球の寿命データ
    1939.6, 1680.3, 1982.1, 2215.6,
    2092.5, 1928.9, 2003.8, 1955.5,
    1800.1, 1659.5, 2066.2, 2107.2,
    2085.5, 1878.6, 2007.6, 1816.1
  )
mean(bulb)

[1] 1951.194

sd(bulb)

[1] 154.4567

先ほどの分布にこのデータのヒストグラムを重ねると、次のようになります。

# ヒストグラムの作成
bulb_df <- tibble(lifetime = bulb)
bulb_df |> 
  ggplot(aes(x = lifetime)) +
  geom_histogram(
    aes(y = ..density..),
    bins = 8,
    fill = "lightblue",
    color = "black",
    alpha = 0.7
    ) +
  stat_function(
    fun = dnorm,
    args = list(mean = averate, sd = standard_dev),
    color = "red",
    size = 1
    ) +
  theme_minimal()

そのとき、この手元にある16個の電球の寿命データの平均値が、母平均1800の周りでどの範囲に入るかを95％の信頼水準で推定してみましょう。

# t検定の実施
bulb_ci <- t.test(bulb)
bulb_ci$conf.int

[1] 1868.890 2033.498
attr(,"conf.level")
[1] 0.95

# 信頼区間の出力（デフォルトで 95％信頼水準）

# コード5-6
qnorm(0.025, lower.tail = FALSE)

[1] 1.959964

# 標準正規分布における上側 2.5％点の分位点を求める。

# コード5-7
n <- length(bulb) # 標本サイズの計算
z <- qnorm(0.025, lower.tail = FALSE) # 上側2.5％点の分位点
xbar <- mean(bulb) # 標本平均の計算
sigma <- 180 # 母標準偏差の設定

# 信頼区間の計算
upper <- xbar + z * (sigma / sqrt(n))
lower <- xbar - z * (sigma / sqrt(n))

# 結果のまとめと出力
ci.bulb <- matrix(c(lower, upper), nrow = 1)
colnames(ci.bulb) <- c("ci.lower", "ci.upper")
knitr::kable(
  ci.bulb, 
  caption = "Bulb data CI (95％)", 
  align = "cc"
  )

Bulb data CI (95％)
ci.lower	ci.upper
1862.995	2039.392

# コード5-8
mean(bulb)

[1] 1951.194

# コード5-9
t.test(bulb, alternative = "two.sided", mu = 1800)


    One Sample t-test

data:  bulb
t = 3.9155, df = 15, p-value = 0.001377
alternative hypothesis: true mean is not equal to 1800
95 percent confidence interval:
 1868.890 2033.498
sample estimates:
mean of x 
 1951.194

# コード5-10
n <- length(bulb) # 標本サイズの計算
z <- qnorm(0.025, lower.tail = FALSE) # 上側2.5％点の分位点
xbar <- mean(bulb) # 標本平均の計算
sigma <- 180 # 母標準偏差の設定
mu <- 1800 # 帰無仮説の平均値 
# Z値の計算
Z <- (xbar - mu) / (sigma / sqrt(n))
Z

[1] 3.359861

# コード5-11
# データの読み込みと前処理
firmdata <- readxl::read_xlsx("data/MktRes_firmdata.xlsx")
firm2018 <- firmdata |>
  filter(fyear == 2018) |>
  mutate(
    ad_dummy = ifelse(adint > median(adint), 1, 0)
    )
# t検定の実施
t.test(sales ~ ad_dummy, data = firm2018)


    Welch Two Sample t-test

data:  sales by ad_dummy
t = -3.3989, df = 85.686, p-value = 0.001029
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
 -1674283.1  -438496.1
sample estimates:
mean in group 0 mean in group 1 
       725009.7       1781399.3

# コード5-12
var.test(sales ~ ad_dummy, data = firm2018, ratio = 1)


    F test to compare two variances

data:  sales by ad_dummy
F = 0.10863, num df = 74, denom df = 71, p-value < 2.2e-16
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.06821636 0.17258882
sample estimates:
ratio of variances 
         0.1086276

# 効果量d=0.8、検出力0.8、有意水準0.05の場合の
# 必要標本サイズの推定
est_n <- pwr.t.test(d = 0.8, power = 0.8, sig.level = 0.05)
est_n


     Two-sample t test power calculation 

              n = 25.52458
              d = 0.8
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

NOTE: n is number in *each* group

# コード5-15
plot(est_n)

# コード5-16
pwr.t.test(d = 0.8, n = 26, sig.level = 0.05)


     Two-sample t test power calculation 

              n = 26
              d = 0.8
      sig.level = 0.05
          power = 0.8074866
    alternative = two.sided

NOTE: n is number in *each* group

## 第5章基礎統計学復習経営学部の1回生がおおよそ学んでいるであろう統計学の内容を、Rコードとともに復習します。テキストの内容に加えて、コードの効率化や関数化についても解説します。 ```{r setup} # コード5-1 pacman::p_load(tidyverse, readxl, pwr, knitr) ``` まずサイコロを作ります。 1から6の数値からなるベクトル`dice`を作ります。つぎに`sample()`関数を用いて、`dice`の中から1つの数をランダムに取り出します。 :::{.aside} `sample()`関数の引数は以下の通りです。 - `x`: 抽出元のベクトル - `size`: 抽出する数 - `replace`: 置換を許可するかどうか ::: ```{r code5-2} # コード5-2 set.seed(442) # サイコロ dice <- 1:6 # サイコロを1回振る d <- sample(dice, size = 1) # 結果の表示 d ``` 同じコードを3回以上書いているなと思ったら、ループ処理や関数化を検討しましょう。 `for()`関数を使ったり、サイコロを投げる回数を引数にした関数を作成したりする方法があります。 ```{r dice-throw} d1 <- sample(dice, size = 10, replace = TRUE) d2 <- sample(dice, size = 10, replace = TRUE) d3 <- sample(dice, size = 10, replace = TRUE) ``` をまとめて書くと、次のようになります。 ```{r} d <- list() # 空のリストを作成 for (i in 1:3) { # ループ処理 d[[i]] <- sample(dice, size = 10, replace = TRUE) } d # 結果の表示 ``` これをさらに発展させて、後で変更する可能性のある変数である試行回数やサイコロを振る回数を別に設定しておいて、コードを読みやすく、また変更しやすくします。 ```{r} #| code-line-numbers: 1:2 trial_count <- 3 # 試行回数 sample_sizes <- 10 # サイコロを振る回数 # 空の数列を作成 d_mean <- numeric(trial_count) for (i in 1:trial_count) { d_mean[i] <- mean(sample(dice, size = sample_sizes, replace = TRUE)) } d_mean ``` 試行回数`trials`とサイコロを振る回数`sample_sizes`を引数にした関数を作成すると、さらに便利になります。ここでは、自作関数`dice_mean()`を定義してみましょう。デフォルト値として、試行回数を3回、サイコロを振る回数を10回に設定しています。 ```{r dice_mean_function} dice_mean <- function(trial_count = 3, sample_sizes = 10) { d_mean <- numeric(trial_count) # 空のリストを作成 for (i in 1:trial_count) { d_mean[i] <- mean(sample(dice, size = sample_sizes, replace = TRUE)) } return(d_mean) } ``` ```{r use_dice_mean} dice_mean() # デフォルト値で実行 dice_mean(trial_count = 5, sample_sizes = 20) # 引数を指定して実行 dice_mean(5, 20) # 省略形でもOK ``` この関数を使って、サイコロを10回、100回、1,000回振ったときの平均をそれぞれ3回ずつ試行してみましょう。 ```{r} set.seed(352) size = c(10, 100, 1000) d_mean <- list() for (i in seq_along(size)) { d_mean[[i]] <- dice_mean(trial_count = 3, sample_sizes = size[i]) } d_mean_mat <- do.call(rbind, d_mean) # 列名を付与 colnames(d_mean_mat) <- paste0("試行", 1:ncol(d_mean_mat)) # 行名を付与 rownames(d_mean_mat) <- paste0("サイコロを", size, "回振る") # 表として出力 knitr::kable( d_mean_mat, caption = "サイコロの標本平均（各サイズ3回ずつ）", align = "ccc" ) ``` `tidyverse`の`purrr`パッケージを使うと、さらに簡潔に書けます。 ```{r} set.seed(352) trial_count <- 3 # 試行回数 sizes <- c(10, 100, 1000) # サイコロを振る回数 # map_dfr()でデータフレームを直接作成 d_mean_df <- map_dfr(sizes, function(n) { # replicate()で3回繰り返し試行 results <- replicate( trial_count, mean(sample(dice, size = n, replace = TRUE)) ) set_names(results, paste0("試行", 1:trial_count)) }) result <- d_mean_df |> mutate(条件 = paste0("サイコロを", sizes, "回振る")) |> relocate(条件) # 「条件」列を一番左へ移動 # 出力 kable( result, caption = "サイコロの標本平均（各サイズ3回ずつ）", align = "cccc" ) ```  ## 区間推定いま、標準正規分布に従う母集団から抽出した無作為標本$Z_1, Z_2, ..., Z_n$を考えます。標準正規分布の確率分布関数$\phi(z)$を用いると、次のように表せます。 $$ P \left( -z_{\frac{\alpha}{2}} \leq \frac{\bar{Z} - \mu}{\frac{\sigma}{\sqrt{n}}} \leq z_{\frac{\alpha}{2}} \right) = 1 - \alpha $$ ### 電球の寿命データで区間推定例として、次のような白熱電球の寿命データが手元にあるものとしましょう。ここで、以下のことが分かっているものとします。 - 平均寿命が1800時間である。 - 寿命の標準偏差$\sigma $は180時間である。 - 平均寿命は正規分布に従う。図にすると次のようになります。 ```{r normal_dist} averate <- 1800 standard_dev <- 180 rnorm(1000, mean = averate, sd = standard_dev) |> as_tibble() |> ggplot(aes(x = value)) + stat_function( fun = dnorm, args = list(mean = averate, sd = standard_dev), color = "red", size = 1 ) + theme_minimal() ``` この分布から抽出したと想定される16個の電球の寿命データが以下の通りです。 ```{r code5-5} # コード5-5 bulb <- c(# 白熱電球の寿命データ 1939.6, 1680.3, 1982.1, 2215.6, 2092.5, 1928.9, 2003.8, 1955.5, 1800.1, 1659.5, 2066.2, 2107.2, 2085.5, 1878.6, 2007.6, 1816.1 ) mean(bulb) sd(bulb) ``` 先ほどの分布にこのデータのヒストグラムを重ねると、次のようになります。 ```{r} # ヒストグラムの作成 bulb_df <- tibble(lifetime = bulb) bulb_df |> ggplot(aes(x = lifetime)) + geom_histogram( aes(y = ..density..), bins = 8, fill = "lightblue", color = "black", alpha = 0.7 ) + stat_function( fun = dnorm, args = list(mean = averate, sd = standard_dev), color = "red", size = 1 ) + theme_minimal() ``` そのとき、この手元にある16個の電球の寿命データの平均値が、母平均1800の周りでどの範囲に入るかを95％の信頼水準で推定してみましょう。 ```{r t-test_bulb} # t検定の実施 bulb_ci <- t.test(bulb) bulb_ci$conf.int # 信頼区間の出力（デフォルトで 95％信頼水準） ``` ```{r code5-6} # コード5-6 qnorm(0.025, lower.tail = FALSE) # 標準正規分布における上側 2.5％点の分位点を求める。 ``` ```{r code5-7} # コード5-7 n <- length(bulb) # 標本サイズの計算 z <- qnorm(0.025, lower.tail = FALSE) # 上側2.5％点の分位点 xbar <- mean(bulb) # 標本平均の計算 sigma <- 180 # 母標準偏差の設定 # 信頼区間の計算 upper <- xbar + z * (sigma / sqrt(n)) lower <- xbar - z * (sigma / sqrt(n)) # 結果のまとめと出力 ci.bulb <- matrix(c(lower, upper), nrow = 1) colnames(ci.bulb) <- c("ci.lower", "ci.upper") knitr::kable( ci.bulb, caption = "Bulb data CI (95％)", align = "cc" ) ``` ```{r code5-8} # コード5-8 mean(bulb) # コード5-9 t.test(bulb, alternative = "two.sided", mu = 1800) # コード5-10 n <- length(bulb) # 標本サイズの計算 z <- qnorm(0.025, lower.tail = FALSE) # 上側2.5％点の分位点 xbar <- mean(bulb) # 標本平均の計算 sigma <- 180 # 母標準偏差の設定 mu <- 1800 # 帰無仮説の平均値 # Z値の計算 Z <- (xbar - mu) / (sigma / sqrt(n)) Z ``` ```{r code5-11} # コード5-11 # データの読み込みと前処理 firmdata <- readxl::read_xlsx("data/MktRes_firmdata.xlsx") firm2018 <- firmdata |> filter(fyear == 2018) |> mutate( ad_dummy = ifelse(adint > median(adint), 1, 0) ) # t検定の実施 t.test(sales ~ ad_dummy, data = firm2018) ``` ```{r code5-12} # コード5-12 var.test(sales ~ ad_dummy, data = firm2018, ratio = 1) ``` ```{r code5-14} # 効果量d=0.8、検出力0.8、有意水準0.05の場合の # 必要標本サイズの推定 est_n <- pwr.t.test(d = 0.8, power = 0.8, sig.level = 0.05) est_n ``` ```{r code5-15} # コード5-15 plot(est_n) ``` ```{r code5-16} # コード5-16 pwr.t.test(d = 0.8, n = 26, sig.level = 0.05) ```