枝分かれ分散分析　Nested ANOVA on R, SAS and SPSS - よりぬきカツラくん　～　ナンコクの生活

　
ここに書いてあることは古くて現状に即していないと思いますし、私は統計学者でもありませんので、別の資料を参照されたほうがよいと思います。誤りがあればご指摘ください。2012年6月15日
　
　後輩のデータは枝分かれnested（階層的;hierarchical）ANOVAだと思うが、アドバイスしようと下手に思ったおかげで解析に苦労する。

sample data:
V1 V2 V3
1 1 4
1 1 3
1 2 4
1 2 8
2 3 5
2 3 6
2 4 8
2 4 10

V1は処理（固定因子）、V2は個体番号（ランダム因子）、V3は反復データ。

SASならこれでいけそうだ。

proc glm;
class V1 V2;
model V3 = V1 V2(V1)/ss1; /* nested-designのときは例数不揃いでもtypeIを使うらしい*1 */
random V2(V1);　/* 一番上の独立変数もランダム因子ならrandom V1 V2(V1) */
test h=V1 e=V2(V1); /* 固定因子V1の検定は、誤差項をV2(V1)にする。*/
run;

SPSSはシンタックス・エディタから。

glm V3 by V1 V2
/METHOD = SSTYPE(1)
/INTERCEPT = INCLUDE
/CRITERIA = ALPHA(.05)
/random = V2
/design = V1, V2 within V1.

最後のピリオドに注意。一番上がランダム因子なら、もちろんV1もrandomの部分に含める（SPSSはほとんど使わないので、どうするかは皆さんで調べてください。すいません。）

V1<-gl(2,4)
V2<-gl(4,2)
V3<-c(4,3,4,8,5,6,8,10)
summary(aov(V3~V1/V2))
#V1（上位のカテゴリー）に関しては、出力のままのP値では間違い。枝分かれの場合は直下の平方和で検定しなくてはならない。
summary(aov(V3~V1+Error(V1:V2))) # significance test for V1 effect
# OR
# より単純であてはまり（尤度）がよい、バランスのとれたモデルはどっちか？（モデル選択）
#V1をランダム要因と見なす場合
library(nlme)
lme1 <- lme(V3~1, random=~1|V2)
lme2 <- lme(V3~1, random=~1|V1/V2)
anova(lme1)
anova(lme2)
anova(lme1, lme2)
# AICの小さいほうを採用する。
VarCorr(lme2) #最尤推定された分散成分
#V1を固定因子と見なす場合
library(nlme)
lme1 <- lme(V3~V1, random=~1|V1/V2,method="ML")
lme2 <- lme(V3~1, random=~1|V1/V2,method="ML")
anova(lme1, lme2)
V1が固定因子の場合は、制限付き最尤法(REML)は使えない。最尤法(ML)を使う。しかし、この方法による検定は甘くなるという(Pinheiro & Bates [isbn:0387989579], pp. 87-)。

・この場合ではV1は２つのカテゴリーしかないが、３つ以上のカテゴリーがある場合がある。その場合の事後検定（多重比較）について、青木先生のページの掲示板で聞いたところ、Underwood,A.J. 1997. Experiments in ecology.[isbn:0521556961]に参考になる記述があるという情報をいただいた（未確認）。
その内容
http://hosho.ees.hokudai.ac.jp/~kato/e_ecology/chap9/9.1_9.6.html
*2
・Nested-anovaはSokal & Rohlfの"Biometry"（初版）には詳しい解説があったが、最新版でもあるんだろうな。(⇒ある。05.7.14追記）
　
ヘボサンプル(R用） 05.7.14追記

# シミュレーション
# 父親オスのそれぞれに母メスを交配させ、それぞれの母メスが数匹の子を産むとする
# 父親・母親が子の表現型に及ぼす効果を調べる。
#
# nested構造のデータ生成
fa=10 #父親の数
mo=5 #それぞれの父親と交配させる母親の数
n1=5 #各母親から生まれる子の数
n=fa*mo*n1 # トータルのデータ数
#
father<-rep(1:fa,rep(n/fa,fa)) # 父のカテゴリー生成
mother<-rep(1:(fa*mo),rep(n1,fa*mo)) # 母のカテゴリー生成
### 母親の番号は全体の実験を通じて個体別につけなくてもよい。父親の中で個体別につけていれば同じ結果を返す。
y<-c(rnorm(n,mean=0,sd=1)) # 父母の効果がないときのデータ
#
# 父だけが子の表現型に効果を持つ場合のデータ
#
r1=rnorm(fa,mean=1,sd=1)
father.effect<-rep(r1,rep(n/fa,fa))
y.father=y+father.effect
#
# 母だけが効果を持つ場合のデータ
#
r2<-rnorm(fa*mo,mean=1,sd=1)
mother.effect<-rep(r2,rep(n1,fa*mo))
y.mother<-y+mother.effect
#
#両親が効果を持つ場合のデータ
#
y.parents<-y+mother.effect+father.effect
#
# データフレーム生成
tab<-data.frame(father,mother,y,y.father,y.mother,y.parents,father.effect,mother.effect)
tab$father<-as.factor(tab$father)
tab$mother<-as.factor(tab$mother)
#
# 解析
# (注）父親・母親の効果は直下の平均平方で検定する。そのままのP値では間違い.
#
# 両親とも効果がある場合
summary(aov(y.parents~father/mother,tab)) #このままでは父親のほうのP値は正しくない。母親のぶんは、直下の平均平方（誤差項）で検定しているのでOK。
summary(aov(y.parents~father+Error(father:mother),tab)) #これが正しい
#
# 父親だけに効果がある場合
summary(aov(y.father~father/mother,tab))
summary(aov(y.father~father+Error(father:mother),tab))
#　母親だけに効果がある場合
summary(aov(y.mother~father/mother,tab))
summary(aov(y.mother~father+Error(father:mother),tab))
#　両親とも効果がない場合
summary(aov(y~father/mother,tab))
summary(aov(y~father+Error(father:mother),tab))

*1:例えば以下参照；http://www.statsoft.com/textbook/stglm.html#s2

*2:この本に対する批判的書評http://www.nature.museum.city.fukui.fukui.jp/gakugei/iso/argo/nlindex.html#no1。この本に限らず、全ての調査に統計を厳密に適用しろという見解には、私もやや原理主義的なものを感じます。ただ、正確な分析法や計画法を知っておくのは無駄ではないし、いろいろなデザインに応用のきくANOVAは有効なツールであると思います。