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

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

統計学備忘録 since2016

Rを使用した統計学

二元配置分散分析 分散分析表の理解

データは仮想です

繰り返し数がすべて等しい場合 (=5)

   

因子B

 

 

B1

B2

B3

因子A

A1

7

9

21

11

11

22

8

10

19

10

9

18

9

11

22

A2

23

12

9

21

16

8

26

16

11

27

17

10

22

15

9

エクセルにペーストして、青文字のみをコピー


Data<-c(x$B1,x$B2,x$B3)x<-read.table("clipboard",header=T)


以下赤文字のみRで実行

inshiA1<-c(rep("A1",5),rep("A2",5))

inshiA<-c(rep(inshiA1,3))

inshiB<-c(rep("B1",10),rep("B2",10),rep("B3",10))

summary ( aov ( Data ~inshiA* inshiB) )

              Df Sum Sq Mean Sq F value   Pr(>F)   

inshiA         1   67.5    67.5   21.89 9.40e-05 ***

inshiB         2   73.3    36.6   11.88 0.000259 ***

inshiA:inshiB  2  850.2   425.1  137.87 6.94e-14 ***

Residuals     24   74.0     3.1                    

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

分散分析表の計算方法を忘れないように記載しておきます!

 

Df

SumSq

MeanSq

Fvalue

inshiA

1

67.5

67.5

21.89

inshiB

2

73.3

36.6

11.88

inshiA:inshiB

2

850.2

425.1

137.87

Residuals

24

74

3.1

 

 

Df:自由度、SumSq:平方和、MeanSq:平均平方(平方和÷自由度)、

Fvalue:F値(…の平均平方÷水準内変動の平均平方)

(平均平方、F値は平方和から算出可能.よって、平方和を調べるべし!)

 

http://yoshida931.hatenablog.com/entry/2017/02/19/144239

 

inshiA:因子Aの水準間の変動を調べる

inshiB:因子Bの水準間の変動を調べる

inshiA:inshiB:因子Aと因子Bの交互作用の変動を調べる

Residuals:水準内の変動を調べる

(変動:各測定値から平均値へのずれ(バラつき)が、「0」にならないように平均値からの変位の二乗和をとったもの.分散分析では平方和.)

 

Df(自由度)

inshiA:a-1

inshiB:b-1

inshiA:inshiB:(a-1)×(b-1)

Residuals:a×b×(n-1)

 

因子Aの水準数a、因子Bの水準数b、繰り返し数n、データ総数a*b*n

総数= length(data)、総平均=mean(data)とします
a=2
b=3
n=5

 

因子Aの水準A1・A2間の変動(平方和)

inshiASumSq=67.5         
A1<-c(Data[1:5],Data[11:15],Data[21:25]) #(A1,B1),(A1,B2),(A1,B3)全てのデータ
mean(A1)
A2<-c(Data[6:10],Data[16:20],Data[26:30])  #(A2,B1),(A2,B2),(A2,B3)全てのデータ
mean(A2)
n*b*( (mean(A1)-mean(Data))^2+(mean(A2)-mean(Data))^2)

 

f:id:yoshida931:20170310132759p:plain

 

因子Bの水準B1・B2・B3間の変動(平方和)
inshiB SumSq73.3
B1<-c(Data[1:10])  #(A1,B1),(A2,B1)全てのデータ
mean(B1)
B2<-c(Data[11:20])  #(A1,B2),(A2,B2)全てのデータ
mean(B2)
B3<-c(Data[21:30])  #(A1,B3),(A2,B3)全てのデータ
mean(B3)
n*a*( (mean(B1)-mean(Data))^2+(mean(B2)-mean(Data))^2+(mean(B3)-mean(Data))^2)

 

交互作用の変動
inshiA:inshiBのSumSq=850.2
水準内平均‐水準A平均-水準B平均+総平均
n*(mean(Data[1:5])-mean(A1)-mean(B1)+mean(Data))^2+
 n*(mean(Data[6:10])-mean(A2)-mean(B1)+mean(Data))^2+
 n*(mean(Data[11:15])-mean(A1)-mean(B2)+mean(Data))^2+
 n*(mean(Data[16:20])-mean(A2)-mean(B2)+mean(Data))^2+
 n*(mean(Data[21:25])-mean(A1)-mean(B3)+mean(Data))^2+
 n*(mean(Data[26:30])-mean(A2)-mean(B3)+mean(Data))^2

f:id:yoshida931:20170310142941p:plain

 

各水準の平均値をプロットしたグラフ
x<- data.frame(A1=c(9,10,20.4),
A2=c(23.8,15.2,9.4))
plot(x$A1,type="o",lty=3,xaxt="n",xlim = c(0.8,3.2),ylim=c(7,25),ylab="",xlab = "",bty = "l")
lines(x$A2,pch=19,type="o")
axis(side=1,at=1:3,labels=c("B1","B2","B3"))
legend(2.7,25, legend = colnames(x) , pch = pchs, lty = ltys)

f:id:yoshida931:20170310172329p:plain

 

