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

読者です 読者をやめる 読者になる 読者になる

統計学備忘録 since2016

Rを使用した統計学

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

(2017.3.13更新)

赤文字の部分だけRにペーストして実行

一元配置分散分析とは、一つの因子に基づく“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)

f:id:yoshida931:20170313105833p: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

 

全変動=∑∑(X-総平均)^2
水準間変動=n∑(水準平均―総平均) ^2
水準内変動=∑∑(X―水準平均) ^2

 

水準間変動
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


水準内変動

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=水準間平均平方÷水準内平均平方
{∑(水準平均―総平均) ^2 ÷(水準数-1)} ÷ {∑∑(X―水準平均) ^2 ÷ n(水準数-1)}
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

回帰係数の区間推定

赤文字をRで実行

 

サンプル)統計学入門(基礎統計学),東京大学出版会,1991,p258

次のデータをコピーしてください
東京 福岡
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

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

 

回帰分析
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
σ二乗 誤差の分散
s二乗    誤差の分散の不偏推定量

 


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

 

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

 

推定値の標準誤差
s<-sqrt(s二乗)

 

ほとんどの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,σ二乗)
t=(β-B)÷s/√Txx ~t(n-2)

 

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

 

単回帰分析

赤文字のみRで実行

 

線形回帰分析
説明変数と目的変数を直線関係で傾向を示す。
説明変数と目的変数との関係を直線でモデル化する回帰分析。

非線形回帰分析
非線形関係でモデル化する回帰分析。
説明変数と目的変数を非直線的関係で傾向を示す。

 

ここでは単回帰分析について解説します
以下のように定義します
観測値 ( xi , yi )
観測値の平均値 ( ya , xa )
予測式 y'= b0 + b1*xi

最小二乗法を用いて回帰係数を求めます
残差の二乗和(残差平方和)を最小にする回帰係数をもとめます

残差ei ( yi - y' )
残差平方和 Q( b0 , b1 )

Q(b0,b1) = ∑(yi-y')^2 = ∑{yi-(b0+b1*xi)}^2    最小二乗基準

Q(b0,b1) を最小にするb0とb1を求めることで回帰直線を得る

Q(b0,b1) はパラメータb0b1の二次関数で最小値が存在します
そこでb0b1偏微分し、0に等しいとします

 

f:id:yoshida931:20170317200530p:plain
Y=b0+b1*x1
この式をyのx上への回帰方程式、あるいは回帰直線といいます

b0 : y切片
b1 : 偏回帰係数

 

例)
統計学入門 ,東京大学出版会,p42,1991より
yi<-c(114,124,143,158,166)
xi<-c(35,45,55,65,75)

 

平均を求めます
(ya<-mean(yi))
141
(xa<-mean(xi))
55

 

偏回帰係数
B<-sum( (yi-ya)*(xi-xa) )/sum( (xi-xa)^2)
1.38

 

Y切片
A<-ya-1.38*55
65.1

y=65.1+1.38x

 

平方和の分解
次に平方和を分解してみます
観測値の平方和(または総平方和)をSyとします
Sy = ∑(yi-ya)^2
   = ∑{(yi-y')+(y'-ya)}^2
   = ∑(yi-y')^2+∑(y'-ya)^2+∑(yi-y')*(y'-ya)
となります
ここで∑(yi-y')*(y'-ya)を考えてみます
予測値と残差の相関は0になります(証明は省略)
ry'ei=Sy'ei/Sy'*Sei=0
Sy'ei=∑(yi-y')*(y'-ya)=0となります

したがって
     Sy     = ∑(yi-y')^2 + ∑(y'-ya)^2
     Sy     =    Sei      +    SR
観測値の平方和 =  残差平方和   +  回帰による平方和

 

決定係数
回帰直線で説明できる部分の割合
R^2 = SR / Sei
Se=0,R^2=1 完全に説明できる
SR=0,R^2=0 全く説明できない

さきほどの例の決定係数を求めてみます

決定係数<-sum( (65.1+1.38*xi-ya)^2)/sum( (yi-ya)^2)
0.9836777

 

Rの関数を使って確認してみます
summary(lm(formula=yi~xi) )

 

Call:
lm(formula = yi ~ xi)

Residuals:
   1    2    3    4    5
 0.6 -3.2  2.0  3.2 -2.6

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  65.1000     5.8284   11.17 0.001538 **
xi            1.3800     0.1026   13.45 0.000889 ***
---
Signif. codes: 
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.246 on 3 degrees of freedom
Multiple R-squared:  0.9837,    Adjusted R-squared:  0.9782
F-statistic: 180.8 on 1 and 3 DF,  p-value: 0.0008894

 

参考) 豊田秀樹 (著, 編集);回帰分析入門 (Rで学ぶ最新データ解析) ,東京図書 ,2012

