ねこすたっと

ねこの気持ちと統計について悩む筆者の備忘録的ページ。

幾何平均（geometric mean）の信頼区間を求める（Gmean関数）[R]

maths/theory bootstrap R neco

はじめに

平均と言えば大体は「算術平均（arithmetic mean, AM）」のこと。全部足して、要素の個数Nで割る。

\begin{aligned} AM &= \frac{x_1 + x_2 + ... + x_N}{N} \end{aligned}

これに対して、「幾何平均（geometric mean, GM）」は全部掛け合わせて、N乗根をとる（1/N乗する）もの。

\begin{aligned} GM &= \sqrt[N]{x_1 x_2 ... x_N} \end{aligned}

幾何平均の信頼区間

両辺の対数を取ると、

\begin{aligned} \log(GM) &= \frac{\log(x_1) + \log(x_2) + ... + \log(x_N)}{N} \end{aligned}

となるので、算術平均の信頼区間を正規分布をもとに求めるのと同じ要領で、

\begin{aligned} Limits &= \log(GM) \pm t_{1-\alpha/2, df=N-1}SE \\ SE &= \frac{SD(\log(x))}{\sqrt(N)} \end{aligned}

で計算された両端の値を、指数変換して求められる。

Rで95%信頼区間を手計算してみる。

set.seed(1234)
N <- 5
x <- runif(n=N, min=0, max=1)

GM <- exp(mean(log(x)))
log_low <- mean(log(x)) - qt(p=0.975, df=N-1)*sd(log(x))/sqrt(N)
log_upp <- mean(log(x)) + qt(p=0.975, df=N-1)*sd(log(x))/sqrt(N)

ここでqt(p, df)は自由度dfのt分布において、累積確率がpになるt値。Nが大きくなるとdfも大きくなって、標準正規分布に近づく（つまりqt→1.96）。

> GM
[1] 0.4708253

> exp(log_low)
[1] 0.1728172

> exp(log_upp)
[1] 1.282723

Gmean( )を使って計算する

DescToolsパッケージ*1のGmean( )を使ってみる。平均を求めるxはさっきとおんなじ。指定する引数は次のとおり。

  • method："classic"とすると前述の手計算と同じ方法。"boot"とするとブートストラップ（多分）。
  • conf.level：何パーセントの信頼区間を求めるか
  • sides：両側（"two.sided"）か片側（"left"または"right"）を指定
  • na.rm：欠測を除去するかどうか
library(DescTools)
Gmean(x, 
      method = "classic", 
      conf.level = 0.95,
      sides = c("two.sided"), 
      na.rm = TRUE)

     mean    lwr.ci    upr.ci 
0.4708253 0.1728172 1.2827226

手計算と同じ境界値になっている。

ちなみにブートストラップ法を選ぶと少し変わる。

Gmean(x, 
      method = "boot", 
      conf.level = 0.95,
      sides = c("two.sided"), 
      na.rm = TRUE)

     mean    lwr.ci    upr.ci 
0.4708253 0.2956773 0.9874273

おわりに

  • ワクチンの抗体価の記述でgeometric mean titer（GMT）を求められたので勉強してみました。
  • 猫がベッドに上がってくる季節になりました。

ランキング参加中
統計

参考資料

  • どれくらいサンプルがあればt分布を標準正規分布とみなしてよいかについてまとめてくださってます。

blog.apar.jp

*1:Tools for Descriptive Statisticsより