統計のコツのこつ(35)

ここしばらくは、各種検定での効果量の求め方をご紹介したいと思います。ここでは、単に、効果量計算の手法のみですので、統計学的な内容については、
下記URL(統計学入門：杉本典夫先生)をご参考になさって下さい。

http://www.snap-tck.com/room04/c01/stat/stat01/stat0106.html

 

それでは、
前回に引き続き、対応のある２群の有意差検定(ノンパラメトリック)での効果量を「R」で求めてみましょう。

ノンパラメトリックでの代表的な方法として、
「ウイルコックスンの符号順位和検定の方法(Wilcoxon signed rank test)

 

があります。

 

「すぐに役立つ統計のコツ」(オーム社刊)の 22ページを見て下さい。

例題は表 3.11(23ページ)を使ってやって見ましょう。
まずは、
情報統計研究所のホームページ(Topページ)にアクセスし、Top ページの「著書の正誤表と例題」→「Excel Samples」→「Excel Samples(1)」からデータをダウンロードして下さい。

Sheet名「表3.11」を開き、項目名を含む「F2：G21」を選択しコピーして下さい。

「R」へのデータの取り込みは前回と同じ要領です。

 

「R」による方法

***
dat<- read.delim("clipboard", head=T)
dat

x<- dat$前
y<- dat$後
n<- length(dat$前)
n
x
y

res<- wilcox.test(x, y, paired=TRUE)
res

p<- res$p.value
z<- qnorm(1-(p/2))
z

es<-  z/sqrt(n)
es

出力結果：

> res

Wilcoxon signed rank test with continuity correction

data:  x and y
V = 182, p-value = 0.0004905
alternative hypothesis: true location shift is not equal to 0

>
> p<- res$p.value
> z<- qnorm(1-(p/2))
> z
[1] 3.485892
>
> es<-  z/sqrt(n)
> es
[1] 0.7997184　# 効果量

***

 

次回は、
多標本の検定として、一元配置分散分析を統計量(例数、平均値、標準偏差)で行う方法をご紹介する予定です。

「すぐに役立つ統計のコツ」(第4章：29ページ)をご覧下さい。

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

 

 

統計のコツのこつ(34)

前号までは、パラメトリック検定での効果量を統計量から求める方法でした。今回は、ノンパラメトリック検定での求め方をご紹介します。
それでは、
「すぐに役立つ統計のコツ」(オーム社刊)の 14ぺージ(ウイルコックスンの順位和検定の方法)を見て下さい。
本書では、
Excelでの方法を紹介していますが、この例題「表 3.2」におけるＧＯＴの検定結果は「有意差なし」ですので、効果量を求めても、有意差に対する効果量としてはあまり意味がありません。そこで、
「表 3.2」のＧＰＴで効果量を求めてみましょう。

すなわち、
ＧＰＴにおける、対応のない(独立)２群の「Wilcoxon rank sum test」による効果量です。

計算は、
データ解析環境「R」で行って見ましょう。

まずは、
情報統計研究所のホームページ(Topページ)にアクセスし、Top ページの「著書の正誤表と例題」→「Excel Samples」→「Excel Samples(1)」からデータをダウンロードして下さい。
そして、
Sheet名「表 3.2」を開き、項目名を含めた「CグループとHグループ」のデータ(D2：E22)を選択しコピーして下さい。

 

「R]での方法：

***

# R Console に下記のコマンドを書き実行する。
dat<- read.delim("clipboard", heade=T)

# データを確認する。
dat

 

# ファイル→新しいスクリプトをクリックしRエディターに次のコマンドを書き実行する。

x<- dat$Cグループ
y<- dat$Hグループ
x　# Cグループ・データの確認
y　# Hグループ・データの確認

 

出力結果：
> x
 [1]  71  56  34  85  29  32  24 123 100  35 153  51  35  61  30 114  41  67
[19]  44  87
> y
 [1]  59 250  53 277  94 146 165  84  53  97 252  48  57 145 107  55  57 121
