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

統計学備忘録 since2016

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

確率点と累積分布

確率点と累積分 

上側確率が0.025の場合

lower.tail
logical; if TRUE (default), probabilities are P[X ≤ x] otherwise, P[X > x].
                         Rのhelpより

正規分布の確率点(quantile)     qnorm( 確率 )

上側確率が0.025になるZ値
( p=0.975とおくことで上側確率は0.025となる )

qnorm ( 0.975 )  = 1.959964
                        = qnorm ( 0.025, lower.tail=TRUE ) 
下側確率が0.025になるZ値( 上側確率は0.975 )
qnorm ( 0.025 )  = - 1.959964

正規分布の累積分布     pnorm( 確率点 )
上側確率が0.025になる累積分
pnorm ( 1.96 )  =  0.9750021 
≒ pnorm ( qnorm ( 0.975 ) )
下側確率が0.025になる累積分
pnorm ( 1.96 )   =  0.0249979   ≒  pnorm ( qnorm ( 0.025 ) )

t分布の確率点( 自由度10 )
上側確率が0.025になるt値
qt ( 0.975 ,10  ) 

= 2.228139
下側確率が0.025になるt値
qt ( 0.025 ,10 )                       
=  - 2.228139

適合度の検定 (ピアソンの 𝜒^2 適合度検定)

適合度の検定 (ピアソンの 𝜒^2 適合度検定)
投稿日2016.11.4

更新日2017.7.21

理論上の確率分布から得られる期待度数 ( Expected frequency ) を求めることが前提となります.その期待度数を利用して観測値の度数 ( Observed frequency ) が適合するかどうか(当てはまりの良さ)を検定します.
f:id:yoshida931:20170714165249p:plain

帰無仮説の表現いろいろ

  • 観測度数の割合は期待度数の割合に適合している.
  • 各階級の発生確率は母集団の確率に適合する.
  • データは想定した確率分布に従う.

k.ピアソンの適合度基準
𝜒^2 =  Σ ( ( ( 観測度数 - 期待度数 )^2 ) / 期待度数 )

       =  sum ( ( Ok - Ek )^2 / Ek )
 

例)次のデータはメンデルの法則 1:2:1 の確率に適合しているか?
                        AA   Aa  aa   n(合計)
                観測度数    52  102  46  200
メンデルの法則(理論確率)     1/4  1/2  1/4         


観測度数 <-c( 52 , 102 , 46 )

理論確率 <- c ( 1/4 , 1/2 , 1/4 )        
期待度数 <-  期待確率 * sum ( 観測度数

50 100  50

     AA   Aa  aa   n(合計)
観測度数  52  102  46 

期待度数  50  100  50 

 

k.ピアソンの適合度基準
𝜒^2 =  Σ ( ( ( 観測度数 - 期待度数 )^2 ) / 期待度数 ) 
  =((52-50)^2/50)+((102-100)^2/100)+((46-50)^2/50)
        =0.44

p値
pchisq ( 𝜒^2  , 2 , lower.tail = FALSE )      # r-de-r様からの助言により修正!
pchisq ( 0.44 , 2 , lower.tail = FALSE )    
=0.8025188

lower.tail
logical; if TRUE (default), probabilities are P[X ≤ x] otherwise, P[X > x].

結果
p値=0.8025188となり帰無仮説は棄却されない
したがって、1:2:1の確率に適合していると言える

Rの関数で確認します
Ok< -c( 52 , 102 , 46 )          #観測度数 
Pk < - c ( 1/4 , 1/2 , 1/4 )      #理論確率 

chisq.test(Ok, p=Pk)

Chi-squared test for given probabilities

data: Ok
X-squared = 0.44, df = 2, p-value = 0.8025

回帰係数の区間推定

投稿日2017.2.3
更新日2017.7.21

回帰係数の区間推定

 

標本回帰係数の不偏推定量から母回帰係数の信頼区間を求めてみます.

東京 福岡
1019.4 1018.4
1005.7 1007.6
1002 1006.2
1006.7 1009.9
1005.1 1010.8
1010.1 1013.2
1016.7 1016.2
1011 1009.1
999.5 1003.1
1006.9 1012.5
1001.9 1006.4
1007.5 1006.3
1014.4 1012.2
1014.3 1015
1014.6 1017.4
1009 1016.5
1006.7 1012.1
1009.4 1008.7
1011.8 1009.2
1009.4 1009.2
    サンプル)統計学入門(基礎統計学),東京大学出版会,1991,p258

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

 Rstudioで確認

f:id:yoshida931:20170721144822p:plain

回帰分析 ( 東京=y、福岡=x )
lm(x$東京~x$福岡)

Call:
lm(formula = x$東京 ~ x$福岡)

Coefficients:
(Intercept)       x$福岡 
    16.0876       0.9822 

