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

統計学備忘録 since2016

Rを使って統計学を勉強するブログです

二項分布のグラフ

こんな説明でどうでしょうか?
  f:id:yoshida931:20170524161843p:plain

コイン投げ(いかさまなし)
表の出る確率確率1/2
3回コイン投げを実施してみた場合・・・
表が出る回数(x)は0回、1回、2回、3回
xの確率は𝐵𝑖(3 , 0.5) に従う.
その確率分布を二項分布といいます.

    f:id:yoshida931:20170524161933p:plain

f:id:yoshida931:20170524161906p:plain

 

 

村上 正康,安田 正実; 統計学演習, 培風館,1989,p51のグラフをRで作ってみます
B(10,0.2),B(10,0.5),B(10,0.8)を重ねたグラフです
赤文字のみRで実行してください

 (foo.bar@baz.com様よりコメントいただき修正しております。簡単で理解しやすい方法です。)

t0.2 <- 0:6
P0.2 <- dbinom(t0.2, 10, 0.2)
t0.5 <- 1:9
P0.5 <- dbinom(t0.5, 10, 0.5)
t0.8 <- 4:10
P0.8 <- dbinom(t0.8, 10, 0.8)
plot(t0.2, P0.2, type="b", lty=2, xlab="", ylab="", xlim=c(0,10), ylim=c(0,0.4))
lines(t0.5, P0.5, type="b", lty=2)
lines(t0.8, P0.8, type="b", lty=2)

f:id:yoshida931:20170121213325p:plain

 

 

 

途中経過を忘れないために残しておきます

 
n=10 、確率0.2 
t0.2 <- 0:6;p<-0.2;n<-10
y<-factorial(n)*(p^t0.2)*(1-p)^(n-t0.2)
r<-factorial(t0.2)*factorial(n-t0.2)
y/r

n=10、確率0.5
t0.5 <- 1:9;p<-0.5;n<-10 
c<-factorial(n)*(p^t0.5)*(1-p)^(n-t0.5)
d<-factorial(t0.5)*factorial(n-t0.5)

n=10、確率0.8
t0.8 <- 4:10;p<-0.8;n<-10 
e<-factorial(n)*(p^t0.8)*(1-p)^(n-t0.8)
f<-factorial(t0.8)*factorial(n-t0.8)
e/f


p=0.2,0.5,0.8の結果をそれぞれ算出
間違わないように名前付け
P0.2<-(y/r)
P0.5<-(c/d)
P0.8<-(e/f)

linesでグラフを重ねる

plot(t0.2,P0.2,type="b",lty=2,xlab="",ylab="",xlim=c(0,10),ylim=c(0,0.4))
lines(t0.5,P0.5,type="b",lty=2)
lines(t0.8,P0.8,type="b",lty=2)

 

グラフに色をつける

library(RColorBrewer)
#RColorBrewerパッケージのサンプル
display.brewer.all()
 
ヒストグラムを塗ってみる
どの色セットを使用するかを指定する
cols <- brewer.pal(8,"Pastel1")   # brewer.pal(何色、パレット名)

y<-c(1,2,3,4,5,6,7)
p<-c(2,3,4,5,4,3,2)
q<-c(2,3,4,5,4,3,2)
par(mfrow = c(3,3),mar = c(5, 4, 1, 4)) #余白 底辺、左、上、右の順 
pos.x <- barplot(q,ylim=c(0,6),col=cols[1])
pos.x <- barplot(q,ylim=c(0,6),col=cols[2])
pos.x <- barplot(q,ylim=c(0,6),col=cols[3])
pos.x <- barplot(q,ylim=c(0,6),col=cols[4])
pos.x <- barplot(q,ylim=c(0,6),col=cols[5])
pos.x <- barplot(q,ylim=c(0,6),col=cols[6])
pos.x <- barplot(q,ylim=c(0,6),col=cols[7])
pos.x <- barplot(q,ylim=c(0,6),col=cols[8])
pos.x <- barplot(q,ylim=c(0,6),col=cols[9])

 

f:id:yoshida931:20170114232836p:plain

 

棒グラフの中央に散布図をプロットする方法 (pos.x)

y<-c(1,2,3,4,5,6,7)
p<-c(2,3,4,5,4,3,2)
q<-c(2,3,4,5,4,3,2)
par(mar = c(5, 4, 1, 4)) #余白 底辺、左、上、右の順 
pos.x <- barplot(q,ylim=c(0,6))
points(pos.x, p)

 

f:id:yoshida931:20170113120942p:plain

Rことはじめ

1、Rのインストール方法と使い方

2、以下のデータの母平均の95%信頼区間を求めてみます