[19] 319 317

間違いなければ、次のコマンドを実行して下さい。

 

# 検定結果
res<- wilcox.test(x, y, paired=FALSE)
res

 

出力結果：本書(17ページ)参照。
> res
        Wilcoxon rank sum test with continuity correction
data:  x and y
W = 89, p-value = 0.002795
alternative hypothesis: true location shift is not equal to 0

 

そして、効果量は次の通りです。

 

# 効果量(es)
p<- res$p.value
z<- qnorm(1-(p/2))
z
es<-  z/sqrt(N)
es

 

出力結果：
> p<- res$p.value
> z<- qnorm(1-(p/2))
> z
[1] 2.989456

 

＃ 効果量(es)
> es<-  z/sqrt(N)
> es
[1] 0.47267

 

次回は、
Wilcoxon test(対応のある場合)です。対応のある２群の検定は、単に、
　res<- wilcox.test(x, y, paired=TRUE)

 

とするだけですが、一応、本書の22ページを｢R」で行ってみます。

最近、

医学での統計学的検定においても、ｐ値と合わせて効果量の記載が多く見受けられる様になって来ました。
ここしばらくは、効果量の計算についてご紹介したいと思います。
それでは又！

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

 

 

 

 

 

 

統計のコツのこつ(33)

前々回の「統計のコツのこつ(31)」では、標本効果量の筆算手順をExcel関数で示しました。
ここは一つ、
データ解析環境「R」での方法をご紹介しておきましょう。

それでは、
前前号「統計のコツのこつ(31)」の例題を見て下さい。

例題は、

(A)年齢 30～39才, 30名の平均値 Ｘa＝122.5 mmHg, 標準偏差 SDa＝ 10.85 mmHg
(B)年齢 40～49才, 20名の平均値 Ｘb＝133.4 mmHg, 標準偏差 SDb＝ 12.24 mmHg

でした。

「R」のパッケージ「lessR」をインストールし、以下のコマンドを実行して見ましょう。

「R」での方法（その１）

****
library(lessR)
tt.brief(n1=30, m1=122.5, s1=10.85, n2=20, m2=133.4, s2=12.24)

 

出力結果：

Compare Y across X levels Group1 and Group2
--------------------------------------------------------------

Y for X Group1:  n = 30,  mean = 122.50,  sd = 10.85
Y for X Group2:  n = 20,  mean = 133.40,  sd = 12.24

---
t-cutoff: tcut =  2.011
Standard Error of Mean Difference: SE =  3.297

Hypothesis Test of 0 Mean Diff:  t = -3.306,  df = 48,  p-value = 0.002

Margin of Error for 95% Confidence Level:  6.629
95% Confidence Interval for Mean Difference:  -17.529 to -4.271

Sample Mean Difference of Y:  -10.900
Standardized Mean Difference of Y, Cohen's d:  -0.954
***

 

