myscratchpad

Rで変量要因を投入した分析を
行うにはlme4パッケージのlmer関数を使うらしい
これでランダム切片/傾きモデルを検討できる。
あと、glmmというのもあるらしい

…そのうち調べよう

aov関数の方がいろんな用途に使えて便利そうなのだが、タイプ1平方和を計算する。
なので、carパッケージのAnova関数を勉強しようと思っていた。
タイプIIかタイプIIIかは議論の分かれるところらしく、別に個人的にはタイプIIでもよいと思う。
しかし、上述のとおりaovではタイプIしか出ない…と思っていたのだが、基本に返って (今頃) Jonathan Baronの解説を読んでみると、SASやSPSSで平方和の指定をデフォルトのタイプ3からタイプ2にすれば、RとSASやSPSSは同じ結果になるよ、と書いてある。
「RでSASやSPSSと同じ出力が欲しい！」という人はそういう回答を求めてないと思うのだが…

ということは、非釣合型の分散分析でもaovでOKじゃないか、と思ったので上記HPの解説とanovakunのタイプ2指定の出力を比較してみた。

上記HP、6.8.1 Example 1: Mixed Models Approach to Within-Subject Factors (Hays, 1988, Table 13.21.2, p. 518)のUnbalanced (Nonorthogonal) Designs セクションより
以下転載

#データ生成
data1<-c(49,47,46,47,48,47,41,46,43,47,46,45,48,46,47,45,49,44,44,45,42,45,45,40,49,46,47,45,49,45,41,43,44,46,45,40,45,43,44,45,48,46,40,45,40,45,47,40)
### 交互作用ありにしたいときは以下を使う (適当にデータをいじったもの)
###data1<-c(49,10,46,47,48,47,41,46,43,5,46,5,48,46,47,45,100,44,44,45,42,45,45,40,49,46,47,45,49,45,41,43,44,46,100,40,45,43,44,45,48,46,40,45,10,45,10,40)
## 見やすくした (上記HPのコードを少し修正した。HPのままこぴぺするとエラーになる)
matrix(data1, ncol= 4, dimnames = list(paste("subj", 1:12), c("Shape1.Color1", "Shape2.Color1" ,"Shape1.Color2", "Shape2.Color2"))) 
## anovakunでも検証するためにデータフレームにした
datkun <- data.frame(matrix(data1, ncol= 4, dimnames = list(paste("subj", 1:12), c("Shape1.Color1", "Shape2.Color1" ,"Shape1.Color2", "Shape2.Color2"))))
# aov用のデータフレーム作成
Hays.df <- data.frame(rt = data1, subj = factor(rep(paste("subj", 1:12, sep=""), 4)),　shape = factor(rep(rep(c("shape1", "shape2"), c(12, 12)), 2)),　color = factor(rep(c("color1", "color2"), c(24, 24))))
# 非釣合型のデータにするため、grp変数を付け加える
Hays.df$grp <- factor(rep(c(1,1,1,1,1,1,1,1,2,2,2,2), 4))
## anovakun用のデータフレームにも付け加える
grp <- c(1,1,1,1,1,1,1,1,2,2,2,2)
datkun <- data.frame(cbind(grp, datkun))
# aovで分析
summary(aov(rt ~ grp*color*shape + Error(subj/(color+shape)), data=Hays.df))
## HPの記述は誤差項のカッコをつけわすれている; summary(aov(rt ~ grp*color*shape + Error(subj/shape+color), data=Hays.df))
# anovakunで分析
anovakun(datkun, "AsBC", 2, 2, 2, type2=T)

出力は全く同じになった。

結論: もうaovでやろう。lm, glmを使ったとき、気が向いたらAnovaを使おう。あと、そろそろまじめに平方和について勉強しよう
最近知ったのだが、anova関数はそもそも尤度を使ってモデル比較をするための関数らしい。知らないことが多いものだ

誤差項はスラッシュを使うよりも以下の形にしたほうが見やすい出力になる
被験者*対応がある要因をs:aとする。
後は足していくだけ
s:a + s:b + s:a:b

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水準

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