×
[PR]上記の広告は3ヶ月以上新規記事投稿のないブログに表示されています。新しい記事を書く事で広告が消えます。
Jonathan Baronの分散分析の解説でSPSSで平方和をタイプ2と指定すると
aov関数での結果と同じになる、と書いてあり、疑問に思ったがよくわからなかったので
Rjpwiki と 青木先生の掲示板で聞いてみた。感謝。
周辺度数が比例しているとき ("要因"が直交していると) 、タイプ1とタイプ2は同じになるようだ。
Jonathan Baronのもよく読んでみると
平方和の疑問が解決したので、ずっと下書き状態だった混合3要因分散分析のコードを書いておこう。 基本はJonathan Baronのもの。少し原典のコードの誤りを修正した。
#データ生成
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)
## 見やすくした (上記HPのコードを少し修正した。HPのままこぴぺするとエラーになる)
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))
# aovで分析
aovres <- aov(rt ~ grp*color*shape + Error(subj/(color+shape)), data=Hays.df)
summary(aovres)
## HPの記述は誤差項のカッコをつけわすれている; summary(aov(rt ~ grp*color*shape + Error(subj/shape+color), data=Hays.df))
aov関数での結果と同じになる、と書いてあり、疑問に思ったがよくわからなかったので
Rjpwiki と 青木先生の掲示板で聞いてみた。感謝。
周辺度数が比例しているとき ("要因"が直交していると) 、タイプ1とタイプ2は同じになるようだ。
Jonathan Baronのもよく読んでみると
The example aov() analysis below can be compared with the results of SPSS using SSTYPE(2)
「以下のaovでの分析例」では…と言っている。私の無知である。ジョナサンッッ！すまないッッ！！
でもその前にSPSS produces the same output as R if the user tells SPSS to calculate the Type I SS (SSTYPE(1)) or Type II SS (SSTYPE(2)) instead of the default SSTYPE(3).って言ってるよね…
平方和の疑問が解決したので、ずっと下書き状態だった混合3要因分散分析のコードを書いておこう。 基本はJonathan Baronのもの。少し原典のコードの誤りを修正した。
#データ生成
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)
## 見やすくした (上記HPのコードを少し修正した。HPのままこぴぺするとエラーになる)
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))
# aovで分析
aovres <- aov(rt ~ grp*color*shape + Error(subj/(color+shape)), data=Hays.df)
summary(aovres)
## HPの記述は誤差項のカッコをつけわすれている; summary(aov(rt ~ grp*color*shape + Error(subj/shape+color), data=Hays.df))
PR
dat <- read.table("http://personality-project.org/R/datasets/R.appendix1.data", header=T)
colnames(dat) <- c("x", "y") # 列名を短くする
aovres <- aov(y ~ x, dat)
summary(aovres)
pairwise.t.test(dat$y, dat$x, p.adj = "holm")
# p値だけで自由度もt値も出ないとはどういうことだ
以下参考
http://www.personality-project.org/r/#anova
colnames(dat) <- c("x", "y") # 列名を短くする
aovres <- aov(y ~ x, dat)
summary(aovres)
pairwise.t.test(dat$y, dat$x, p.adj = "holm")
# p値だけで自由度もt値も出ないとはどういうことだ
以下参考
http://www.personality-project.org/r/#anova
最近全然やらないので書き忘れていたが、t検定をメモしておこう
対応のない場合 (独立2標本) のときと対応のあるときで記述が違うので注意
等分散を仮定するときはvar.equal=TRUEとする。しないけど
対応のない場合
t.test(y ~ x, dat)
## xが独立変数、yが従属変数。xは因子変数にしなくてもよい
対応のある場合
t.test(datp[,1], datp[,2], paired=TRUE)
## ここで、pairedを指定しないと対応なしの検定になり、同じ結果を返す
## pairedの前で片側検定を指定できる (デフォルトは両側) 。まあ両側やっときゃいいだろう
## "greater"は前者が後者より大きい (datp[,1]>[,2]) 、"less"だとその逆。y~xの形式でどう指定するかは不明。というかできない (たぶん)
## データフレームは以下のような感じ
# このブログのエディタはエクセルから貼り付けたら勝手にtableつくってくれた。便利
対応のない場合 (独立2標本) のときと対応のあるときで記述が違うので注意
等分散を仮定するときはvar.equal=TRUEとする。しないけど
対応のない場合
t.test(y ~ x, dat)
## xが独立変数、yが従属変数。xは因子変数にしなくてもよい
対応のある場合
t.test(datp[,1], datp[,2], paired=TRUE)
## ここで、pairedを指定しないと対応なしの検定になり、同じ結果を返す
## pairedの前で片側検定を指定できる (デフォルトは両側) 。まあ両側やっときゃいいだろう
## "greater"は前者が後者より大きい (datp[,1]>[,2]) 、"less"だとその逆。
通りすがりさん、ありがとうございます
t.test(y~x, dat, alternative="greater")
ところで、グループ1 > グループ2 を検定しているが、
xのカテゴリーに割り当てた数値 (e.g., 0と1) のどちらがグループ1かわからない (たぶん数値の若い方をグループ1としてるんだと思うが)
# これで調べよう
aggregate(dat[,2], list(dat[,1]), mean, na.rm=TRUE) # 平均
t.test(y~x, dat, alternative="greater")
ところで、グループ1 > グループ2 を検定しているが、
xのカテゴリーに割り当てた数値 (e.g., 0と1) のどちらがグループ1かわからない (たぶん数値の若い方をグループ1としてるんだと思うが)
# これで調べよう
aggregate(dat[,2], list(dat[,1]), mean, na.rm=TRUE) # 平均
## データフレームは以下のような感じ
| 対応なし dat |
対応あり datp |
|||
| x | y | y1 | y2 | |
| 1 | 7 | 7 | 3 | |
| 1 | 5 | 5 | 2 | |
| 1 | 5 | 5 | 1 | |
| 1 | 9 | 9 | 4 | |
| 1 | 5 | 5 | 1 | |
| 0 | 3 | |||
| 0 | 2 | |||
| 0 | 1 | |||
| 0 | 4 | |||
| 0 | 1 | |||
# このブログのエディタはエクセルから貼り付けたら勝手にtableつくってくれた。便利
Error(...) で誤差項を指定しないなら
drop1(aovmodel, ~., test="F")
でできる。
誤差項を指定すると
以下にエラー switch(mode(x), `NULL` = structure(NULL, class = "formula"), :
無効なモデル式です
というエラーメッセージ
つまり、結局これならcarパッケージのAnova関数を使っても同じ
参考URL
http://myowelt.blogspot.com/2008/05/obtaining-same-anova-results-in-r-as-in.html
drop1(aovmodel, ~., test="F")
でできる。
誤差項を指定すると
以下にエラー switch(mode(x), `NULL` = structure(NULL, class = "formula"), :
無効なモデル式です
というエラーメッセージ
つまり、結局これならcarパッケージのAnova関数を使っても同じ
参考URL
http://myowelt.blogspot.com/2008/05/obtaining-same-anova-results-in-r-as-in.html
エヴェリット (2007) を読んでのメモ書き
lme, もしくはlmer関数で実行される
変量効果 (random effect) として個体 (~1|sub) や個体と反復測定 (~time|sub) を指定できる
前者は反復測定ANOVAでランダム切片モデルといわれるもの
後者は混合ANOVAでランダム切片・傾きモデルといわれる …ような気がする
lmerの方が誤差項指定をいろいろ複雑にできるらしい
"|"は前者で後者をネストしたもの (？)
http://shiriuskun.srv7.biz/toukei_hosoku/GLM/GLM_contents.htm
lme, もしくはlmer関数で実行される
変量効果 (random effect) として個体 (~1|sub) や個体と反復測定 (~time|sub) を指定できる
前者は反復測定ANOVAでランダム切片モデルといわれるもの
後者は混合ANOVAでランダム切片・傾きモデルといわれる …ような気がする
lmerの方が誤差項指定をいろいろ複雑にできるらしい
"|"は前者で後者をネストしたもの (？)
http://shiriuskun.srv7.biz/toukei_hosoku/GLM/GLM_contents.htm