標本回帰方程式
Y=16.09+0.98X


定義
yi   実測値   x$東京
xi   実測値   x$福岡
B 回帰係数(最小二乗法より)   ・・・ この信頼区間を求める
β 回帰係数の不偏推定量 
ei 誤差項
Ei   (回帰)残差
Txx    xの偏差平方和      ∑(xi-xi平均値)^2
Txy xyの偏差平方和    ∑(yi-yi平均値)(xi-xi平均値)
Tyy    yの偏差平方和      ∑(yi-yi平均値)^2
σ^2 誤差の分散
s^2    誤差の分散の不偏推定量

回帰残差

母回帰方程式  yi = α + βxi + Ei  Ei~( 0 , σ^2 )  (i=1,2,...,n)
回帰残差  Ei = yi - α - βxi
Ei<- ( x$東京 - 16.09 - 0.9822 * x$福岡 )

誤差の分散の不偏推定量
回帰残差の平方和を ( n-2 ) で割る
回帰方程式の当てはまりの良さがわかります
s^2 <- ( sum (Ei^2) ) / ( n - 2 ) 
s^2 <- ( sum ( Ei^2 ) ) / 18

推定値の標準誤差
s<-sqrt ( s^2 )

ほとんどの ei  0 ± 2s の中に入入ります
残差の散布図で確認(言訳;グラフは初心者です)
plot ( ei , type="p" , ylim = c(-6,6) )
abline ( h=s , col='black' , lty="twodash" )
abline ( h = -s , col='black' , lty="twodash" )
abline( h= 2*s , col='black' ,lty="dashed" )
abline( h = -2*s , col='black' , lty="dashed" )


graphics.off()

偏回帰係数の標本分布
σは未知なのでsに置き換えます
Z = ( β - B ) ÷ σ / √Txx ~N ( 0 , σ^2 )
t = ( β - B ) ÷ s / √Txx    ~t ( n - 2 )

t値 
自由度18、p値0.025の確率点は
t値 = qt(0.975,18) = 2.101

Bの信頼区間
β - t × s / √Txx ≦ B ≦ β + t × s / √Txx
Txx<- sum ( ( x$福岡 - mean( x$福岡 ) ) ^2 )
0.98 - 2.101 * s / sqrt ( Txx )
0.6189584
0.98 + 2.101 * s / sqrt ( Txx )
1.341042
     0.62 ≦ B ≦ 1.34

Rの関数で確認します
confint ( lm ( x$東京 ~ x$福岡 ) , level = 0.95 )

              
(Intercept) -348.9125481 381.087689
x$福岡         0.6211872   1.343239

 

二項分布からの最尤推定

投稿日2017.6.13
更新日2017.7.20

二項分布から最尤推定に挑戦

二項分布のグラフ - 統計学備忘録 since2016

例)試行回数3回・確率1/2の二項分布
       𝐵𝑖(3 , 0.5)

f:id:yoshida931:20170613163424p:plain

 

表が出る確率をp=1/2とすると、次のような分布が考えられます.
確率pのときに表が〇回出る確率(同時確率

    f:id:yoshida931:20170613163444p:plain

これらの標本から得られる同時確率は、右に記載している確率です.
L(P) = 1回目表の確率×2回目表の確率×3回目表の確率


最尤原理
表が出る確率をP=?とした場合に、今現在起きている事象が起きる確率(尤度)が最も大きなる値が、その母数の推定値として「最も尤もらしい」とする考え方です.尤度が最大になる母数の値を最尤推定といいます.またそのような値を与える推定値を最尤推定といいます.

例)表の出る確率Pが未知で、表、表、ウラという事象を考えてみます.
L(P) = p × p × (1-p) = p^2 × (1-p) となります
ここで p=0.2 と p=0.8 を仮定してみます.
L(P) = p×p×(1-p) = p^2×(1-p)
L(0.2)=0.032
L(0.8)=0.128
最尤原理から考えると、L(0.2)<L(0.8)となり、p=0.2にくらべるとp=0.8が尤もらしい値となります.つまり、今起きている「表、表、ウラ」の尤度は、p=0.8が推定値として尤もらしいということになります.「表、ウラ、ウラ」となった場合には、p=0.2が尤もらしい推定値になります.

 L(P) 尤度関数 ( maximum function )
pのいろいろな値における尤もらしさを表す関数です.尤度関数L(P)を最大にする値が最尤推定値、関数としては最尤推定量になります.

求めるためには微分が必要になります.
L(P) = p×p×(1-p) = p^2×(1-p) を微分して=0としたときのpが最尤推定値となります.
dL(P)/dp = (p^2 – p^3)’ = 2p-3p^2 = 0
最尤推定値p=2/3
    ↓
