Rで標本分散から母分散の信頼区間を求めたい
東京大学出版会の統計学入門（基礎統計学Ⅰ）を眺めていてp.226-227、242の母分散の信頼区間を求めたり、母分散のχ二乗検定をRでやりたいと思った。
最初に自作関数（パッケージ）を使うパターンを示すが、よくわからんパッケージをインストールするのも微妙だろうて、その場合はgithubのコードを追っかけてもらうか、その次に示すステップバイステップで計算する方法を実行するとよいと思う。
Githubにあげている自作関数（パッケージ）をつかう
計算は単純だが、この計算を行う組み込み関数はない？様子なので適当に関数を作ってgithubの{infun}パッケージにvar_()として入れてみた。ついでに入力された値が母集団である場合の母分散も計算するようにした（これも組み込み関数はなさそう）。
計算は教科書通り計算しているつもり。
もしかしたらパッケージを探せばあるかもしれない。
install.packages("remotes") remotes::install_github("indenkun/infun")
でインストールできるようにしている。
例としてirisのSepal.Lengthの母分散の信頼区間を求める。
infun::var_(iris$Sepal.Length)
##
## Chi-squared test
##
## data: iris$Sepal.Length
## X-squared = Inf, df = 149, p-value < 2.2e-16
## alternative hypothesis: true population variance is not equal to 0
## 95 percent confidence interval:
## 0.5531973 0.8725029
## sample estimates:
## sample_variance population_variance
## 0.6856935 0.6811222
χ二乗検定を行って、任意の分散（デフォルトでは0）との検定を行いつつ、信頼区間（デフォルトで95％）を計算する。
sample
estimeteのところに標本分散（普通にvar()したときの答え）に加え、入力された値が母集団とした場合の母分散を表示されるようにしている。
返す答えは"htest"クラスのリストなので、他の検定系の関数と同じように必要な値だけを簡単に取り出すことができる。
# 標本分散（不偏分散）を取り出す infun::var_(iris$Sepal.Length)$estimate[["sample_variance"]]
## [1] 0.6856935
# 入力された値を母集団とした母分散を取り出す infun::var_(iris$Sepal.Length)$estimate[["population_variance"]]
## [1] 0.6811222
# 標本分散から計算した母分散の信頼区間の上限下限を取り出す infun::var_(iris$Sepal.Length)$conf.int[[1]]
## [1] 0.5531973
infun::var_(iris$Sepal.Length)$conf.int[[2]]
## [1] 0.8725029
母分散が0を帰無仮説とする検定をχ二乗検定を行っているが、これは基本的に意味がない。
分散が0ということは、すべての値が同じということなので、一つでも異なる値があれば分散は0にならないので0ということはありえなし、χ二乗値を計算するときに0で割る計算をしているので、標本分散が0でない場合はRが計算でInfを返しているのでそのままつかって計算しているだけなので意味はない。
意味があるのはそう滅多に無いが母分散が既知の場合にその母分散と有意な差があるかどうか（東大出版会の統計学入門の例でいうと例年の分散と今年度の分散がことなるか検定したい場合など）に使える。
たとえば、あやめのデータでがくの長さの分散が0.6で有ることが既知である場合には次のようにすることで、今回のデータがそれと分散が異なるかどうかを検定することができる。
infun::var_(iris$Sepal.Length, pv = 0.6)
##
## Chi-squared test
##
## data: iris$Sepal.Length
## X-squared = 170.28, df = 149, p-value = 0.2237
## alternative hypothesis: true population variance is not equal to 0.6
## 95 percent confidence interval:
## 0.5531973 0.8725029
## sample estimates:
## sample_variance population_variance
## 0.6856935 0.6811222
ちなみに、標本分散が0となるデータの場合に、母分散が0を帰無仮説とする検定を行おうとする場合には、χ二乗値の計算で0÷0がでてくるので、p値がNAになるようになっている。標本分散が0で母分散0と有意な差があるか検定することに意味があるのかどうかはわからない。エラーを返さないようにするための措置。
infun::var_(rep(1, 10))
##
## Chi-squared test
##
## data: rep(1, 10)
## X-squared = NaN, df = 9, p-value = NA
## alternative hypothesis: true population variance is not equal to 0
## 95 percent confidence interval:
## 0 0
## sample estimates:
## sample_variance population_variance
## 0 0
あと、他のt.test()などの検定に倣って片側検定もできるようにしている。
infun::var_(iris$Sepal.Length, pv = 0.5, alternative = "greater")
##
## Chi-squared test
##
## data: iris$Sepal.Length
## X-squared = 204.34, df = 149, p-value = 0.001779
## alternative hypothesis: true population variance is greater than 0.5
## 95 percent confidence interval:
## 0.5724186 Inf
## sample estimates:
## sample_variance population_variance
## 0.6856935 0.6811222
infun::var_(iris$Sepal.Length, pv = 0.9, alternative = "less")
##
## Chi-squared test
##
## data: iris$Sepal.Length
## X-squared = 113.52, df = 149, p-value = 0.01368
## alternative hypothesis: true population variance is less than 0.9
## 95 percent confidence interval:
## 0.0000000 0.8389097
## sample estimates:
## sample_variance population_variance
## 0.6856935 0.6811222
このとき、信頼区間は"less"のとき0～となるようにしているが、これはt検定で"less"で検定したときに信頼区間が必ず-Infとなるようなもので意味はない。先に述べたように、標本分散が0でない時点で母分散が0ということはありえないが、分散は負の値を取り得ないので-Infとするのは明らかにおかしいので0としているだけ。
ステップバイステップで計算する
まずは教科書の例にならって、例年平均50点で分散36だったテストで、今年のランダムに選ばれた25人に実施したら平均53点で分散48だった場合に、この児童の揃い方が例年と異なるかどうかを確認したい場合を考える。
異なるかどうかなので両側の分散のχ二乗検定（帰無仮説：例年の分散（母分散）と今年の母分散が等しい）を行う。教科書では有意水準10％（α = 0.10）なのでそれに倣って行っていく。
とりあえず使いやすいように、使う値をすべて変数に代入しておく。平均点が例示されているが、今回は使わない。
# 標本数 n <- 25 # 母分散（例年の分散） pv <- 36 # 標本分散（今年の分散） ss <- 48 # 有意水準10％ alpha <- 0.10
ここで統計量χ二乗値（χ2）は標本数n、標本分散s2、母分散σ2としたとき
[tex: χ2 = (n - 1)s2/σ2]
となる（うまくTex記法が評価されない（わからん）。めんどくさいからそのままベタ打ちしておく。）ので、
# χ二乗値(統計量）を計算する statistic <- (n - 1) * ss /ps statistic
## [1] 32
とχ二乗値は32となる。
これが自由度（n-1）が24のχ二乗値の90％信頼区間に入っているかどうかを確認する。
両側なので確認したい確率分布のパーセント点は上側が
[tex: χ2_{1-α/2}(n - 1) ]
下側が
[tex: χ2_{α/2}(n - 1) ]
となる（棄却域が上下合わせて10％となるようにα/2となっている）のでそれをqchisq()計算する。
# 上側 chi_upper <- qchisq(1 - alpha / 2, n - 1) chi_upper
## [1] 36.41503
# 下側 chi_lower <- qchisq(alpha / 2, n - 1) chi_lower
## [1] 13.84843
となるので、χ二乗値の90％信頼区間は[13.8 - 36.4]と区間推定でき、32はこの範囲内なので、10％有意水準で帰無仮説は棄却されない、となる。ということで有意な差があるとはいえないとなる。
ここで「今年」の母分散の90％信頼区間を今年の標本分散から区間推定したいと思う。
母分散の信頼区間は標本分散から
[ tex : \frac{(n-1)s2}{χ2_{α/2}(n - 1)}, \frac{(n-1)s2}{χ2_{1-α/2}(n - 1)}\ ]
となるので、これを計算する。
先程求めた値を使って10％信頼区間を求めていく。
# 下限 pv_lower <- (n - 1) * ss / chi_upper pv_lower
## [1] 31.63529
# 上限 pv_upper <- (n - 1) * ss / chi_lower pv_upper
## [1] 83.18635
標本分散から計算される、母分散の10％信頼区間は[31.6 - 83.2]となり、この点からも例年の分散と有意な差がみとめられないということがわかる。
当然、これはベクトルデータの場合でも同様にデータから標本分散を計算することで母分散の信頼区間を計算できる。
例としてirisのSepal.Lengthの母分散の95％信頼区間を求めてみる。
# 標本分散を求める ss <- var(iris$Sepal.Length) ss
## [1] 0.6856935
# 標本数を確認する n <- length(iris$Sepal.Length) n
## [1] 150
# 今回は95％信頼区間（5％有意水準）で計算する alpha <- 0.05
となる。
先程と同様にχ二乗値の上限下限を計算する。
# 上側 chi_upper <- qchisq(1 - alpha / 2, n - 1) chi_upper
## [1] 184.687
# 下側 chi_lower <- qchisq(alpha / 2, n - 1) chi_lower
## [1] 117.098
χ二乗値の95%信頼区間が[117.1 - 184.7]とわかる。
あとはこの値を用いて母分散の95％信頼区間をもとめる。
# 下限 pv_lower <- (n - 1) * ss / chi_upper pv_lower
## [1] 0.5531973
# 上限 pv_upper <- (n - 1) * ss / chi_lower pv_upper
## [1] 0.8725029
ということで母分散の95％信頼区間は[0.55 - 0.87]と区間推定できる。
入力されたデータが母集団とした母分散の計算
Rでvar()で計算される分散は不偏分散なので「各観測値と平均との差を2乗したものの総和」を「サンプル数-1で割る」ので、入力した値が母集団すべての値の場合は母分散が「各観測値と平均との差を2乗したものの総和」を「サンプル数で割る」ので、計算された不偏分散にサンプル数-1をかけ、サンプル数で割ると、母分散がでる。
つまり、
var(iris$Sepal.Length) * (length(iris$Sepal.Length) - 1) / length(iris$Sepal.Length)
## [1] 0.6811222
となる。