73.0 75.8 74.0 72.3 76.9 82.5 69.1 80.9 74.9 71.6 73.5 84.6 81.3 82.5
78.1 71.3 76.3 74.9 74.8 72.5 78.5 71.5 70.0 81.1 76.6 77.9 75.5 70.8
81.9 75.1 85.4 82.0 73.4 76.9 78.5 77.0

 以下のファイルをご覧ください

Rの使い方.pdf - Google ドライブ

 

 

サンプル

対馬栄輝研究室 - 弘前大学医学部保健学科・大学院保健学研究科

http://www.hs.hirosaki-u.ac.jp/~pteiki/

統計解析→Rコマンダー

 

級内相関係数 ICC

修正日 2017.3.2
修正日 2017.6.21


ICC (Intraclass correlation coefficients)

ICCは信頼性の指標の一つとして利用されています.
ICCには3つの分類があり、6Caseに分類されています.

f:id:yoshida931:20170621175544p:plain

注)ICCの記載方法は統一されていない

ICCは分散分析の結果を活用しています
一元配置分散分析の重要な仮定
水準Aiに対してXij = μ+ Ai + Eijというモデルを考えた場合
μ+ Ai=水準Aiの母平均、Eij=誤差
条件1:標本は無作為抽出
条件2:誤差は母平均0の正規分布に従う(標本は正規分布に従う)
条件3:水準間の母分散はすべて等しい(等分散性の検定が必要)
ICCの場合、条件3は必ず保証されなければいけないものではない(対馬

一般可能性理論
一般可能性理論は一般化可能性研究(generalizability study:G研究)と決定研究(decision study:D研究)からなります.G研究は、検査の測定誤差に着目しており、その変動要因の成分とバラツキの大きさの推定値を求めます.そして各変動要因や交互作用について検討していく研究方法です.D研究では、検査の信頼性を検討します.G研究で得られた分散分析の推定値を利用し、一般化可能性係数を算出します.そして、どの程度の評価項目や評価者が必要なのかを検討する研究です(山西2005).分散分析では、水準の分散と誤差分散の推定値の比を見ることで、信頼性を推定します。これは一般化可能性理論と考えられます.

f:id:yoshida931:20170621195330p:plain

 

ICCは分散分析から求められる平均平方和の期待値を利用します.
全ては下記の検査結果の一覧表から始まります.

f:id:yoshida931:20170302184239p:plain


各自由度
被験者間:n—1
被験者内:n(k-1)
検者間:k-1
残差:(n-1)(k-1)

Case1の例
1名の検者で複数の被験者にそれぞれ数回検査を実施します.検者の誤差は考慮しません.ICC ( 1,1 ) は一回の測定に関する検者内信頼性、ICC ( 1,k ) はk回繰り返し測定した平均値に対する検者内信頼性. 

f:id:yoshida931:20170626153155p:plain

Case2の例
複数の検者で複数の被験者にそれぞれ数回検査を実施します.検者の誤差が発生します.ICC ( 2,1 ) はランダム効果による、つまり任意の検者に対する検者間信頼性.ICC ( 1,k ) の1回目、2回目、3回目を検者に置き換えると、k人の検者がn人の患者を検査した2元配置 ( 繰り返しなし ) の結果となる.

f:id:yoshida931:20170626153022p:plain
ICC ( 2,k )  は、複数(k人:3人)の検者が、複数(5人)の被験者を複数回(4回)繰り返した平均値に対する検者間信頼性.

f:id:yoshida931:20170626153428p:plain

Case3の例
複数の検者で複数の被験者にそれぞれ数回検査を実施します.検者が混合モデルとなります.ICC ( 3,1 ) は、固定効果による、つまり測定に参加した検者群そのものが対象となる.測定の精度(一致性:agreement)というよりも整合性(一貫性:consistency)を確かめるもの.
f:id:yoshida931:20170626161532p:plain

 

Rによるケース2(Shrout&Fleiss,1979)の算出方法
赤文字のみRにペーストして実行するのみ
二要因変量モデル

ICC(2,1)のRによる出力は以下のような分類になります McGraw & Wong (1996)
ICC(A,1) 一致性を示す級内相関係数(検者間の信頼性)
ICC(C,1) 一貫性を示す級内相関係数(検者内の信頼性)

準備
行列の配置を確認
行:被験者(測定の対象)
列:検者(評定者)
irrを読み込む
library(irr)

村山航;級内相関係数についての覚書p4のサンプル
例1
x<-c(1,2,3,4,5)
y<-c(1,2,3,4,5)
例2
x<-c(1,2,3,4,5)
y<-c(5,6,7,8,9)
例3
x<-c(1,2,3,4,5)
y<-c(3,5,7,9,11)

行列形式に変換cbind()してiccを算出

 

ICC(A,1) 一致性を示す級内相関係数
icc(cbind(x,y), model ="twoway",type = "agreement")
ICC(C,1) 一貫性を示す級内相関係数
icc(cbind(x,y), model ="twoway",type = "consistency")

  

ICCが0.7以上であれば信頼性は良好と考える
ICCは点推定であるので区間推定を提示することが重要

0.0-0.20    slight
0.21-0.40   fair
0.41-0.60   moderate
0.61-0.80   substantial
0.81-1.00   almost perfect

 

 

例題
対馬栄輝;頼性指標としての級内相関係数p.23の膝屈曲のサンプル

A B C D
126 122 131 125
137 143 141 141
113 119 115 105
153 143 135 144
146 157 150 149
161 157 160 160
110 109 105 113
145 151 152 156
126 141 132 122
114 126 130 125

x<-read.table("clipboard",header=T)

 

Case1 :一要因モデル (A=1回目、B=2回目、C=3回目、D=4回目とする)
検者内信頼性 ( Intra-rater reliability )

膝角度計測の検者信頼性を求める
一人の理学療法士が10人の患者に対してそれぞれ4回膝の角度を計測し信頼性を求める

ICC ( A,1 ) = ICC ( 1 , 1 ) 一致性を示す級内相関係数

icc(x, model ="oneway",type = "agreement")

Single Score Intraclass Correlation

Model: oneway
Type : agreement

Subjects = 10
Raters = 4
ICC(1) = 0.909

F-Test, H0: r0 = 0 ; H1: r0 > 0
F(9,30) = 40.9 , p = 2.06e-14

95%-Confidence Interval for ICC Population Values:
0.788 < ICC < 0.973


参考)

