統計ブログはじめました!

各専門分野の統計技術、方法、テクニックなどを気ままに分かり易く例題をもとに解説します。

統計のコツのこつ(41)

2017-04-18 18:03:17 | 日記・エッセイ・コラム
ここしばらくは、「すぐに役立つ統計のコツ」では触れていない効果量(ES)や検出力(Power)の求め方をご紹介して来ました。
ES や Power の筆算でのキーワードは「非心t分布」です。これを筆算するのはやっかいですが、幸いにも"CASIO" の電卓??(下記URL)を使えば楽ちんです。
 
http://keisan.casio.jp/keisan/service.php
 
やってみましょう!
上記URL画面より、
[専門的な計算]→[統計関数]→[非心t分布]→[非心t分布]を選択
 
***
http://keisan.casio.jp/exec/system/1161228870
 
非心t分布の確率密度、下側累積確率、上側累積確率を求めます。
***
 
図1 未入力画面 
 
 
図1の画面で次により計算した値を入力します。
・パーセント点(x):t(α, df)
・自由度(v):n1+n2-2
・非心度(λ):es × √(n1*n2/(n1+n2))
 
例えば、
2つの平均値の差(対応なし)の検定データから、次の値から Power を求めて見ましょう。
n1 = 30             
n2 = 20              
効果量d = 0.954      
有意水準α = 0.05
 
まずは、下記の筆算を行います。
 
***
● 自由度(v)=30+20-2=48
● パーセント点(x)=t(0.05, 48)=2.0106:(注1)
● 非心度(λ)=0.954×√30×20/(30+20)=3.3047
 
注1:
t分布のパーセント点は、
[専門的な計算]→[統計関数]→[スチューデントのt分布(パーセント点)]→
[累積確率:0.025, 自由度:48]→[2.0106]
 
として求めれば良いでしょう。
***
 
そして、
上記の値をそれぞれ図2の様に入力し[計算](6桁に変更)をクリックすると、図2の結果が得られます。
 
図2 入力と計算結果の画面 
 
 
上側累積確率の値(0.8994)が Power(1-β)=89.94% となります。
 
では、
次の「対応あり」のデータをやって見ましょう。
 
n = 101                           
効果量d = 0.322      
有意水準α = 0.05
 
***
●自由度(v)=101-1=100
●パーセント点(x)=t(0.05, 100)=1.984
●非心度(λ)=0.322×√101=3.236
***
 
上側累積確率の値(0.893) が Power(1-β)=89.3% となります。
 
*******************
統計のコツのこつ(40)に対するコメントをご紹介します。
 
杉本典夫 先生のコメント要旨
 