t検定のまとめ

Rでは次の式を使用します
t.test(x, alternative=c("two.sided","less","greater"),
mu=0, paired=FALSE, var.equal=FALSE, conf.level=0.95)

引数
alternative:両側検定("two.sided")か片側検定("less","greater")か
mu=0:母平均が0であるかどうかを検定する
paired=TRUE:対応のあるt検定
var.equa=TRUE:等分散を仮定する
conf.level:信頼区間

統計量
t=(標本平均-母平均)/√(不偏分散/n) 
=(標本平均-母平均)/標準誤差  
自由度n-1のt分布t(n-1)に従う 

 

1標本問題(1つの母集団の母数に関する検定)
例)次のデータを図で概観し,母平均=160.0を有意水準5%で検定せよ
x = c(153.0,148.0,161.0,166.0,160.0,173.0,173.0,165.0,153.0,168.0,171.0)

箱ひげ図
boxplot(x)
ドット図
stripchart(x, pch=16, at=0, method="stack")
ヒストグラム
hist(x, right=FALSE, col="gray")

t検定を実行
t.test(x)

One Sample t-test
data: x
t = 62.79, df = 10, p-value = 2.554e-14
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
157.0405 168.5959
sample estimates:
mean of x
162.8182

標本xの平均は162.82
95%信頼区間[157.04, 168.60]で「160.0」が含まれている
したがって帰無仮説は5%有意水準で採択される(p値<0.01)

もし、この検定の有意水準が1%であった場合は
t.test(x,conf.level=0.01)

One Sample t-test
data: x
t = 62.79, df = 10, p-value = 2.554e-14
alternative hypothesis: true mean is not equal to 0
1 percent confidence interval:
162.7849 162.8515
sample estimates:
mean of x
162.8182

99%信頼区間[162.78, 162.85]で「160.0」が含まれていません
したがって帰無仮説は1%有意水準では棄却されることになります
研究前のαの設定は重要です

 

対応のある2標本の場合
同じ原理です
治療前後の比較を考えてみましょう
治療前のデータPRE,治療後のデータPOSTとします
治療の前後差を有意水準1%で検定せよ
PRE<-c(25,27,32,37,38,39,44,45,47,48,52,54,59,60,65,65,67,68,68)
POST<-c(41,28,32,60,50,48,44,55,95,48,61,67,60,83,66,77,70,81,87)
対応あるt検定の帰無仮説は「前後差=0」と設定します
def<-POST-PRE
つまり1標本となります
t.test(def)

One Sample t-test
data: def
t = 4.1537, df = 18, p-value = 0.0005965
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
5.540338 16.880715
sample estimates:
mean of x
11.21053

差の平均は11.21となり有意水準1%で差を認める
(当然この11.21という数値は検査内容により実質的な意味合いが異なります)

 

2標本問題(2つの母集団の母数に関する検定)
ここでは2標本の平均の差の検定のみを扱います
確率変数Xからのm個の標本x, 確率変数Yからのn個の標本y

 

母分散が既知の場合
標本x,標本yは母分散が等しい正規分布N(μ1,σ^2),N(μ2,σ^2)に従う
標本x平均~N(μ1,σ^2/m)、標本y平均~N(μ2,σ^2/n)
d=標本x平均-標本y平均 も正規分布に従う
E[d]=μ1-μ2、V[d]=σ^2/m+σ^2/n
母分散が既知の場合には、等しくなくてもdを利用して検定可能
例)有意水準5%の棄却域
|d|≧1.96√(σ^2/m+σ^2/n)

 

母分散が未知で等しい場合
推定量はプールした分散s^2(合併した分散)を使用
t=(標本xの平均-標本yの平均)/√(s^2/m+s^2/n)
自由度m+n-2のt分布t(m+n-2)
t=(標本x平均-標本y平均)/{s*(√1/m+1/n)}
例)x若年群の片足立ち(sec)、y高齢群の片足立ち(sec)
x<-c(39,47,48,59,65,68,74,79,80,82,87,88,89,94,99,
100,101,106,107,108,111,120,122,129,130,136,138,139,
147,157,158,179,189,197,213)
y<-c(25,27,32,37,38,44,45,52,54,60,65,67,68,69,71,73,
82,88,89,92,99,100,101,116,126,131,151,156,174)

t.test(x, y, var.equal=T)

Two Sample t-test
data: x and y
t = 2.9234, df = 62, p-value = 0.004829
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
9.671654 51.500760
sample estimates:
mean of x mean of y
111.00000 80.41379

有意水準1%で差があるという結果になるのだが、
95%信頼区間は[9.67,51.50]と範囲が大きい
約9秒の差、50秒の差をどのようにとらえるかは研究者が十分に考察しなければならない