この「lessR」を使えば簡単に「独立２群のｔ検定」と「効果量(Cohen's d)」を求めることが出来ます。

もう一つの方法は、

「R」のパッケージ「rpsychi」をインストールし、以下のコマンドで実行して見ましょう。

 

「R」での方法（その２）

***
library(rpsychi)
ind.t.test.second(m=c(122.5,133.4), sd=c(10.85,12.24), n=c(30,20))

 

出力結果：
> ind.t.test.second(m=c(122.5,133.4), sd=c(10.85,12.24), n=c(30,20))
$samp.stat
    m1    sd1     n1     m2    sd2     n2
122.50  10.85  30.00 133.40  12.24  20.00

$raw.difference
mean.diff     lower     upper       std
  -10.900   -17.529    -4.271     3.297

$standardized.difference
    es  lower  upper    std
-0.939 -1.536 -0.343  0.304
***

ここでの効果量(es)は Hedges'g値，g値の信頼区間とその標準誤差となっています。

よって、
「R」での方法(その1)の効果量と一致しませんのでご注意下さい。

最近では、 Hedges'g を用いる場合も多い様ですので、専門書等を参考にして下さい。

 
参考図書：
　すぐに役立つ統計のコツ(オーム社刊)

 

 

 

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

 

 

統計のコツのこつ(32)

このブログは「すぐに役立つ統計のコツ」(オーム社）に紹介されている内容に沿って書いています。
本書を見ながら、このブログをお読みいただければ分かり易いかと思います。

今回も、例数、平均値、標準偏差を知って、対応のある(ペア)２群のｔ検定の手順をご紹介しましょう。

まずは、
本書(21ページ)の「４. ２つのデータの差を比較する　～対応する２標本の有意差検定～」を開いて下さい。

 

例題は、
情報統計研究所のホームページ(やさしい医学統計手法)から引用します。
次の URLにアクセスして下さい。

URL：
http://kstat.sakura.ne.jp/medical/med_016.htm

 

[例題22] (引用先：やさしい医学統計手法)
高血圧患者 10名に対する降圧剤の投与前と後の血圧値を次にに示します．
投与後の血圧は投与前より有意に低いと云えるかどうかを検定します。

本例題の統計量は次の通りです。
------------------------------------------
N10 　　　　　   ：ペアの例数
diff_Mean 34.8  ：差の平均値
diff_SD  20.094 ：差の標準偏差
diff_SE  6.354   ：差の標準誤差
------------------------------------------
なお、sd＝√V、se＝sd／√N となります。

 

それでは、
この統計量を用い、データ解析環境「R」でｔ検定、効果量、検出力を求めて見ましょう。

 

「R」による方法：

検出力を求めるために、事前に、パーケージ「pwr」をインストールしておきます。そして、
-----------------------------------------------------------
Ndiff_Meandiff_SDdiff_SE

 # t-value
t_valround(t_val,3)

 # p-value
round(2*pt(-abs(t_val), N-1),5)

 # 95% CI
tLowerUpperround(Lower,3); round(Upper,3)

 # Effect size
esround(es,3)

 

 # Power
library(pwr)
power.t.test(n=N, delta=es, sig.level=0.05, type='paired',
  strict=T)$power

 

***

出力結果

>  # t-value
> t_val> round(t_val,3)
[1] 5.477
>
>  # p-value
> round(2*pt(-abs(t_val), N-1),5)
[1] 0.00039
>
>  # 95% CI
> t> Lower> Upper>round( Lower,3); round(Upper,3)
[1] 20.426
[1] 49.174
>
>  # Effect size　(効果量)
> es> round(es,3)
[1] 1.732
>
>  # Power　(検出力)
[1] 0.9979495　→(99.8%)
-----------------------------------------------------------

 

次回も、例数、平均値、標準偏差から色々な検定の計算をご紹介する予定です。

 

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

 

 

統計のコツのこつ(31)

 

このブログは「すぐに役立つ統計のコツ」(オーム社）に紹介されている内容に沿って書いています。
本書をお読みいただければ分かり易いかと思います。

今回は前回と同じように、
２つの標本の大きさ(n1,n2)、平均値(m1,m2)、標準偏差(s1,s2)だけが分かっている場合の"独立２群の平均値の差"の検定をやって見ましょう。

● 本書の第３章（10ページ）を参考にして下さい。

例題は、
「やさしい医学統計手法」(情報統計研究所：下記URL)から引用します。

URL:

http://kstat.sakura.ne.jp/medical/med_014.htm

 

例題は次の様になっています。

(A)年齢 30～39才, 30名の平均値 Ｘa＝122.5 mmHg, 標準偏差 SDa＝ 10.85 mmHg
(B)年齢 40～49才, 20名の平均値 Ｘb＝133.4 mmHg, 標準偏差 SDb＝ 12.24 mmHg
(C)年齢 50才～～, 10名の平均値 Ｘc＝139.0 mmHg, 標準偏差 SDc＝ 20.40 mmHg

 

筆算の過程は「やさしい医学統計手法」を見て下さい。

ここでは、
Excelの関数を使った方法のご紹介です。

図１の様なフォーマットを作ります。

 

図１　ｔ検定のためのフォーム

 

 

そして、以下の関数式を各セルに入力しておいて下さい。

 

・平均値の差 Ｃ列５行(1)：=C2-C3
・95%CI_Lower Ｄ列５行(2):=ROUND(C5-T.INV.2T(0.05,B2+B3-2)*E10,3)
・95%CI_Upper Ｅ列５行(3)：=ROUND(C5+T.INV.2T(0.05,B2+B3-2)*E10,3)
・自由度　
　Ａ群 Ｃ列７行(4)：=B2-1
　Ｂ群 Ｄ列７行(5)：=B3-1
  全体 Ｅ列７行(6)：=B2+B3-2
・併合した分散
　Ａ群 Ｃ列８行(7)：=D2^2*(B2-1)
　Ｂ群 Ｄ列８行(8)：=D3^2*(B3-1)
　全体 Ｅ列８行(9)：=(C8+D8)/E7
・併合した標準偏差 Ｅ列９行(10)：=SQRT(E8)
・平均値差の標準誤差 Ｅ列10行(11)：=E9*SQRT(1/B2+1/B3)
・ｔ値 Ｅ列11行(12)：=C5/E10
・ｐ値 Ｅ列12行(13)：=T.DIST.2T(ABS(E11),E7)
・標本効果量 Ｃ列13行(14)：=ROUND(ABS(E11*SQRT((B2+B3)/(B2*B3))),5)
　S_e.s. Ｅ列13行(15)：=SQRT((B2+B3)/(B2*B3)+C13^2/(2*E7))
・e.s. 95%CI Ｃ列14行(16)：=ROUND(C13-1.96*E13,3)
　e.s. 95%CI Ｄ列14行(17)：=ROUND(C13+1.96*E13,3)

 

ここで、
例数(n)、平均値(m)、標準偏差(sd)を入力すれば、次の結果(Student's t-test)が出力されます。

 

・平均値の差の95% 信頼区間
・ｔ値
・ｐ値
・効果量(e.s.)
・効果量の95％信頼区間

 

それでは、黄色のセルに例題(A)の数値を入力して下さい。

 

検定結果は、
図２で示した緑色セルに次のように出力されます。

 

図２　独立２群(対応のない)ｔ検定の結果

 

 

・平均値の差＝-10.9（95% CI＝-17.529～-4.271）
・ｔ値＝-3.3062（ｐ値＝0.0018）
・効果量(e.s.)＝0.954（e.s.の95％CI＝0.357～1.552）

なお、
分散比、及び、その検定と Welch method による t 値、自由度、p 値 は 図３の関数式で求められます。

 

図３ Welch method による関数式

　

 

この例題での分散比は 1.273(ｐ＝0.545)で等分散と見なされます。

もし、
等分散でなければ「Welch method」の結果の採用を検討します。

詳しくは、
「すぐに役立つ統計のコツ」(18ページ)を見て下さい。

 

また、
検出力は図４、図５の関数式で近似値を求めることが出来ます。

 

図４　Cohen による近似値を求める関数式

 

 

図５　Mathews による近似値を求める関数式

 

 

・検出力（1-β）：
    Cohen近似値＝0.897926(89.8%)、Mathews近似値＝0.89963(89.9%)

 

・「R」による正確な近似値＝0.8993719（89.9％）

 

「R」による正確な検出力は次により実行して下さい。
---------------------------------------------------
library(pwr)
pwr.t2n.test(n1=30, n2=20, d=0.95, sig.level=0.05)

出力結果  
t test power calculation            
n1 = 30             
n2 = 20              
d = 0.954      
sig.level = 0.05          
power = 0.8993719    
alternative = two.sided
---------------------------------------------------

最近の論文等の記述では、効果量の併記が多く見られます。

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