UnsplashのBob Brewerが撮影した写真
前回は、cor.test()関数を使って病院の数の変化幅と、開始年の病院の数の相関を検定しました。結果は、相関係数は-0.725と負の強い相関があり、95%信頼区間は-0.838 ~ -0.553でした。
今回は、inferパッケージを使って、この病院の数の変化幅と開始年の病院の数の相関を推論してみます。
まず、specify()関数とcalculate()関数で相関係数を計算します。
-0.725と(当たり前ですが)、cor.test()関数の結果と同じになりました。
specify(chg ~ y1999)としているので、これから、chgがy1999と関係しているかどうかを調べますと設定しています。そして、calculate(stat = "correlation")としているのでchgとy1999の相関を計算しています。
次に、specify()関数とhypothesize()関数とgenerate()関数とcalculate()関数で、もしもchgとy1999が相関が無ければ、と仮定した世界での相関係数を何回も計算します。
hypothesize(null = "independence")としているのが、chgとy1999が独立しています、相関は無いです、という帰無仮説の設定の部分です。
そして、generate(reps = 1000, type = "permute")で1000回、chgとy1999をそれぞれ独立にシャッフルした組み合わせを1000セット作っています。
calculate(stat = "correlation")のところで1000セットの相関係数を計算しています。
replicateの値が1000セットの番号で、statがそれぞれの相関係数です。
visualize()関数でこの1000個の相関のヒストグラムを描き、shade_p_value()関数でp値の部分を重ねます。
赤い垂線が、実際の相関係数の-0.725の水準です。こうしてみると、-0.725という相関係数は、かなり起こりにくい(実際、1000回の試行ではおきませんでした)値だとわかります。
get_p_value()関数でp値を計算します。
p値は0です。cor.test()関数では、p値は8.177e-09と表示されていました。これは、0.000000008177ということなので、ほとんどゼロというころですね。
こんどは信頼区間を算出します。はじめに、ブートストラップ法で相関係数を1000回計算します。
statの列の値を見ると、-0.628, -0.740, -0.652などと実際の相関係数の-0.725に近い値なっています。
visualize()関数でヒストグラムを描いて視覚化してみましょう。
右側のほうの裾野が広いヒストグラムですね。
信頼区間を計算しいます。get_confidence_interval()関数を使います。はじめはパーセンタイル法の信頼区間を計算します。
-0.863 ~ -0.465です。
パーセンタイル法の信頼区間は、ブートストラップ法で作成した(今回は1000個の)相関係数の下位2.5%の値と上位2.5%の値を信頼区間の下限と上限にしています。
もうひとつ、スタンダードエラー法があります。これは、ブートストラップ法で作成した相関係数のバラツキ、スタンダードエラーを計算式で求めて、これを利用する方法です。ほんとうは、この方法は、ブートストラップ法で作成した相関係数のヒストグラムが左右対称の正規分布になっていないと使うのは適切ではいそうですが、とりあえずやってみます。
-0.929 ~ -0.521とブートストラップ法の信頼区間よりも左よりですね。
shade_confidence_interval()関数で２つの信頼区間をさきほどのヒストグラムに重ねます。
ブートストラップ法の信頼区間のほうが適切な感じですね。
最後にcor.test()関数の結果と今回のinferパッケージを使ったp値、信頼区間の結果を比較してみます。
相関係数、p値は同じ値、信頼区間は、３つそれぞれ微妙に違います。
今回は、
https://infer.netlify.app/articles/observed_stat_examples
を参考にしました。
今回は以上です。
次回は、
です。
初めから読むには、
です。