確率密度関数Fp(y)をもつ母集団から、大きさ3の確率標本を取り出します
{ 1回目 , 2回目 , 3回目 }
同時確率密度:Fp(y1,y2,y3) = p×p×(1-p)

これをpの関数とみなしたものが尤度関数です
L(p) =  fp( y1 ) * fp( y2 )* ,…, * fp( yn )


つまり次ようになります
θを未知の母数とした確率密度関数Fθ(y)の母集団から
大きさ n の確率標本 { Y1 , Y2 ,…, Yn } を取り出した場合 
同時確率分布:Fθ( yn )
同時確率密度:Fθ( y1 , y2 ,…, yn ) = Fθ( y1 )* Fθ( y2 )* ,…, * Fθ( yn )
これをθの関数とみなしたものが尤度関数です
L(θ) = ΠFθ(y) = Fθ( y1 )* Fθ( y2 )* ,…, * Fθ( yn )

{ Y1 , Y2 ,…, Yn } を観測して { a1 , a2 ,…, an } が得られた場合、
尤度関数を最大にするθのことをY1= a1、Y2= a2、…、Yn= anが観測されたときのθの最尤推定値( maximum likelihood estimate )といいます.


例)
表の出る確率が p で{ Y1=表、Y2=表、Y3=ウラ }が観測された場合
標本から得られる確率は L(p) = p^2 * (1-p)
これを微分して=0としたpの値が最尤推定値となります
dL(p)/dp= (p^2 – p^3)’ = 2p-3p^2 = 0
p=2/3
表が2回観察されたときのpの最尤推定値は2/3


対数尤度 ( log likelihood function )

尤度関数を考えるとき、積の形なので数学的に扱うのが不便であるため対数をとります
log L(θ) = ∑ log { fθ(y1) *  fθ(y2) * ,…, *  fθ(yn) }

二項分布の例 
L(p) = Π (p^ai × (1-p)^n-ai  ) = p^T × (1-p) ^n-T    (T=∑ai)
logL(p)= log { p^T × (1-p) n-T ) }  =  T× logp + ( n-T ) × log( 1-p )
最尤推定値を求めるためにやはり微分して=0とします
d logL(p) / dp =  { T×logp + (n-T) × log(1-p) } = T / p – ( n - T ) / ( 1 - p ) = 0
p = T / n = ( 1/n ) *∑ai (相対頻度
例)p = ( 1/3 )*( 1+1 ) = 2 / 3
pの最尤推定値は相対頻度

参考文献
柳川 堯 , 荒木 由布子; バイオ統計の基礎―医薬統計入門,近代科学社 ,2010
東京大学教養学部統計学教室 (編集); 統計学入門 (基礎統計学), 東京大学出版会, 1991

Rstudioでの data.frame 表示の不具合

データフレームを表示したときに、ラベルが見えなくなりました

f:id:yoshida931:20170720121311p:plain

 

RStudio-1.0.143 から RStudio-1.0.153 へバージョンアップすることで不具合は解消しました.
https://www.rstudio.com/products/rstudio/download/preview/

f:id:yoshida931:20170720121437p:plain

 

Rstudioのサポートページの掲示板で助言を受けました!

https://www.rstudio.com/

f:id:yoshida931:20170720121521p:plain

ファイ係数

連関係数
 クロス集計表における2つの変数間の関連性の程度を表す指標として,連関係数が提案されています.

連関係数<|0.2| 連関はほとんどない
|0.2|≦連関係数<|0.4| やや連関がある
|0.4|≦連関係数<|0.7| 強い連関がある
|0.7|≦連関係数 かなり強い連関がある
   対馬栄輝; よくわかる医療統計 -「なぜ?」にこたえる道しるべ, 東京図書, 2015

 

ファイ係数 (φ係数)
 クロス集計表における行要素と列要素の関連の強さを示す指標. 1と0の2つの値からなる変数に対して計算される相関係数です。つまり、2つの2値変数をそれぞれ1と0で数値化した上で算出した積率相関係数です。1と0の2つの値とは、「好き、嫌い」、「はい、いいえ」など、2択で表現されるデータです. 選択肢Aに0を、選択肢Bに1を割り付けることで、ファイ係数を求めることができます.


例)2つの質問 ( q1 , q2 ) に関する答え ( yes , no ) のベクトルを考えます.
q1<-c("yes","yes","yes","yes","yes","yes","yes","no","no","no")
q2<-c("yes","yes","yes","no","no","no","no","yes","yes","no")
q101<-ifelse(q1=="yes",1,0)        # 0,1のベクトルに変更
q201<-ifelse(q2=="yes",1,0)       # 0,1のベクトルに変更

ピアソンの積率相関係数を2行×2列のクロス集計表に適用したもの
   ピアソンの積率相関係数

           f:id:yoshida931:20170718173706p:plain