対応の無いt検定(一元配置分散分析)と対応の有るt検定の違いは、前者が2群の平均値の差が0かどうかを検定するのに対して、後者が変化量平均値が0かどうかを検定する点です。そのため「一般的には有意差の出やすい方法を採用」するのではなく、時期変動を検討するための評価指標として、平均値の差が妥当なのか、それとも変化量平均値が妥当なのかで使い分けるべきです。僕のサイト(http://www.snap-tck.com/room04/c01/stat/stat03/stat0303.html)にこのことを端的に表す例を挙げてあるので参考にしてください。
 
情報統計研究所はここから!
 
 
 
 
 
 

この記事をはてなブックマークに追加

統計のコツのこつ(40)

2017-04-11 17:06:26 | 日記・エッセイ・コラム
前回の「統計のコツのこつ(39)」でご紹介した例題は誤解を招く恐れがありますので、言い訳みたいですが、説明させて頂きます。
統計のコツのこつ(36)~(39)までは、一元配置分散分析(ANOVA)における効果量の求め方についてご紹介してきました。
しかし、
例題は同一被験者の薬剤投与前後における値ではないか・・、だとすれば、これは「対応のある検定」を適用するべきだ・・、と言うことです。
すなわち、
各被験者の「投与前と投与後」の"平均値の差"ではなく、"差の平均値"が有意であるかどうかです(一般的には有意差の出やすい方法を採用します)。
ごもっともですが、
この例題では、対応を考慮せず(前後の変化量=前後差を比較せず)に平均値差の検定と効果量の求め方をご紹介しました。
対応する2標本の有意差検定については、
「すぐに役立つ統計のコツ」の「第3章:4.対応する2標本のの有意差検定(21ページ)」や「第7章:二元配置分散分析(100ページ)」、あるいは、「やさしい医学統計手法」(http://kstat.sakura.ne.jp/medical/med_016.htm)を参考にして下さい。
 
この様な誤解を避けるためにも、論文などでは「対応のあるt検定(Paired t-test)」と記載しておくと良いでしょう。
 
前回の「統計のコツのこつ(39)」では、効果量計算サイトで、
「2.Comparison of groups with different sample size(Cohen's d, Hedges' g)」
 
を選択しましたが、
「対応のあるt検定」を用いるときは、
「4.Calculation of d from the test statistics of dependent and independent t-tests」
 
を用いて下さい。
 
その結果は次の通りです。
***
Mode of testing [dependent]
Student t value [6.1052     ] # Paired t-test の値
n1                   [6              ]
r                     [0.472        ] # 相関係数
 
Effect Size d     [2.561       ]
 
なお、
「G*Power」から求めたPower(1-βerror prob)=0.997(99.7%)でした。
***
 
貴重なるご指摘に感謝致します。
 
 

この記事をはてなブックマークに追加

統計のコツのこつ(39)

2017-04-10 13:41:10 | 日記・エッセイ・コラム
前回の続きです。
今回は無料のWebサイト(効果量計算サイト)を利用し、「すぐに役立つ統計のコツ」(31ページ、表4.1)の「投与前と4か月後」を例題にしてみましょう。
データは次のように要約されています。

    投与前    4か月後
  n           6             6
mean   221.3  150.7
 sd          28.68      26.39
 
計算サイトの1例として、下記のURLにアクセスして下さい。
 
https://www.psychometrica.de/effect_size.html#cohenb
 
Top 画面:
 
 
効果量計算の選択(No.2 を選択):
 
 
統計量の入力画面:
 
 
効果量(Effect Size d_cohen)=2.562
効果量の95%信頼限界=1.035~4.089
 
なお、
「G*Power」から求めたPower(1-βerror prob)=0.9784(97.84%)でした。
 
その他にも色々なサイトがありますの、自己責任でご利用されたらいかがでしょうか・・・。
 
その他、Webサイトの1例:
 
http://www.uccs.edu/~lbecker/
http://www.socscistatistics.com/effectsize/Default3.aspx

このブログは、統計手法の一端をご紹介しているに過ぎません。
詳しい統計学的な内容は専門書などで学習して欲しいものです。

情報統計研究所はここから!
 

この記事をはてなブックマークに追加

統計のコツのこつ(38)

2017-04-04 19:02:35 | 日記・エッセイ・コラム
このブログは「すぐに役立つ統計のコツ」をもとに書いています。
統計分析は「習うより慣れろ・・!」ではないかと思っています。難しそうな統計学を一から始めるよりも、まずはやって見ましょう。
最近は、色々な無料の統計ソフトや著名な商用統計ソフトの試用版などが容易に使えますので、本書の内容をなぞってみては如何ですか。
このブログでは、無料の統計ソフトとMS-Excelでの方法を中心にご紹介しています。
それでは本題です。
「すぐに役立つ統計のコツ:表4.9(37ページ)」を見て下さい。
 
Kruskal Wallis テスト(K-W test)において、厳密にはバラツキが等しいかを調べます。参考までにその「R」での方法をご紹介しておきましょう。
前回のデータ(一元配置分散分析)の出力結果(縦長にしたデータ)を用いて試して下さい。
 
「R」での方法
***
fligner.test(VALUE ~ FACTOR, data=dat) # Fligner-Killeen(フリグナー・キリーン)
kruskal.test(VALUE ~ FACTOR, data=dat)
 
出力結果:
 
Fligner-Killeen test of homogeneity of variances
data:  VALUE by FACTOR
Fligner-Killeen:med chi-squared = 0.003643, df = 2, p-value = 0.9982
 
# p-value から、有意と言えませんので「バラツキは等しい」と判断します。
# 参考までに、Kruskal-Wallis test の結果は次の通りです。
# 「すぐに役立つ統計のコツ:表4.9(37ページ)」参照。
 
Kruskal-Wallis rank sum test
data:  VALUE by FACTOR
Kruskal-Wallis chi-squared = 10.312, df = 2, p-value = 0.005765
***
 
さて、
ここで一元配置分散分析(ANOVA)であれ、K-W test であれ、有意であればすべての組み合わせか、あるいは、どれかの組み合わせにおいて有意差があると言えます。それを調べるのが多重比較です。
「すぐに役立つ統計のコツ」(第4章、29ページ)では、ノンパラメトリック法の「ホルム・ボンフェローニ」とパラメトリック法の「チューキのHSD」について説明しています。
 
本書では記載していませんが、多重比較でも効果量の記載が求められる様になってきましたので、例題でご紹介しておきましょう。例題は、
前回と同じ「すぐに役立つ統計のコツ:例題(表4.1、31ページ)」です。
まずは、
「統計のコツのこつ(36)」でご紹介した検出力分析ソフト「G・Power」でやってみましょう。
 
●「GPower」による手順「統計のコツのこつ(36)参照)」
例えば、
本書(表4.1)の「投与前と2か月後」の効果量を求めて見ましょう。
 
test family→t-test→Mean:Wilcoxon-Mann-Whitney(two group)→Post hoc
 [Tail(s):Two、Sample size group 1 = 6、Sample size group 2 = 6]→Determine→

 ◎ n1=n2[Mean group 1 = 221.3, Mean group 2 = 195.3, SD group 1 = 28.68, SD group 2 = 29.67]

→Calculate and transfer to main window→Calculate
 
出力結果:
 Effect size d = 0.891、Power = 0.2738(27.38%)
 
注釈:
効果量(Effect size)はLarge(0.80以上)なのに、検出力(Power)は小さい。すなわち、効果量に見合ったデータ数で無い様です。下記のMann-Whitney U-test からも有意でないことが p-value から分かります。
 
***
 Mann-Whitney U-test  
 U = 11.5 , E(U) = 18 , V(U) = 38.86364  
 z-value = -1.043  
 p-value = 0.297107  
***
 
効果量を表す指標は色々あり、正直に言ってどの指標を用いるべきか迷う時もあります。しかし、大きく分ければ
 ・Cohen's d :平均値の差を標準化した効果量(d famiky:d族)
 ・r :相関関係の強さ(association:r family:r族)
 
の2種類が多用されていると思います。
 
長くなりましたので、「r family」の計算は次回にします。
それから、
この例題は「パラメトリック検定」のものを用いて紹介しています(念のため・・・)。
 
 
情報統計研究所はここから!
 
 
 
 

この記事をはてなブックマークに追加

統計のコツのこつ(37)

2017-03-27 18:36:42 | 日記・エッセイ・コラム
さて、
前回は「G*Power」によって一元配置分散分析(ANOVA)での効果量の求め方をご紹介しました。
今回は、もう少し詳しいANOVAでの出力結果を求めるために「データ解析環境 R」でやって見ましょう。
前回と同じ、
「すぐに役立つ統計のコツ:例題(表4.3)」で分散分析にける効果量など・・などを求めて見ましょう。
まず、
「R」を立ち上げ、パッケージ(rpsychi)をインストールしておいて下さい。
そして、
下記のコマンドを書き実行して下さい。
 
***
library(rpsychi)
value<- c(
 224,235,220,204,265,180,
 186,190,235,181,225,155,
 147,121,175,166,177,118
)
factor=factor(rep(c("A", "B", "C"), c(6, 6, 6)))
dat<- data.frame(VALUE=value, FACTOR=factor)
dat
ind.oneway(VALUE~ FACTOR, data=dat)
 
出力結果:# 縦長にしたデータ
   VALUE FACTOR
1    224      A
2    235      A
3    220      A
4    204      A
5    265      A
6    180      A
7    186      B
8    190      B
9    235      B
10   181      B
11   225      B
12   155      B
13   147      C
14   121      C
15   175      C
16   166      C
17   177      C
18   118      C
 
> ind.oneway(VALUE~ FACTOR, data=dat)
 
$anova.table # ANOVA表
                      SS     df      MS          F
Between (A) 15330     2   7664.9   9.586
Within          11994   15   799.6     
Total            27324   17
 
# omnibus.es returns a omnibus effect size which is a η2, and it's confidence interval
 
$omnibus.es # イータ二乗とその95%信頼限界
etasq    etasq.lower   etasq.upper 
0.561       0.127           0.717
 
# 前回の「η^2(イータの2乗)=SSa/SSt=0.561」と同じです。
# F検定の標本効果量 f <- sqrt( SSa / SSe ) =√(グループ間の変動/グループ内の変動)=√(15330/11994)=1.131
# 回帰分析の標本効果量 η2 = SSa / SSt=グループ間の変動/全体の変動=15330/27324=0.561(56.1%)
   = R2=決定係数
# 確かに、ややこしいので論文などでは注意して下さい。
 
$raw.contrasts
       mean.diff   lower       upper        std
1-2    26.000    -8.798     60.798    16.326 # 投薬前と2カ月後の平均値の差
1-3    70.667    35.869   105.464    16.326 # 投薬前と4カ月後の平均値の差
2-3    44.667      9.869     79.464    16.326 # 2カ月後と4カ月後の平均値の差
 
# std=標準誤差と呼ばれ、「sqrt(MS_within*(1/n1+1/n2))=799.6*(1/6+1/6)=16.326」です。
 
$standardized.contrasts 
         es      lower     upper      std
1-2  0.919   -0.311   2.15      0.577
1-3  2.499    1.268   3.73      0.577
2-3  1.580    0.349   2.81      0.577
 
# コントラスト(Hedges's g)の標準化平均差、母集団標準化平均差の近似信頼区間、および標準誤差
# Hedgesのg = 平均値の差/√MSe=26/sqrt(799.6)=0.919
 
$power
 small   medium     large
 0.059    0.125      0.263
 
# 検出力は上記の基準から「large」と言えます。

情報統計研究所はここから!
 
 
 
 
 

この記事をはてなブックマークに追加