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

統計学備忘録 since2016

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

相関係数の区間推定

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

 

身長のデータ
大学生
x<-c(172,167,184,175,176,175,170,180,170,179,167,175,174,162,165,163,170,169,165,175)
父親
y<-c(165,165,178,176,150,171,172,175,170,156,163,170,165,160,163,170,163,165,160,172)

 

cor.test関数を使用します
同じ数の要素をもつベクトル2つの相関関係、また相関係数の 95% 信頼区間,p値を求めます.Pearsonの相関係数をもとめるのであれば次の式を実行するのみでOKです.


  cor.test(x,y)

 

Pearson's product-moment correlation

data: x and y
t = 1.4587, df = 18, p-value = 0.1619
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.1371017 0.6711054
sample estimates:
cor
0.3251458

cor.test(x,y)の引数
x, y:同じ大きさの数字ベクトル
alternative(hypothesis): "two.sided", "greater", "less".
method : "pearson(相関係数cor)", "kendall"(相関係数tau), "spearman"(相関係数rho).

   (r-de様のコメントを参考に修正しました)

 

次に上記関数結果の95%信頼区間[-0.1371017 0.6711054]の求め方を忘れないように残しておきます

大学生の標準偏差
Vx<-(sum( (x-mean(x) )^2) )/20
Sx<-sqrt(Vx)
Sx

父親標準偏差
Vy<-(sum( (y-mean(y) )^2) )/20
Sy<-sqrt(Vy)
Sy

大学生と父親の共分散
Sxy<-(sum((x-mean(x) )*(y-mean(y) ) ) )/20

標本の相関係数r
Sxy/(Sx*Sy)
r=0.325

 

相関係数ρの区間推定

 

標本の相関係数rの分布を考えてみます.
r=0の場合
t=r*(sqrt(n-2)/sqrt(1-r^2) )は自由度n-2のt分布に従う

