理学療法士がまとめたノート

統計学備忘録(R言語のメモ)

since2016 ときどきTEXのメモ

有意差とは・・・?

乱数を発生させて、set.seed( )で記憶させてシミュレーションしてみます
乱数なので再現できませんが、set.seed( )を使用することで再確認できます

set.seed(1)     #もう一度確かめたいときはset.seed( )で乱数を記憶させておきます.( )の中は何でもOK.
x20 <- rnorm(20,3.25,2.25)  #平均3.25,標準偏差2.25のデータ(乱数)20個
set.seed(2)
y20 <- rnorm(20,2.12,1.95)  #平均2.12,標準偏差1.95のデータ(乱数)20個
t.test(x20,y20)  #独立した2群の差の検定を実施してみます

        Welch Two Sample t-test

data:  x20 and y20
t = 1.808, df = 37.999, p-value = 0.07852
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.1409049  2.4959643
sample estimates:
mean of x mean of y 
 3.678679  2.501149 

♯ p-value =  0.07852 危険率0.05でも有意差は認められませんでした・・・

次に同じ平均値、同じ標準偏差で100個ずつ用意して、検定してみます

set.seed(3)
x100 <- rnorm(100,3.25,2.25)
set.seed(4)
y100 <- rnorm(100,2.12,1.95)
t.test(x100,y100)

        Welch Two Sample t-test

data:  x100 and y100
t = 3.6836, df = 196.81, p-value = 0.0002971
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.4491119 1.4841008
sample estimates:
mean of x mean of y 
 3.274830  2.308224 

P = 0.0002971  危険率0.01でも有意差あり!

統計量(この場合はt値)を考えれば当然の結果なのですが、
まだこのことに気づかずに有意差のみで議論している理学療法士も多いようです.
独立した2群であれば、それぞれのサンプルを増やすことで有意差を出すことが可能になります!

では、箱ひげ図で確認してみます

boxplot(x20,y20,x100,y100,xaxt="n")   
name<-c("x20","y20","x100","y100")
axis(side=1,at=c(1,2,3,4),labels=name) 

f:id:yoshida931:20181217225224p:plain

x20とy20に有意差がなくて、x100とy100に有意差がある・・・ようには見えませんね!
平均値と標準偏差が等しい正規分布からの乱数ですので、同じような箱ひげ図になるのは当然です.
x20とy20に有意差がなくて、x100とy100に有意差があるという結果の理解に苦しみます.
ここに数字のマジックが隠れています.

どのように有意差検定を実施しているのか、理解することが必要になります.
以下のページで、n増加→統計量大→p値小の構図が理解できると思います.

yoshida931.hatenablog.com

さて次はどうでしょうか?

# x20,y20の平均値は同じで、標準偏差を1/10にしてみましょう
set.seed(5)     
x202 <- rnorm(20,3.25,0.225)  
set.seed(6)
y202 <- rnorm(20,2.12,0.195) 
t.test(x202,y202)

    Welch Two Sample t-test
data:  x202 and y202
t = 15.657, df = 37.994, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.8960426 1.1621642
sample estimates:
mean of x mean of y 
 3.186874  2.157771 

# p-value = 2.2e-16  危険率0.01でも有意差あり!

箱ひげ図で確認してみましょう.
f:id:yoshida931:20180124195732p:plain:w400

平均値が同じでも、ばらつき(分散、標準偏差)に違いがあれば統計結果は異なります.

有意差に一喜一憂するのではなく、サンプル数や分散、またそのデータがもつ性質をよく考えて、差を考察するべきです.P値よりも信頼区間や効果量などが比較の参考になるでしょう.
yoshida931.hatenablog.com

平均

算術平均 arithmetin mean (=相加平均)

1回 102回 153回 8人
平均回数は
(10*1+15*2+3*8)/(10+15+8)

度数分布からの平均
真の平均の近似値なので多少のズレが生じます
最初と最後の階級が少ない場合には無視して求めます(年収平均など)
無視できない数であれば、最初の下限は0にして、最後の上限は直線の範囲と同じ範囲にします


幾何平均 geometoric mean (=相乗平均)

