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

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

統計学備忘録 since2016

理学療法、医療、福祉、小児

一元配置分散分析

平方和の分解から分散分析を理解する

よくみかける一元配置のデータを分解してみます

 

f:id:yoshida931:20161102130649j:plain

 

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

例)繰り返し数が5の場合
治療後の満足度調査
20名を無作為に分けて治療A,B,C,D後にVASで調査した
治療=要因、A,B,C,D=水準(群)

A B C D
10 10 5 9
6 5 4 5
10 5 11 2
9 12 4 3
10 4 6 1

(注:上記一覧と行列が逆になります。また少々数値がずれてますが、Rに落とせば問題ありません。)

 

自分でも何やってるのか分からなくなるのでエクセルイメージで追跡していきます

こんな感じです

f:id:yoshida931:20161104191605p:plain

(次からはセルのみのコピペです)
行数が繰り返し数
列数が水準数


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

 

#各群の平均(各水準の平均)

(群平均<-colMeans(x))
...................................
A B C D
9.0 7.2 6.0 4.0
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

 f:id:yoshida931:20161123211908p:plain

 

データのみを抜き出してベクトルにします

全データ<-c(x$A,x$B,x$C,x$D)

[1] 10 6 10 9 10 10 5 5 12 4 5 4 11 4 6 9 5 2 3 1

 

水準をベクトルに変換します

もちろん 水準<-c("A","A","A" ,"A", "A", "B" ,"B" ,"B" ,"B" ,"B" ,"C" ,"C" ,"C" ,"C" ,"C" ,"D" ,"D" ,"D", "D", "D") と入力しても構いませんが、次のようにrep関数を使用します。

水準<-c(rep("A",5),rep("B",5),rep("C",5),rep("D",5))

 

さらに要因型のベクトルに変換します

(水準2<-factor(水準)) 

この時点で水準が数値に変換される…str(水準2)で確認

.............................................................................

[1] A A A A A B B B B B C C C C C D D D D D
Levels: A B C D
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

 

rbind(全データ,水準2)

#rbind() 行方向にまとめる