繰り返しのある一元配置分散分析
       Df Sum Sq Mean Sq F value Pr(>F)

被験者内 treat 3 76 25.4 0.894 0.457
被験者間 pat 9 10319 1146.6 40.421 2.25e-13 ***

ICC(A,1) =
    ( 群間平均平方ー 郡内平均平方 ) ÷ ( 群間平均平方 + ( 繰り返し数 - 1 )×郡間平均平方 )
     = ( 1146.6 - 25.4 ) ÷ ( 1146.6 + ( 4 - 1 ) × 25.4  ) = 0.916912

 

4回づつ測定した10名分の平均値から検者内信頼性 ( Intra-rater reliability )を求める
ICC ( C,1 )ICC ( 1,k ) 一貫性を示す級内相関係数
icc(x, model ="oneway",type = "consistency")
結果

Single Score Intraclass Correlation

Model: oneway
Type : consistency

Subjects = 10
Raters = 4
ICC(1) = 0.909

F-Test, H0: r0 = 0 ; H1: r0 > 0
F(9,30) = 40.9 , p = 2.06e-14

95%-Confidence Interval for ICC Population Values:
0.788 < ICC < 0.973

 

 Case2 :二要因モデル (4人の理学療法士をA,B,C,Dとする)
検者間信頼性 ( Inter-rater reliablity )

膝角度計測の検者信頼性(ばらつきの程度)を求める
4人の理学療法士が10人の患者に対してそれぞれ1回膝の角度を計測し信頼性を求める
ICC ( 2 , 1 ) 一致性を示す級内相関係数
icc(x, model ="twoway", type = "agreement")

Single Score Intraclass Correlation

Model: twoway
Type : agreement

Subjects = 10
Raters = 4
ICC(A,1) = 0.909       #一致性を示す級内相関係数

F-Test, H0: r0 = 0 ; H1: r0 > 0
F(9,30) = 40.4 , p = 2.46e-14

95%-Confidence Interval for ICC Population Values:
0.788 < ICC < 0.973

 

膝角度計測の検者信頼性(ばらつきの程度)を求める
4人の理学療法士が測定した10人の患者の平均値から信頼性を求める
ICC ( C,1 ) = ICC ( 2 , 1 ) 一致性を示す級内相関係数
icc(x, model ="twoway", type = "consistency")

Single Score Intraclass Correlation

Model: twoway
Type : consistency

Subjects = 10
Raters = 4
ICC(C,1) = 0.908         #一貫性を示す級内相関係数

F-Test, H0: r0 = 0 ; H1: r0 > 0
F(9,27) = 40.4 , p = 2.25e-13

95%-Confidence Interval for ICC Population Values:
0.782 < ICC < 0.973


ICC全てを一気に計算する
psychをインストールする
library ( psych )

ICC ( x )                        #大文字

Call: ICC(x = kneef)

Intraclass correlation coefficients
type ICC F df1 df2 p lower bound upper bound
Single_raters_absolute  ICC1   0.91  41  9  30  2.1e-14  0.79  0.97
Single_random_raters ICC2   0.91  40   9  27   2.3e-13   0.79   0.97
Single_fixed_raters ICC3   0.91   40   9   27  2.3e-13   0.78   0.97
Average_raters_absolute ICC1k   0.98   41   9   30   2.1e-14   0.94   0.99
Average_random_raters ICC2k   0.98   40   9   27   2.3e-13   0.94   0.99
Average_fixed_raters ICC3k   0.98   40   9   27   2.3e-13   0.93   0.99