mean ( q101 )
=0.7
mean ( q201 )
=0.5
sum ( ( q101-0.7 ) * ( q201 - 0.5 ) )
= - 0.5
sum ( ( q101 - 0.7 )^2 )
=2.1
sum ( ( q201 - 0.5 )^2 )
=2.5
r = - 0.5 / sqrt ( 2.1 * 2.5 )
  = - 0.2182179


Rの関数で確認してみます

cor(q101,q201)    
           
#相関係数をcor( )で算出
= - 0.2182179


次に例題をテーブルにして眺めてみます
tq<- table ( q1,q2 )
 
   

f:id:yoshida931:20170718172336p:plain

t<- tq [ c(2,1) , c(2,1) ]   
#これは不要ですがラベルの並べ替えです.yesから始めたいので挿入しておきます(1番目と2番目の入れ替え)

f:id:yoshida931:20170718172305p:plain

理解のためにデータフレームに変換します
t< - data.frame ( t ) 
     

f:id:yoshida931:20170720121016p:plain

 
     q1   q2    Freq
1   yes  yes   3         #  t $ Freq[1]
2   no    yes   2        #  t $ Freq[2]
3   yes   no    4        #  t $ Freq[3]
4   no    no     1        #  t $ Freq[4]

ファイ係数φその1 カイ二乗値より
φ=√(カイ二乗値/n)
chisq.test ( tq , correct=F )
X-squared = 0.47619
φ=sqrt ( 0.47619 / 10 )
= 0.2182178

ファイ係数φその2  観測値と周辺度数より

f:id:yoshida931:20170718213849p:plain

φ =  ( a*d - b*c ) / √ ( a+b )*( c+d )*( a+c )*( b+d )*
a<- t$Freq[1] * t$Freq[4] - t$Freq[2] * t$Freq[3]   # a=対角線上の観測値を掛けたものの差
b<-( t$Freq[1] + t$Freq[3] )*                                 # b=全ての周辺度数を掛けたもの
  ( t$Freq[2] + t$Freq[4] )*
  ( t$Freq[1] + t$Freq[2] )*
  ( t$Freq[3] + t$Freq[4] )
φ=a/sqrt(b)
= - 0.2182179

 

Rの関数で確認してみます
install.packages ( "psych" )
library ( psych ) 
phi ( tq , digits=8 )
= - 0.2182179

イェーツの補正

イェーツの補正 / イェーツの連続修正(Yate's continuity correction)
 離散型分布を連続型分布に近似させて統計的検定を行う際に使用する修正です.2×2分割表のデータに対して行われ、より正確な検定が可能になります.「連続的なカイ二乗分布(自由度1のχ二乗分布)」に近似するために、2×2分割表の中で観察される「二項分布型度数の離散型の確率」を補正します.各々の観測値と期待度数の差の絶対値より0.5を差し引くことによりカイ二乗検定の式を調整します.計算の結果得られるカイ二乗値を減らすことになりp値が増加することになります.

 

例)
x<-c(116,244)
y<-c(76,44)
xy<-data.frame(x,y)

 

期待度数
sum(x) * sum(116+76) / sum(xy)
sum(y) * sum(116+76) / sum(xy)
sum(x) * sum(244+44) / sum(xy)
sum(y) * sum(244+44) / sum(xy)
xe<- c( sum(x) * sum(116+76) / sum(xy) , sum(x) * sum(244+44) / sum(xy) )
ye<- c( sum(y) * sum(116+76) / sum(xy) , sum(y) * sum(244+44) / sum(xy) )
xye<-data.frame(xe,ye)

 

カイ二乗値
( 116 - xye[1,1] ) ^2 / xye [1,1] +
 ( 76 - xye [1,2] ) ^2 / xye [1,2] +
 ( 244 - xye [2,1] ) ^2 / xye [2,1] +
 ( 44 - xye [2,2] ) ^2 / xye [2,2]
 = 35.01157

 

イエーツの補正 を実施したカイ二乗値
 ( abs ( 116 - xye[1,1] ) - 1 / 2 ) ^2 / xye[1,1] +
  ( abs ( 76 - xye[1,2] ) - 1 / 2 ) ^2 / xye[1,2] +
  ( abs ( 244 - xye[2,1] ) - 1/2 ) ^2 / xye[2,1] +
  ( abs ( 44 - xye[2,2] ) - 1/2 ) ^2 / xye[2,2]
 = 35.01157

 

Rの関数で確認します
chisq.test(xy)

Pearson's Chi-squared test with Yates' continuity correction
data:  xy
X-squared = 35.012, df = 1, p-value = 3.278e-09

 

サンプルは日本統計学会 (編集); 日本統計学会公式認定 統計検定2級対応 統計学基礎, 東京図書, 2012,p156の例7