そもそも「幾何」とは、「いくばく」と読むので、わずか・すこしという意味があります 例1)利回りの例で考えてみます

10%, 15%, -5%, 13% と年々変化する利回り
平均は
 ( 1.1 + 1.15 + 0.95 + 1.13 ) / 4 
=1.0825 ではありません!
正解は・・・
( 1.1 * 1.15 * 0.95 * 1.13 )^(1/4)
= sqrt( sqrt( 1.1 * 1.15 * 0.95 * 1.13 ) )
=1.079501 
年平均利率は7.95%が正解です

奇数の場合は対数で計算します

例2)極端な例ですが、消費税で考えます

100円の品物の買う場合に
昨年5%になり、翌年15%に引きあがりました
この2年間の平均上昇率は何%でしょうか?
( 5 + 15 ) / 2 = 10% ではありません
(1.05 * 1.15 ) ^ (1/2) =1.098863
平均9.88%が正解です


調和平均 harmonic mean

距離÷時間=速度を例にとって考えてみます

片道100kmを往復します
往路は10km/h, 復路は15km/hかかりました
平均速度は(10+15)/2=12.5km/hではありません
往路の時間は100/10時間、復路の時間は100/15時間で合計50/3時間になります
したがって平均時速は
(100*2)/( (100/10)+(100/15) ) = 12 km/h

データフレームからの抽出

例)ChickWeightからDietの1と3だけ抜き出す

subset(ChickWeight,subset = Diet==c(1,3))

例)iris アヤメ

head(iris,5)
   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1           5.1         3.5          1.4         0.2  setosa
2           4.9         3.0          1.4         0.2  setosa
3           4.7         3.2          1.3         0.2  setosa
4           4.6         3.1          1.5         0.2  setosa
5           5.0         3.6          1.4         0.2  setosa

Sepal.Length(へたの長さ)とSpecies(種)
2列のみを抽出したフレームにします

irisSL <- subset(iris,select = c(Sepal.Length, Species))
head(irisPW,5)

  Sepal.Length Species
1          5.1  setosa
2          4.9  setosa
3          4.7  setosa
4          4.6  setosa
5          5.0  setosa

つぎに列を種(3種類)にしてへたの長さを提示します

iun2 <- unstack(irisSL,formula = Petal.Width ~ Species)
head(iun2,10)

  setosa versicolor virginica
1    5.1        7.0       6.3
2    4.9        6.4       5.8
3    4.7        6.9       7.1
4    4.6        5.5       6.3
5    5.0        6.5       6.5

Sepal.Length(へたの長さ)とSpecies(種)の2列に戻します

x<-c(iun2$setosa,iun2$versicolor,iun2$virginica)
Species <- c(rep("setosa",50),rep("versicolor",50),rep("virginica",50))
xy <- data.frame(Sepal.Length=x, Species)

  Sepal.Length Species
1          5.1  setosa
2          4.9  setosa
3          4.7  setosa
4          4.6  setosa
5          5.0  setosa
6          5.4  setosa


irisのvirginicaだけを抜き出す方法

ise <- split(iris,iris[5])     #irisの5列目の項目で分割
head(ise$virginica)

    Sepal.Length Sepal.Width Petal.Length Petal.Width   Species
101          6.3         3.3          6.0         2.5 virginica
102          5.8         2.7          5.1         1.9 virginica
103          7.1         3.0          5.9         2.1 virginica
104          6.3         2.9          5.6         1.8 virginica
105          6.5         3.0          5.8         2.2 virginica
106          7.6         3.0          6.6         2.1 virginica



Sepal.Lengthの順に並べ換えます

Sepal.Lengthの昇順に並び替え
ise <- split(iris,iris[5])     #irisの5列目の項目で分割
isev <- ise$virginica
isevo <- isev[order(isev$Sepal.Length),]
head(isevo)

    Sepal.Length Sepal.Width Petal.Length Petal.Width   Species
107          4.9         2.5          4.5         1.7 virginica
122          5.6         2.8          4.9         2.0 virginica
114          5.7         2.5          5.0         2.0 virginica
102          5.8         2.7          5.1         1.9 virginica
115          5.8         2.8          5.1         2.4 virginica
143          5.8         2.7          5.1         1.9 virginica
````