r≠0の場合(FisherのZ変換
rの分布はρが0から離れると非対称性が強くなります.
Z={log((1+r)/(1-r) ) }/2と変換した場合、Zの分布は対称になります.
またZは近似的にN({log( (1+ρ)/(1-ρ) ) }/2,1/(n-3) )に従うことが知られています.

 

r=0.325をZ={log((1+r)/(1-r) ) }/2に代入してZ=0.337を求めます.

{log((1+ρ)/(1-ρ) )=Rとおくと

(0.337-R)=±1.96×標準偏差
0.337-1.96/√(17)<R<0.337+1.96/√(17)
2×(0.337-1.96/√(17))<log((1+ρ)/(1-ρ)<2×(0.337+1.96/√(17))
‐0.138×2< log((1+ρ)/(1-ρ) )<0.671×2
したがってρの95%信頼区間
(exp(-0.276)-1)/(exp(-0.276)+1)
-0.137
(exp(1.626)-1)/(exp(1.626)+1)
0.671

 

また勉強して更新していきますのでコメントお願いします.

等分散性の検定

準備としてカイ二乗分布とF分布が必要になります

カイ二乗分布
確率変数Z1,Z2,…Znが互いに独立で標準正規分布にしたがうとき
W=Z1^2+Z2^2+…+Zn^2
の従う分布を自由度nのカイ二乗分布とよび、χ^2(n)で表します.
期待値E[W]=n, 分散V[E]=2n

カイ二乗分布の性質
① 二つの確率変数W1,W2が独立にχ^2(n1), χ^2(n2)に従うとき
W1+ W2は、χ^2(n1+n2)に従います.(カイ二乗分布の再生性)
② 母集団N(μ,σ^2)からの無作為標本X(X1,X2,…X3)について
Xiを標準すると Zi= (Xi-μ)/ σとなる
Ziは互いに独立で標準正規分布に従うので
∑{(Xi-μ)/ σ}^2 はカイ二乗分布χ^2(n)に従うことになります
③ 標本不偏分散
s^2={∑{(Xi-X平均) ^2 }/(n-1) より
(n-1) ×s^2=∑(Xi-X平均) ^2
母集団N(μ,σ^2)からの無作為標本Xi=(X1,X2,…X3)は互いに独立でN(μ,σ^2)に従うので
W=∑{(Xi-X平均) ^2}/ σ^2 = {(n-1) ×s^2}/σ^2
母平均μをX平均で置き換えると自由度がn-1になります
したがって{(n-1) ×s^2}/σ^2は自由度n-1のカイ二乗分布χ^2(n-1)に従います.

f:id:yoshida931:20170222184024p:plain

http://tips-r.blogspot.jp/2015/01/r.htmlより)


2標本の母分散の比の検定
散らばりを表す分散は非負の値なので、その大きさを考える場合には差ではなく比で考える方が適切です.

例)単位はmm
A群の身長の標準偏差はσa=5 (分散25)
B群の身長の標準偏差はσb=10 (分散100)
σa/σb=0.5 A群の方がB群に比べると2倍散らばりの小さい集団である
現実的には標準偏差が便利で理解しやすいのですが、推定での取り扱いは上述のような分散を使います.

 

F分布
条件
・確率変数U:自由度mのカイ二乗分布に従う
・確率変数V:自由度nのカイ二乗分布に従う
・UとVは独立である

フィッシャーの分散比 F= (U/m)/(V/n)

Fが従う確率分布を自由度(m,n)のF分布F(m,n)という

2標本の標本分散s1^2, s2^2を考えます
標本不偏分散s^2={∑{(Xi-X平均) ^2 }/(n-1) 
標本s1^2:{(m-1) ×s1^2}/σ1^2は自由度m-1のカイ二乗分布χ^2(m-1)に従います
標本s2^2:{(n-1) ×s2^2}/σ2^2は自由度n-1のカイ二乗分布χ^2(n-1)に従います
s1^2とs2^2は独立です

F={(m-1) ×s1^2}/(σ1^2)*(m-1)÷{(n-1) ×s2^2}/(σ2^2)*(n-1)
= (s1^2/σ1^2)÷(s2^2/σ2^2)
= (s1^2/s2^2)*(σ2^2/σ1^2)
の従う分布は、自由度(m-1,n-1)のF分布F(m-1,n-1)となります

 

 

それでは始めます
赤文字のみRで実行

2標本の場合

母分散の比の検定になります
二つの母集団からの標本をXi~N(μx,σx^2),Yi~N(μy,σy^2)とします.
帰無仮説:σx^2=σy^2、(σx^2)/(σy^2)=1
対立仮設:σx^2≠σy^2、(σx^2)/(σy^2)≠1

サンプルは「石村卓夫;分散分析のはなし,東京図書,1998,p69」
Xi<-c(12.2,18.8,18.2)
Yi<-c(26.4,32.6,31.3)
σx二乗<-sum((Xi-mean(Xi) )^2)/2;σx二乗
σy二乗<-sum((Yi-mean(Yi) )^2)/2;σy二乗

T<-σx二乗/σy二乗
1.25
F分布を使用します

qf(0.025,2,2,lower=F) # F 分布の上側5%点
39
qf(0.975,2,2,lower=F) # F 分布の下側5%点
0.026

T=1.25は棄却域には入らないので帰無仮説は棄却できない
標本をXiとYiの母分散は等しいと仮定できることになります

 

SPSSではlenene検定が使用されていますので以下に算出しておきます
Levene's Test もF分布を使用していますが、上記方法とは異なります
leveneTest(x,group)

Levene's Test for Homogeneity of Variance (center = median)
    Df  F value  Pr(>F)
group  1  0.0031  0.9585
    4

 

 

 

 

3標本以上の場合

 

繰り返し数が等しいとき:Hartleyの検定
各水準の分散(Si^2 )について
Si^2の最大分散÷最小分散が1に近づくことで等分散性を仮定します
最大分散/最小分散=max{Si^2}/min{Si^2}

 

f:id:yoshida931:20170222184102p:plain

http://www2.vmas.kitasato-u.ac.jp/lecture0/statistics/hartley.pdf


検定表の見方:Fmax(k , f)(k) = Fmax(水準数,繰り返し数-1)(水準数)
帰無仮説:各水準の分散は等しい

サンプルは「石村卓夫;分散分析のはなし,東京図書,1998,p108」
A<-c(12.2,18.8,18.2)
B<-c(22.2,20.5,14.6)
C<-c(20.8,19.5,26.3)
D<-c(26.4,32.6,31.3)
E<-c(24.5,21.2,22.4)

x<-c(A,B,C,D,E)

SuppDistsをインストトールします
library(SuppDists)
factor()カテゴリーを要素としたベクトル
group=factor(rep(c("A", "B", "C", "D", "E"), c(3, 3, 3, 3, 3)))
levels()でグループ化をstr()で構成を確認できます
fmax=max( by(x, group, FUN=var) ) / min( by(x, group, FUN=var) )  
5.70

by関数の引数(data,INDICES,FUN=,simplify=)
data:data.frameやmatrixもOK
INDICES:factorをsubsetに指定
FUN :引数dataに適用する関数を指定。ここでは分散varを指定。
simplify :結果のデータ形式。TRUE:scalar、FALSE:list or array。

 

積分布pmaxFratio
pmaxFratio( fmax, 繰り返し数の自由度, 水準数)
1-pmaxFratio(fmax, 2, 5)  #p値
0.8154272

帰無仮説は棄却されず、等分散性を仮定することができる

 

 

繰り返し数が異なるとき:Bartlettの検定
Rヘルプより
繰り返し数は異なりませんが同じサンプルで試行してみます
以下実施方法のみ記載
A<-c(12.2,18.8,18.2)
B<-c(22.2,20.5,14.6)
C<-c(20.8,19.5,26.3)
D<-c(26.4,32.6,31.3)
E<-c(24.5,21.2,22.4)

x<-c(A,B,C,D,E)
group=factor(rep(c("A", "B", "C", "D", "E"), c(3, 3, 3, 3, 3)))

bartlett.test(x,group)


Bartlett test of homogeneity of variances
data: x and group
Bartlett's K-squared = 1.2291, df = 4, p-value = 0.8733

 


Leveneの検定
最もよくみかける等分散性の検定だと思われます
水準の繰り返し数に縛りはありません

A<-c(12.2,18.8,18.2)
B<-c(22.2,20.5,14.6)
C<-c(20.8,19.5,26.3)
D<-c(26.4,32.6,31.3)
E<-c(24.5,21.2,22.4)

統計量(Wikipediaより)
W={(N-k)/(k-1)}×{∑Ni(Zi.-Z..)^2/∑∑(Zij-Zi.)^2}

数値はsample
k:比較する群の数水準
水準数i  1,2,3,4,5
繰り返しj  1,2,3
N:総サンプル数
Ni:第i群のサンプル数
自由度N-k=15-5=10
自由度k-1=5-1=4
Z.. 全データの平均値
Zi. i番目の水準内の平均値 
Yij:第i群のj番目の変数
Zij:| Yij – Yiの平均|、| Yij – Yiの中央値|


carパッケージをインストールします
x<-c(A,B,C,D,E)
group=factor(rep(c("A", "B", "C", "D", "E"), c(3, 3, 3, 3, 3)))

leveneTest(x,group)

Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 4 0.1256 0.9698
10

p値=0.9698で帰無仮説は棄却されず、等分散性を仮定した検定が可能になります

 

引用・参考
バイオインフォマティクス入門 http://bio-info.biz/
石村卓夫;分散分析のはなし,東京図書,1998

一元配置分散分析 (対応なし) F値の算出方法

更新日2017.3.13
更新日2017.6.22


        f:id:yoshida931:20170721165122p:plain

一元配置分散分析とは、一つの因子に基づく“3つ以上の母平均の差の検定”です

因子Factor:実験結果に影響を与えると考えられる要因
水準Lebel:因子をいくつかに分ける基準 (水準=群

例題)
治療因子に基づく、患者20名の満足度データを比較する(繰り返し数=5)
因子Factor=治療、水準Level=治療A, 治療B, 治療C, 治療D
因子<c(水準)としてベクトルを作成

まずイメージを一覧表にしておきます

A<- c ( 10 , 6 , 10 , 9 , 10 )
B<- c ( 10 , 5 , 5 , 12 , 4 )
C<-c ( 5 , 4 , 11 , 4 , 6 )
D<-c ( 9 , 5 , 2 , 3 , 1 )
xx<- data.frame ( A , B , C , D )
xx<- t ( xx )
colnames ( xx )<- c ( 1:5 )

f:id:yoshida931:20170622170844p:plain

繰り返し数(n)=5(順序が入れ替わっても問題はない)

 

ALL<- c ( A , B , C , D )  #全データ 
mean ( ALL )     #総平均
Levels<- c ( rep("A",5) , rep("B",5) , rep("C",5) , rep("D",5) )    # 水準
oneway.test ( ALL ~ Levels , var.equal = T )         #一元配置分散

One-way analysis of means
data: AllD and AllL
F = 2.5789, num df = 3, denom df = 16, p-value = 0.08979

それではF = 2.5789を求めてみます

 

分散分析の理解
分散分析は、観測データにおける変動(誤差変動水準内変動)と各要因による変動(水準間変動)に分解し効果を判定する方法です。(2元配置では交互作用による変動も考えます)。
1、 データの分散成分の平方和(総平方和全変動全平方和)を分解し、誤差による変動と要因の効果による変動を分離します。
2、 平方和を自由度で割った平均平方を算出します。
3、 統計量はF値になります。「要因の効果によって説明される平均平方(水準間変動)」を分子、「誤差によって説明される平均平方(水準内変動)」を分母にとった値が分散分析の統計量です。
4、 F検定。F分布により各効果の有意性を判定します。

以下、水準間変動水準内変動という表現を使用します

 

まず各水準の平均間に差がない場合を考えてみます
赤線は全平均
F = 0.30185,  num df = 3,  denom df = 16,  p-value = 0.8236

f:id:yoshida931:20170219143818p:plain
残っている変動は水準内の変動になります。つまり、「全変動=水準内の変動」ということになります。偏差の2乗和(平方和)を全2乗和、または全変動といいます。

逆に考えると水準間変動>水準内変動の場合には、「各水準間に差がある」と考えることができます。したがって「水準間変動÷水準内変動」が分散分析の統計量になります

ここで一元配置分散分析の数学モデル(X=x1,x2,…)を考えてみます
X = 総平均 + (水準平均―総平均) + (X―水準平均)
(X-総平均) = (水準平均―総平均) + (X―水準平均)

f:id:yoshida931:20170313111752p:plain


ここで差の平方をとって考えていきます
(X-総平均)^2 、 (水準平均―総平均) ^2 、(X―水準平均) ^2
例) 20個       4個         20個 
平方和∑にすると次の式が成り立ちます
n=各水準に含まれている繰り返し数(例題ではn=5)

∑(X-総平均)^2=n∑(水準平均―総平均) ^2 + ∑(X―水準平均) ^2

f:id:yoshida931:20170313112921p:plain

f:id:yoshida931:20170622172124p:plain

 

水準間変動  ( A ~ D で4水準なので df = 4 - 1 = 3  )
y<-5*{(mean(A)-mean(ALL))^2+
(mean(B)-mean(ALL))^2+
(mean(C)-mean(ALL))^2+
(mean(D)-mean(ALL))^2}
=66.15


群間平方和 (
水準間変動の平均平方和 )
水準間変動:n∑(水準平均―総平均) ^2
各水準平均は正規分布N( 総平均+(水準平均―総平均)、σ^2/n)に従います。
したがって∑(水準平均―総平均) ^2÷(σ^2/n)は自由度(水準数-1)のχ二乗分布に従います。
水準間平均平方=∑(水準平均―総平均) ^2 ÷ (水準数-1)
66.15÷(4-1) = 22.05


水準内変動   (A~Dで4水準, 繰り返い数は5回なので df=4×(5-1)=16 )

sum( (A-mean(A))^2)+
sum( (B-mean(B))^2)+
sum( (C-mean(C))^2)+
sum( (D-mean(D))^2)
=136.8

郡内平方和、残差平方和 ( 水準内変動の平均平方和 )
水準内平方和:∑(X―水準平均) ^2
それぞれの水準における∑(X―水準平均) ^2 ÷ σ^2 はn-1のχ二乗分布に従う統計量になります。したがって∑∑(X―水準平均) ^2 ÷ σ^2は自由度、水準数×(n-1)のχ二乗分布に従います。
水準内平均平方=∑∑(X―水準平均) ^2 ÷水準数×(n-1)
136.8÷{4☓(5-1)} = 8.55


ポイント
誤差を正規分布に従う確率変数と考えます。したがって例題の水準A標本(10,6,10,9,10)は、正規母集団(6.55+9.0、σ^2)からの大きさ5の標本ということになります。∑(A-6.55)^2/σ^2 は自由度4のχ二乗分布に従う統計量となります。

分散分析の統計量F
    水準間平均平方 ÷ 水準内平均平方
= 22.05 ÷ 8.55 = 2.5789


Rで分散分析表を求めてみます
anova(aov(ALL~Levels))

Analysis of Variance Table
Response: ALL
       Df Sum Sq Mean Sq F value Pr(>F)
Levels 3 66.15 22.05 2.5789 0.08979 .
Residuals 16 136.80 8.55
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

f:id:yoshida931:20170313123012p:plain

正規分布のグラフ

グラフ初歩の初歩

手探り状態でやってます

赤文字のみRで実行

 

正規分布
#1正規分布の乱数を100個生成
x<-rnorm(100)
#ヒストグラム,freq=Fで確率密度,ylimでy軸範囲
hist(x,freq=F,ylim=c(0,0.6))
#枠の処理、lty線種は実践、btyで左と下のみ
box(lty=1,bty = "l")
#確率密度のグラフ,add=Tで重ねる、x軸範囲-3~3
curve(dnorm(x),add=T,-3,3)


これらをオリジナルの関数として記憶させます
mynorm<-function(n,g)
{
x<-rnorm(n)
hist(x,freq=F,ylim=c(0,0.6),xlab=g)
box(lty=1,bty = "l")
curve(dnorm(x),add=T,-3,3)
}


mynorm(100,"n=100")

 

f:id:yoshida931:20170212153358p:plain

 

#四つのグラフを2行2列のまとめてみます
par(mfrow = c(2,2))
mynorm(30,"n=30")
mynorm(100,"n=100")
mynorm(300,"n=300")
mynorm(1000,"n=1000")

f:id:yoshida931:20170212153327p:plain

二項分布のグラフ

こんな説明でどうでしょうか?
  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