水準内の変動
ResidualsSumSq74
sum( (Data[1:5]-mean(Data[1:5]))^2)+
 sum( (Data[6:10]-mean(Data[6:10]))^2)+
 sum( (Data[11:15]-mean(Data[11:15]))^2)+
 sum( (Data[16:20]-mean(Data[16:20]))^2)+
 sum( (Data[21:25]-mean(Data[21:25]))^2)+
 sum( (Data[26:30]-mean(Data[26:30]))^2)

f:id:yoshida931:20170310140533p:plain


注意)水準内変動はいわゆる各水準内のバラつきです.
F値は水準内平均平方を分母にとるため、バラつきが
小さいと交互作用や因子内の水準間差が出やすくなります

例)
交互作用の変動:F=交互作用の変動の平均平方 ÷ 水準内変動の平均平方

 

相関係数の区間推定

サンプル:日本統計学会 (編集); 日本統計学会公式認定 統計検定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

 

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

なぜt検定を繰り返してはいけないのか?

sample は「石村卓夫;分散分析のはなし,東京図書,1992,p137」

A<-c(12,14,16)
B<-c(13,15,17)
C<-c(15,17,19)
D<-c(17,19,21)

x<-c(rep("A",3),rep("B",3),rep("C",3),rep("D",3))
y<-c(A,B,C,D)
stripchart(y~x,vertical=T)

f:id:yoshida931:20170224165734p:plain

視覚的にはAD間、BD間には差があるように見えます

 

重要な前提です
Aの母平均をμa
Bの母平均をμb
Cの母平均をμc
Dの母平均をμd とします
「μa=μb=μc=μd」を検定する場合にt検定を繰り返してはいけないのです.
ある特別なペアの違いなどを確認したい場合にはt検定を使用します.研究内容により異なるので間違いないように.

 

参考 http://www.gen-info.osaka-u.ac.jp/MEPHAS/ave.html
2種類の既存薬AとBを組み合わせた配合薬Cの配合効果を評価する場合.この場合、既存薬A,Bのそれぞれの効果と配合薬Cの効果を比較します.ここで いいたいのはCが既存薬A、Bの両方よりも効果があるということです.「CがAよ りも優れている、かつCがBよりも優れている」ということを 示します.つまり帰無仮説はC=AかつC=Bとなり、A=B=Cとは異なります.このような場合には2標本t検定を繰り返して用いるほうが適切だと考えられます.

 

一元配置分散分析
oneway.test(y~x,var.equal = T)
#F = 3.6875, num df = 3, denom df = 8, p-value = 0.06214
#有意水準5%で差が「ない」という結果になりました

 

#全ての群の母分散は等しいと仮定してt検定を繰り返してみます
AB<-t.test(A,B,paired=FALSE, var.equal=T, conf.level=0.95)
AC<-t.test(A,C,paired=FALSE, var.equal=T, conf.level=0.95)
AD<-t.test(A,D,paired=FALSE, var.equal=T, conf.level=0.95)
BC<-t.test(B,C,paired=FALSE, var.equal=T, conf.level=0.95)
BD<-t.test(B,D,paired=FALSE, var.equal=T, conf.level=0.95)
CD<-t.test(C,D,paired=FALSE, var.equal=T, conf.level=0.95)

 

#それぞれの検定結果のP値のみを取り出してみます
str(AB3);str(AC3);str(AD3);str(BC3);str(BD3);str(CD3)
AB 0.573
AC 0.14
AD 0.0376
BC 0.288
BD 0.0705
CD 0.288

 やはりAD間には有意水準5%で「差がある」という結果になりました.しかし分散分析では差がないという結果でした.

 

分散分析の帰無仮説
有意水準0.5において「μa=μb=μc=μd」と仮定する

t検定の帰無仮説
有意水準0.5において「μa=μb」、
かつ有意水準0.5において「μa =μc」、
かつ有意水準0.5において「μa =μd」、
かつ有意水準0.5において「μb=μc」、
かつ有意水準0.5において「μb=μd」、
かつ有意水準0.5において「μc¬=μd」と仮定した検定になります.

したがってt検定を使った「μa=μb=μc=μd」の確率は(1-0.05)^6となります.少なくとも一組には差がある確率は1-(1-0.05)^6=0.2649081となり、有意水準26.5%の検定になってしまいます。つまりαエラーの危険性が高くなるのです.このような矛盾を解消するための手法が多重比較です(しかし、本当にt検定は繰り返して使ってはいけないのだろうか? 本当に有意水準の高い検定になっているのだろうか? 2群間比較にはt検定しかないのだが…疑問が残ります)

多重比較(Tukey)を行ってみます

TukeyHSD(aov(y ~ x),ordered = TRUE)
Tukey multiple comparisons of means
95% family-wise confidence level
factor levels have been ordered

Fit: aov(formula = y ~ x)

$x
  diff  lwr    upr   p adj
B-A 1 -4.2294189 6.229419 0.9252929
C-A 3 -2.2294189 8.229419 0.3245304
D-A 5 -0.2294189 10.229419 0.0609499
C-B 2 -3.2294189 7.229419 0.6297636
D-B 4 -1.2294189 9.229419 0.1441838
D-C 2 -3.2294189 7.229419 0.6297636

有意水準0.05ではどのペアにも差はありませんでした.

等分散性の検定

準備としてカイ二乗分布と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更新)

赤文字の部分だけ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