ベイズの定理

ベイズの定理とは「観測値(データ)で条件付けられた、母数の分布を与える定理 」です

#Rのサンプル women$height を使用します
x <- women$height
#標本平均
me <- mean(women$height)
=65
#標本標準偏差
sd <- sqrt(sum((x - mean(x))^2) / length(x))
=4.32

  女性の伸長(インチ)の平均65、標準偏差4.32を正規分布の母数と見立ててしまえば、95%予測区間を導き出すことは簡単です.
  しかし、x~N(65,4.3)の分布が真の分布とは言えません.x~N(64,4.5)やx~N(66,4.6)などの母分散も考えられます.つまり、母数(平均や標準偏差)も分布することが考えられます.

確率密度関数 f(母数|観測値) の分布

例として母数を平均値, 標準偏差とします

f (θ|x) = f ( 平均値, 標準偏差 | x )

平均値, 標準偏差 とデータxは独立していないので条件付き分布となります

条件付き分布

f (x,θ) ≠ f(x)*f(θ)
# x が与えられた場合
  f (x,θ) = f(x)*f (θ|x)
# θ が与えられた場合
  f (x,θ) = f (x|θ)*f(θ)   


分布に関するベイズの定理

観測値(データ)で条件付けられた、母数の分布を与える定理

f (x,θ) = f(x)*f (θ|x)より
f (θ|x) = f (x,θ) ÷ f(x)
        = f (x|θ)*f(θ) ÷ f(x)

上記式より

f (θ|x) =事後分布,\ \ f(x|θ)=尤度,\ \ f(θ) =事前分布,\ \ f(x) =正規化定数

f (θ|x) = \frac{f (x|θ)*f(θ)}{ f(x)}



確率に関するベイズの定理

Aを得られた結果とします.Bをその原因とします.
例)5つの袋に赤白の玉が入り混ざっています.どれかの袋から玉が取られたとします.

取り出した玉=結果A
どの袋から取り出したか=原因B

 結果Aが得られた原因Bを推定します.袋が5つなのでB1, B2, B3, B4, B5 が原因となります. ベイズの定理では、結果Aであったときに原因が、B1またはB2またはB3またはB4またはB5 である確率を算出することになります.

ベイズの定理
P( Bi | A ) = \frac{P(Bi)P(A|Bi)}{P(A)} = \frac{P(Bi)P(A|Bi)}{ΣP(Bi)P(A|Bi)}

証明

A = A∩Ω =A∩(B1∪B2∪…∪Bi)
=(A∩B1)∪(A∩B2)∪...∪(A∩B1)
 P(A) = ΣP(Bi)P(A|Bi)

例)ベイズの逆確率
袋B1または袋B2から一つの球を取り出します.
B1には白玉3個、赤玉2個
B2には白玉1個、赤玉3個
取り出した球が赤だった場合(事象A)、袋がB1である確率、B2である確率を求めよ.

f:id:yoshida931:20180117184454p:plain:w300

正規化

f:id:yoshida931:20180316120813p:plain:w300


P(赤)= ΣP(Bi)P(赤|Bi)=P(B1)P(赤|B1)+P(B2)P(赤|B2)=\frac{1}{2}*\frac{2}{5}+ \frac{1}{2}*\frac{3}{4}

# 赤玉を引いた袋がB1である確率(ベイズの逆確率)   
  P(B1|A)= (1/2)*(2/5) / ((1/2)*(2/5) + (1/2)*(3/4)) = 0.3478261
# 赤玉を引いた袋がB2である確率(ベイズの逆確率)   
  P(B2|A)= (1/2)*(3/4) / ((1/2)*(2/5) + (1/2)*(3/4)) = 0.6521739
#正規化しているので
  P(B1|A) + P(B2|A) = 1

参考文献
豊田秀樹; はじめての統計データ分析 ベイズ的<ポストp値時代>の統計学, 朝倉書店, 2016

このブログの参考・引用文献

投稿日 2016-10-31
最終更新日 2017-12-29

文献:著者五十音順 巨人の肩に乗らねば・・・

