Rで標準誤差を体感する

2023年10月26日2023年12月1日

今回は、標準誤差の理解を深めることが目的です。

私たちが何かを調べるとき、ほとんどの場合には全数調査を行うことは不可能です。そこで標本調査を行うことになるわけですが、そうすると標本から得られた推定値は本来の値とはズレます。そのズレがどのくらいかを見積もるものが標準誤差です。

ここでは一番簡単な例として、母平均を標本平均から推定することを考えます。

この式の分母にサンプル数$N$の平方根が入っています。つまりサンプル数が多くなれば標準誤差は小さくなる関係があります。これは、サンプルを増やせばより精確になるだろう、という直観的な感覚にも合致するのではないでしょうか。

さて、今回はこのことを確かめてみます。

ソース・ファイルは以下からダウンロードできます（zip形式です）

まずは繰り返し使うので、関数として定義します。

#標準誤差の実験

f <- function(N) {
  data <- rnorm(100*N,50, 10)
  mat <- array(data, c(N,100))
  means <- apply(mat, 2, mean)
  cat("標本数=",N,"\n")
  cat("標本平均の平均=",mean(means),"\n")
  cat("標本平均の標準偏差=",sd(means)*sqrt((N-1)/N),"\n")
  cat("標準誤差=",10/sqrt(N),"\n\n")
  hist(means,breaks=seq(0.5,99.5,1))
}

Nはサンプル数です。そしてサンプリングを100回行います。

ここでは一度に100✕サンプル数分の1次元配列を乱数で作ってから、それをarray関数でN✕100の表に変形しています。母集団は平均値50、標準偏差10としています。
図で説明すると、N=10を例にとると、最初は1000個の一列に並んだ乱数列を作ります。

array関数はこれを次元の違う配列(array)に変形します。

そして、apply関数で行ごとに（2つ目の引数の2がその指定）mean関数を適用します。

全体としては、「N個のサンプリングをして平均値を取る」作業を100回繰り返したことになります。

これで100個分の平均値が得られました。これがどのくらいばらつくのか、ということを調べます。

次に、catで表示を行います。関数sd()は不偏分散を返しますので標本分散に戻す式を掛けています。
最後にヒストグラムを作成していますが、これはいろいろと区切り方を変えて試してみたのですが、整数を挟むプラスマイナス0.5にしてちょうど母平均の50を含む形にしました。

N=10
h1<-f(N)

N=100
h2<-f(N)

N=1000
h3<-f(N)

これが実行部分で、サンプル数を10、100、1000と変えています。

ここまでの実行結果を示します。まず、コンソールへの出力ですが、

標本数= 10 
標本平均の平均= 50.3196 
標本平均の標準偏差= 2.964488 
標準誤差= 3.162278 

標本数= 100 
標本平均の平均= 49.82651 
標本平均の標準偏差= 1.020339 
標準誤差= 1 

標本数= 1000 
標本平均の平均= 49.95057 
標本平均の標準偏差= 0.3228284 
標準誤差= 0.3162278 

結果はランダムで変わりますので実際に動かしたものは異なります。

ここでは、標本平均の標準偏差と標準誤差がほぼ同じくらいの数字になっていることが確認できます。
つまり、標本平均を何度も観測したとき、その値のブレ具合を表していることが分かると思います。

これをヒストグラムで確認してみましょう。

左からN＝10, 100, 1000のときのヒストグラムになっています。Nが大きいほど中心に集まっている様子が分かります。

これでは高さが比較しにくいので重ねてみましょう。

plot(0, 0, type="n", xlim=c(40,60), ylim=c(0, 100), xlab="", ylab="")
lines(h1$mids, h1$counts, col="blue")
lines(h2$mids, h2$counts, col="red")
lines(h3$mids, h3$counts, col="cyan")

先程のヒストグラムの結果から区間の代表値midsと度数countsを取り出してグラフにしています。Nが少ない方から、青・赤・シアン（水色）となっています。

こうしてみると、Nが少ない方が横に広がっている様子がよく分かると思います。