事後の分析を行ってみます
まずはt値から標本効果量を算出
2.9234*sqrt(64/(35*26))
[1] 0.7752781
次に5%水準で検定力を求めてみます
library(pwr)
pwr.t2n.test(n1=35,n2=29, d =0.7752781 , sig.level =0.05 , power = NULL,alternative = "greater")

t test power calculation

n1 = 35
n2 = 29
d = 0.7752781
sig.level = 0.05
power = 0.9205285
alternative = greater

検定力は0.92となり
検定力の高い検定であったといえます


母分散が未知で等しくない場合(ウェルチの検定)
まず等分散性を検定する
Leveneの検定
帰無仮説:xの分散=yの分散
library(car)

score<-c(100,100,101,101,106,107,108,111,116,
120,89,100,145,85,106,107,90,149,116,78)
group=factor(rep(c("x", "y"), c(10,10)))
leveneTest(score ~ group)

Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 1 6.3953 0.02101 *
18

帰無仮説は棄却され、母分散が等しくないことが判明したため
ウェルチの検定を行う

t=(標本xの平均-標本yの平均)/√( s1^2/m + s2^2/n )
このtは帰無仮説の下で近似的に次の自由度fをもつt分布に従う
g1=s1^2/m ,g2=s2^2/n
f=(g1+g2)^2/{(g1^2/(m-1))-(g2^2/(n-1))}
このtを用いる方法はウェルチの検定(Welch's test)とよばれる

ベクトルscoreよりデータを抽出する
x<-score[1:10]
y<-score[11:20]

ウェルチの補正による2標本のt検定
t.test(x,y,var=F)

Welch Two Sample t-test
data: x and y
t = 0.062759, df = 10.485, p-value = 0.9511
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-17.14089 18.14089
sample estimates:
mean of x mean of y
107.0 106.5

Rで簡単(適合度検定)

例1)次のデータはメンデルの法則 1:2:1 の確率に適合しているか?
AA Aa aa n
52 102 46 200
帰無仮説:標本はメンデルの法則に適合している

観測度数<-c(52,102,46)
確率<-c(1/4,1/2,1/4)      
chisq.test(観測度数, p=確率)

Chi-squared test for given probabilities
data: 観測度数
X-squared = 0.44, df = 2, p-value = 0.8025

結果:帰無仮説は棄却されず、標本はメンデルの法則によく適合していることが確認された。(完全に適合しているという意味ではありません)

 

例2)男性は女性よりも多いと言えるか
男性 女性
229  179
帰無仮説:男性と女性の数は同じである(1/2, 1/2に適合している)

観測度数<-c(229,189)
確率<-c(1/2,1/2)      
chisq.test(観測度数, p=確率)

Chi-squared test for given probabilities

data: 観測度数
X-squared = 3.8278, df = 1, p-value = 0.05041

結果:有意水準5%では帰無仮説は棄却されず、男女数に差があるとは言えない。
(「男女比が1:1に適合している」という結果になります.あくまで有意水準0.05での結果です.)

Rで簡単(独立性の検定)

赤文字のみをRで実行

次のクロス表を検定してみます
問)男性と女性で、yesの解答者数とnoの解答者数の比率に違いはあるか?

   Yes N0
男性 n1 n2
女性 n3 n4

帰無仮説:男性と女性では解答者数の比率に違いがない

例)仮想データ

n1<-15
n2<-11
n3<-4
n4<-16

観測度数を行列で作成して解いてみます
(x<-matrix(c(n1,n2,n3,n4),nrow=2,ncol=2, byrow=T))

カイ二乗検定
chisq.test(x, correct=F)

Pearson's Chi-squared test

data: x
X-squared = 6.6244, df = 1, p-value = 0.01006

結果:帰無仮説は有意水準5%で棄却され、男女間での解答の比率に違いが認められた。

全セル数の20%以上で期待値が5未満の場合には次のいずれかの補正を実施
Yates(イェ-ツ)補正
chisq.test(x)
フィッシャーの正確確率検定
fisher.test(x)


カイ二乗検定の結果を代入する
RES <- chisq.test(x)

期待値
RES$expected

[,1] [,2]
[1,] 10.73913 15.26087
[2,] 8.26087 11.73913

標準化残差(残差÷√期待値)
RES$residuals

[,1] [,2]
[1,] 1.300211 -1.090708
[2,] -1.482468 1.243599

調整済み(標準化)残差
RES$stdres

[,1] [,2]
[1,] 2.573799 -2.573799
[2,] -2.573799 2.573799

ファイ係数
(φ<-(n1*n4-n2*n3)/sqrt(sum(x[1,])*sum(x[2,]) *sum(x[,1])*sum(x[,2])))

 0.379486

結果:男性はYES、女性はNOと解答するという正の連関を認めた。

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