nora_goes_far

ggplot2の（正確にはHmiscの）関数が返す区間を、自力で計算した出力と比較する

R

目的

ggplot2パッケージの関数で、95%信頼区間など、母数の区間推定の実現値を直接可視化できる。この計算が、自力で計算した値と本当に（ほぼ）一致するか調べる。

参考資料

なお本記事は現時点で執筆中。

`mean_cl_normal()`

データが正規分布にi.i.d.に従うと仮定できる場合。以下の2通りのコードはどちらも同じ出力を返すはず。
なお以下のコードではmtcars$wtの母平均を推定する文脈で可視化しているが、このデータに関して上記の仮定が満たされるかどうかは今回は考慮していない（コードを示すことが目的なので）。

library(tidyverse)

# ggplotの関数で95%信頼区間を計算した場合 ---------

g1 = ggplot(data = mtcars,
            mapping = aes(x = factor(am), y = wt)) +
  stat_summary(fun.data = "mean_cl_normal") +
  coord_cartesian(ylim = c(2, 4.5)) +
  labs(x = "am", title = "g1")

g1


# 自力で95%信頼区間を計算した場合 ---------
g2 = mtcars %>% 
  dplyr::group_by(am) %>% 
  dplyr::summarise(mean = mean(wt),
                   upper_cl = mean(wt) +
                     qt(0.975, n() - 1) * (sd(wt) / sqrt(n())),
                   lower_cl = mean(wt) -
                     qt(0.975, n() - 1) * (sd(wt) / sqrt(n()))
  ) %>% 
  ggplot(mapping = aes(x = factor(am), y = mean)) +
  geom_pointrange(mapping = aes(ymax = upper_cl, ymin = lower_cl)) +
  coord_cartesian(ylim = c(2, 4.5)) +
  labs(x = "am", y = "wt", title = "g2")

g2

`mean_cl_boot()`

標本平均を点で、ノンパラメトリック・ブートストラップ法による、ブートストラップ信頼区間をエラーバーで表したグラフ。ブートストラップ法なので、厳密に一致させることはそもそもできないが、ほぼ同じとみてよい？

library(tidyverse)

# ggplotの関数で95%信頼区間を計算した場合 ---------
g3 = ggplot(data = mtcars,
            mapping = aes(x = factor(am), y = wt)) +
  stat_summary(fun.data = "mean_cl_boot") +
  coord_cartesian(ylim = c(2, 4.5)) +
  labs(x = "am", title = "g3")

g3


# 自力で95%信頼区間を計算した場合 ---------
iter = 3000 # ブートストラップ標本数
m = data.frame(
  a_boot_mean = rep(0, iter),
  m_boot_mean = rep(0, iter)
)

set.seed(123)
for(j in 1:2){
  # ここはややこしいが、am == 0とam == 1があるので
  df = mtcars %>% dplyr::filter(am == (j - 1))
  for(i in 1:iter){
    # dfと同じサイズを復元抽出して、wtの平均を計算
    m[i, j] = dplyr::sample_frac(tbl = df, size = 1, replace = TRUE) %>% 
      dplyr::pull(wt) %>% 
      mean()
  }
}

g4 = mtcars %>% 
  dplyr::group_by(am) %>% 
  dplyr::summarise(mean = mean(wt)) %>% 
  dplyr::mutate(upper = c(quantile(m[,1], 0.975), quantile(m[,2], 0.975)),
                lower = c(quantile(m[,1], 0.025), quantile(m[,2], 0.025))
                ) %>% 
  ggplot(mapping = aes(x = factor(am), y = mean)) +
  geom_pointrange(mapping = aes(ymax = upper, ymin = lower)) +
  coord_cartesian(ylim = c(2, 4.5)) +
  labs(x = "am", y = "wt", title = "g4")

g4