×
[PR]上記の広告は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で分析
しかし、反復測定を含む混合要因の分散分析について解説したものが無かったので (ヘルプにはあるけど) 色々ネットを彷徨って試してみた。
比較的わかりやすい例があったのでまねしてみた。
また、ヘルプも参照
しかし、やっぱり割りと面倒なようだ。
被験者内要因の変数"名"だけ代入したデータフレームをつくり、
それから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で分析
PR