Number of subjects = 10 Number of Judges = 4

 

反省
分散分析を復習してから再更新した方がよさそうです

Rで簡単(t検定の検定力分析)

パッケージpwrを使用します
library(pwr)

pwr.t.test(n = , d = , sig.level = , power = , type = c("two.sample", "one.sample", "paired"))
n=標本数
d=効果量()
sig.level=有意水準
power=検定力
type=t検定の形("two.sample=対応なし", "one.sample=1標本", "paired=対応あり")

n, d, power, and sig.level のいずれかを空白(NULL)にします
strict=T   #厳密な検定が可能になる
 (効果量(effect size)のはなし - 六本木で働くデータサイエンティストのブログ



必要な標本数を求めてみます
pwr.t.test(n =NULL, d =0.2 , sig.level =0.05 , power = 0.8, type ="two.sample")

Two-sample t test power calculation

n = 393.4057
d = 0.2
sig.level = 0.05
power = 0.8
alternative = two.sided

 

効果量0.2、有意水準0.05、検定力0.8の場合の標本の大きさは393必要になります
もちろんフリーソフトG*powerを使用した結果と同じになります

f:id:yoshida931:20170121182510p:plain

 


効果量をもとめてみます
pwr.t.test(n =150, d =NULL , sig.level =0.05 , power = 0.8, type ="two.sample")

Two-sample t test power calculation

n = 150
d = 0.3245459
sig.level = 0.05
power = 0.8
alternative = two.sided

 

標本数150、有意水準0.05、検定力0.8の場合の効果量は0.32と低くなります
フリーソフトG*powerを使用した結果

f:id:yoshida931:20170121183731p:plain


Rの予備知識
for(i in c(2,4,6)) でiに2,4,6を順に挿入する命令
print(・・・) 引数が数値ベクトルの場合にはベクトルの中身を表示する
print(c(1,2,3))
合わせて
for(i in c(0.26,0.24,0.22)){print(x<-c(i+1))}
round 四捨五入を導入して・・・
for(i in c(0.26,0.24,0.22)){print(round(x<-c(i+1),1))}

上記を利用して
有意水準が.05,.01,.001)のときの対応なしt検定のpowerを一気に算出してみます
pwr.t.test(n =NULL,d =0.2,sig.level=0.05,power=0.8,type ="two.sample")
for(i in c(.05,.01,.001)){print(round(pwr.t.test(n =394,d =0.2,sig.level=i,
,type ="two.sample")$power,2))}

[1] 0.8
[1] 0.59
[1] 0.31


対応のなし(2群の標本数が異なる場合)
pwr.t2n.test(n1 = , n2= , d = , sig.level =, power = )
例)
pwr.t2n.test(n1 =120, n2=NULL , d =0.5 , sig.level =0.05, power =0.8 )

t test power calculation

n1 = 120
n2 = 43.21666
d = 0.5
sig.level = 0.05
power = 0.8
alternative = two.sided

 

n1=120、効果量0.5、有意水準0.05、検定力0.8、両側検定の場合n1=43.2必要

Rで簡単(t検定:対応あり)

例えばn=14の前後比較・・・

下のデータをコピーして

 

前 後

90 98

90 95

92 77

91 84

83 75

93 88

89 89

81 88

80 86

77 86

78 79

77 90

76 88

79 99

 

Rに貼り付け

x<-read.table("clipboard",header=T)

 

対応のあるデータなので「差」を検定することになる

差<-c(x$後-x$前)

 

帰無仮説の下で、この差は自由度13のt分布に従うことから

t.test(差)

 

One Sample t-test

data: 差
t = 1.2899, df = 13, p-value = 0.2195
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
-2.217176 8.788605
sample estimates:
mean of x
3.285714

 

帰無仮説(差がない)が正しいという仮説のもとでp値は0.2195

「起こりにくい」とは言えないですね!

 

次の式でも同様に求められます

t.test(x$後,x$前,paired=TRUE)

 

99%信頼区間を求めたいなら・・・

conf.level=0.99

t.test(x$後,x$前,paired=TRUE,conf.level=0.99) 

 

> t.test(x$後,x$前,paired=TRUE,conf.level=0.99)

Paired t-test

data: x$後 and x$前
t = 1.2899, df = 13, p-value = 0.2195
alternative hypothesis: true difference in means is not equal to 0
99 percent confidence interval:
-4.387154 10.958582
sample estimates:
mean of the differences
3.285714