[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
全データ 10 6 10 9 10 10 5 5 12 4 5 4
水準2 1 1 1 1 1 2 2 2 2 2 3 3
[,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20]
全データ 11 4 6 9 5 2 3 1
水準2 3 3 3 4 4 4 4 4

(エクセルのイメージ)

f:id:yoshida931:20161124124959p:plain

 

#群平均行列
(群平均行列<-matrix(rep(群平均,5),nrow=5,ncol=4,byrow=TRUE))
.....................................
[,1] [,2] [,3] [,4]
[1,] 9 7.2 6 4
[2,] 9 7.2 6 4
[3,] 9 7.2 6 4
[4,] 9 7.2 6 4
[5,] 9 7.2 6 4
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

f:id:yoshida931:20161123211847p:plain

 

#全データの平均
mean(全データ)
[1] 6.55

 

#(全データ)ー(全データの平均
x-mean(全データ)
...............................
A B C D
1 3.45 3.45 -1.55 2.45
2 -0.55 -1.55 -2.55 -1.55
3 3.45 -1.55 4.45 -4.55
4 2.45 5.45 -2.55 -3.55
5 3.45 -2.55 -0.55 -5.55
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

 f:id:yoshida931:20161123211818p:plain

 

#(群平均行列)ー(全データの平均
群平均行列-mean(全データ)
..................................
[,1] [,2] [,3] [,4]
[1,] 2.45 0.65 -0.55 -2.55
[2,] 2.45 0.65 -0.55 -2.55
[3,] 2.45 0.65 -0.55 -2.55
[4,] 2.45 0.65 -0.55 -2.55
[5,] 2.45 0.65 -0.55 -2.55
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

f:id:yoshida931:20161123212336p:plain

 

#(全データ)ー(群平均)
x-群平均行列
...........................
A B C D
1 1 2.8 -1 5
2 -3 -2.2 -2 1
3 1 -2.2 5 -2
4 0 4.8 -2 -1
5 1 -3.2 0 -3
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

f:id:yoshida931:20161123212620p:plain

 

 

#(全データー全平均)^2  この合計が総合平方和
#各データが全平均に近いと小さく、離れると大きくなる
(x-mean(全データ))^2
..........................................
A B C D
[1,] 11.9025 11.9025 2.4025 6.0025
[2,] 0.3025 2.4025 6.5025 2.4025
[3,] 11.9025 2.4025 19.8025 20.7025
[4,] 6.0025 29.7025 6.5025 12.6025
[5,] 11.9025 6.5025 0.3025 30.8025
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

f:id:yoshida931:20161123212724p:plain 

 

 

#群間^2 この合計が群間平方和
#各群の平均が全平均に近いと小さく、離れると大きくなる
(群平均行列-mean(全データ))^2
.........................................
[,1] [,2] [,3] [,4]
[1,] 6.0025 0.4225 0.3025 6.5025
[2,] 6.0025 0.4225 0.3025 6.5025
[3,] 6.0025 0.4225 0.3025 6.5025
[4,] 6.0025 0.4225 0.3025 6.5025
[5,] 6.0025 0.4225 0.3025 6.5025
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

f:id:yoshida931:20161123212856p:plain

 

#郡内^2 この合計が群内平方和
#各データが群平均に近いと小さく、離れると大きくなる
(x-群平均行列)^2
.......................
A B C D
[1,] 1 7.84 1 25
[2,] 9 4.84 4 1
[3,] 1 4.84 25 4
[4,] 0 23.04 4 1
[5,] 1 10.24 0 9
^^^^^^^^^^^^^^^^^^^^^^^^^^

f:id:yoshida931:20161123213016p:plain

 

#総平方和
sum( (x-mean(全データ))^2)
[1] 202.95

 

#群間平方和、水準間平方和、A間平方和
sum( (群平均行列-mean(全データ))^2)
[1] 66.15

 

#群内平方和、残差平方和

#「群の効果に応じて全員の得点が一律に押し上げられたり、
引き下げられたした」というだけでは説明できないデータのばらつき
sum( (x-群平均行列)^2)
[1] 136.8

 

#自由度の考え方(ncol列数、nrow行数)
#群間の自由度は因子数(列の数)ー1
ncol(x)-1
#群内の自由度は(各群のデータ数ー1)× 因子数
(nrow(x)-1)*ncol(x)
#全体の自由度は全データ数—1
length(全データ)-1

 

#平均平方
群間の平均平方
(sum( (群平均行列-mean(全データ))^2))/(ncol(x)-1 )

[1] 22.05

#群内の平均平方
sum( (x-群平均行列)^2)/( (nrow(x)-1 )*ncol(x))

[1] 8.55

#全体の平均平方(=全体の分散)
(sum( (x-mean(まとめ))^2))/(length(まとめ)-1)

var(まとめ) #上式と同じ

 

#F値 = 群間の平均平方 ÷ 群内の平均平方
(F値=(sum( (群平均行列-mean(全データ))^2))/(ncol(x)-1 ))/
(sum( (x-群平均行列)^2)/( (nrow(x)-1)*ncol(x)))

 [1] 2.578947

 

#p値
pf(F値,3,16,lower.tail=FALSE)  

#R関数の解と違うので、直接F値を挿入する
pf(2.578947,3,16,lower.tail=FALSE)

[1] 0.08978798

 

#結果

F>2.579でp=0.0898となり治療間で満足度の差はないと言えます

 

 

実際に検定する場合には、以下の関数だけでOKです
summary(aov(全データ~水準2))

...............................................................................

                            Df    Sum Sq   Mean Sq   F value   Pr(>F)
ABCDまとめ        3    66.15       22.05       2.579       0.0898 .
Residuals           16   136.80      8.55
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

 

参考)山田 剛史, 杉澤 武俊, 村井 潤一郎; Rによるやさしい統計学, オーム社 , 2008