石井一夫 ;Rとグラフで実感する生命科学のための統計入門, 羊土社 ,2017
石田 基広 ;改訂3版 R言語逆引きハンドブック ,シーアンドアール研究所; 改訂3版,2016

石田 基広 ;Rで学ぶデータ・プログラミング入門 ―RStudioを活用する―,共立出版 ,2012
石村卓夫;すぐわかる多変量解析,東京図書,1997
石村卓夫;入門はじめての時系列分析,東京図書,2012
石村卓夫;分散分析のはなし,東京図書,1992
稲葉 由之;プレステップ統計学I 記述統計学, 弘文堂, 2012
稲葉 由之;プレステップ統計学II 推測統計学, 弘文堂, 2013
北川 源四郎; 時系列解析入門, 岩波書店, 2005
金 明哲; Rによるデータサイエンス-データ解析の基礎から最新手法まで, 森北出版,2007
向後 千春, 冨永 敦子;統計学がわかる,技術評論社,2007
向後 千春, 冨永 敦子;統計学がわかる【回帰分析・因子分析編】,技術評論社,2008
田中 孝文; Rによる時系列分析入門, シーエーピー出版, 2008
東京大学教養学部統計学教室 (編集); 統計学入門 (基礎統計学), 東京大学出版会, 1991
東京大学教養学部統計学教室 (編集); 自然科学の統計学, 東京大学出版会, 1992

対馬栄輝; よくわかる医療統計 -「なぜ?」にこたえる道しるべ, 東京図書, 2015
豊田秀樹 (著, 編集);回帰分析入門 (Rで学ぶ最新データ解析) ,東京図書 ,2012
豊田秀樹; 検定力分析入門, 東京図書,2009
豊田秀樹; はじめての統計データ分析 ベイズ的<ポストp値時代>の統計学, 朝倉書店, 2016
西内 啓; 統計学が最強の学問である, ダイヤモンド社, 2013
西内 啓; 統計学が最強の学問である[実践編]---データ分析のための思想と方法, ダイヤモンド社, 2014
日本統計学会 (編集); 日本統計学会公式認定 統計検定2級対応 統計学基礎, 東京図書, 2012
南風原 朝和; 心理統計学の基礎―統合的理解のために, 有斐閣 , 2002
南風原 朝和; 続・心理統計学の基礎--統合的理解を広げ深める, 有斐閣 , 2014
舟尾 暢男; The R Tips―データ解析環境Rの基本技・グラフィックス活用集, オーム社; 第2版,2009
村上 正康,安田 正実; 統計学演習, 培風館,1989
柳川 堯 , 荒木 由布子; バイオ統計の基礎―医薬統計入門,近代科学社 ,2010
柳川 堯 ; 観察データの多変量解析―疫学データの因果分析,近代科学社 ,2016

山田 剛史, 杉澤 武俊, 村井 潤一郎; Rによるやさしい統計学, オーム社 , 2008

Hadley Wickham (著),石田 基広 (翻訳) ;R言語徹底解説, 共立出版,2016
R.A. フィッシャー(原著)、遠藤 健児、鍋谷 清治(訳); 実験計画法, 森北出版株式会社, 2013
R.A. フィッシャー(原著)、遠藤 健児、鍋谷 清治(訳); 研究者のための統計的方法, 森北出版株式会社, 2013
Peter Dalgaard (著), 岡田 昌史 (監修, 翻訳) ;Rによる医療統計学 原書2版,丸善,2007


ホームページ

おしゃべりな部屋 http://aoki2.si.gunma-u.ac.jp/
竹中明夫のページ  (http://takenaka-akio.org/index.html)

梶山喜一郎;コピペで学ぶ Rでテクニカルデータプレゼンテーション      (http://monge.tec.fukuoka-u.ac.jp/R_analysis/0r_analysis.html
間瀬 茂;R 基本統計関数マニュアル(https://cran.r-project.org/doc/contrib/manuals-jp/Mase-Rstatman.pdf)
続・わしのページ(R-Tips_舟尾 暢男 http://nfunao.web.fc2.com/
Quick-R (http://www.statmethods.net/)
R-Tips (http://cse.naro.affrc.go.jp/takezawa/r-tips/r2.html)
RjpWiki (http://www.okadajp.org/RWiki/)