myscratchpad

Rでの分散分析関数はaov, anovaが用意されている。 被験者内要因があっても誤差項を指定できるので便利。 記述例 (松井先生のサイトを参照)
# 2要因混合釣合型 　aov(result ~ a * b + Error(s:a + s:a:b), SPFp.q)))

ただし、どちらもType I平方和を使っている。 Type II, IIIを使うにはcarパッケージのAnova関数を使う。 ただし、こちらは誤差項は指定できない。 あと、線形回帰なのでlm関数も同じ結果を返す (summary()に渡すオブジェクトは違う)

平方和の解説
中澤先生: http://phi.ypu.jp/swtips/Rnewsarchives.html#SAS-SS
岸本先生: http://mat.isc.chubu.ac.jp/fpr/fpr1999/0051.html

carパッケージのAnova関数を使う。データはiris
library(car)

irissp <- c(1, 2, 3, 4)
irissp <- as.factor(irissp)
idat <- data.frame(irissp)
lmres <- lm(cbind(iris$Sepal.Length, iris$Sepal.Width, iris$Petal.Length, iris$Petal.Width) ~ 1, contrasts=contr.sum)
Aovi <- Anova(lmres, idata=idat, idesign=~irissp, type=3)
 
anovakunの結果と一致
datkun<- iris[1:4]
anovakun(datkun, "sA", 4)

Rで分散分析を行うにはaov関数、もしくはanova関数を用いる…
というのが普通だと思うが、井関先生の作成された関数 "anovakun" を使うのが一番よいと思う。
詳細は上記リンクより

以下、使い方の備忘録
## 入力例
anovakun(df, "AsB", 2, 3, holm=T)
# sの左 (A) が被験者間、sの右 (B) が被験者内
# 2要因被験者間計画は"ABs", 2要因混合は"AsB"
# A、Bの順に水準数を記述。上の例は被験者間要因Aが2水準、被験者内要因Bが3水準

## なぜかわからないが、水準数で並べ替えておかないとうまく計算できなかったことがある

carパッケージのAnova関数はlm, glmのオブジェクトなんかも渡せると昨日知った。
しかし、反復測定を含む混合要因の分散分析について解説したものが無かったので (ヘルプにはあるけど) 色々ネットを彷徨って試してみた。
比較的わかりやすい例があったのでまねしてみた。
また、ヘルプも参照

しかし、やっぱり割りと面倒なようだ。
被験者内要因の変数"名"だけ代入したデータフレームをつくり、
それからlm関数のオブジェクトを渡して混合にするようだ。

sex(2), type(3), 以上被験者間、time(3) 被験者内の
2x3x3の混合3要因分散分析を行う。
##データ生成
set.seed(123) # 乱数の種
# 被験者内要因の水準1
purc<-rnorm(100, mean=50, sd=10)
# 被験者内要因の水準2
oneyr<-purc-20+rnorm(100)
# 被験者内要因の水準3
twoyr<-oneyr-10+rnorm(100)
# 被験者間要因1。2水準
sex<-rep(c(0,1),c(50,50))
#被験者間要因2。3水準
type<-rep(c(0,1,2,0,1,2),c(17,17,16,17,17,16))
#数値操作
oneyr[sex==0 & type==1]<-oneyr[sex==0 & type==1]-10 + rnorm(17)
twoyr[sex==0 & type==1]<-twoyr[sex==0 & type==1]-10 + rnorm(17)

##分析準備
#被験者間要因をfactorにしておく
sex<-as.factor(sex)
type<-as.factor(type)
#carパッケージの読み込み
library(car)

###こっから重要
## 被験者内要因の水準数だけ、factor変数をつくり、それをidatというデータフレームにする
time<-c(0,1,2) # 0がpurc, 1がoneyr, 2がtwoyr。下のlm関数内、cbindで結合するときと同じ順
time<-as.factor(time) # factorにする
idat<-data.frame(time) # データフレームにする
mod1<-lm(cbind(purc,oneyr,twoyr)~sex*type,contrasts=list(sex=contr.sum, type=contr.sum)) # lm関数にcbindで結合した被験者内要因、それを予測する被験者間要因を入れる。被験者間要因は*でつながないと交互作用が検定されない。
#lm関数を実行する際、contrasts=list(sex=contr.sum, type=contr.sum) を入れて対比行列を指定する。こうしないとType III平方和で分析したとき変な結果になっちゃう。
井関先生のサイトを参考にさせていただいた。感謝

##　分散分析の実行
aov1<-Anova(mod1, idata=idat, idesign=~time) # デフォルトのTypeII平方和で分散分析
summary(aov1) 
aov2 <- Anova(mod1, idata=idat, idesign=~time, type="III") # TypeIII平方和を指定して分散分析
summary(aov2)

##以上

###以上の分析をanovakunで行ってみた
datkun <-data.frame(cbind(sex, type, purc, oneyr, twoyr)) # anovakun用のデータフレーム
anovakun(datkun, "ABsC", 2,3,3, type2=T) # anovakunはTypeIIをオプションで指定
anovakun(datkun, "ABsC", 2,3,3) # デフォルトのTypeIIIで分析

Revelle, W.
psychパッケージの開発者。心理学者
http://www.personality-project.org/revelle.html
Using R for psychological research: A simple guide to an elegant package
http://www.personality-project.org/r/

Baron, J.
http://www.sas.upenn.edu/~baron/
Notes on the use of R for psychology experiments and questionnaires
http://www.psych.upenn.edu/~baron/rpsych/rpsych.html

Fox, J.
semパッケージ、carパッケージの開発者。社会学 (?)
http://socserv.mcmaster.ca/jfox/