×
[PR]上記の広告は3ヶ月以上新規記事投稿のないブログに表示されています。新しい記事を書く事で広告が消えます。
反復測定分散分析のチュートリアル…を集めたブログ記事より
Repeated measures ANOVA with r (tutorials)
Good Tutorials
personality-project: Rをつかった分散分析の基本。最後のほうに反復測定の例 (原文: A basic tutorial about ANOVA with R (only the last bit holds some example of repeated measures) on personality-project)
motor control lab: さらに詳しいチュートリアル (原文: A thorough tutorial on motor control lab)
UCLA seminar page: さらに詳しいチュートリアル (原文: A thorough tutorial on UCLA seminar page)
平方和のタイプに関する解説
http://prometheus.scp.rochester.edu/zlab/sites/default/files/InteractionsAndTypesOfSS.pdf
このブログに、数学/統計のグラフを踊ってみた、という大変どうでもいい動画のリンクがあったので貼っておこう。30秒あたりから
Repeated measures ANOVA with r (tutorials)
Good Tutorials
personality-project: Rをつかった分散分析の基本。最後のほうに反復測定の例 (原文: A basic tutorial about ANOVA with R (only the last bit holds some example of repeated measures) on personality-project)
motor control lab: さらに詳しいチュートリアル (原文: A thorough tutorial on motor control lab)
UCLA seminar page: さらに詳しいチュートリアル (原文: A thorough tutorial on UCLA seminar page)
平方和のタイプに関する解説
http://prometheus.scp.rochester.edu/zlab/sites/default/files/InteractionsAndTypesOfSS.pdf
このブログに、数学/統計のグラフを踊ってみた、という大変どうでもいい動画のリンクがあったので貼っておこう。30秒あたりから
PR
UCLAのページを見ていたら、aov関数の誤差項の指定は簡単にできると初めて知った
うーむ、奥が深い
http://www.ats.ucla.edu/stat/R/seminars/Repeated_Measures/repeated_measures.htm
# 松井先生のサイトの例からコピペさせていただく。感謝
# 対応がある1要因分散分析(RBデザイン) (p.91)
RB <- data.frame(
a = factor(c(rep(1,8), rep(2,8), rep(3,8), rep(4,8))),
s = factor(rep(1:8, 4)),
result = c(9,7,8,8,12,11,8,13, 6,5,6,3,6,7,10,9,
10,13,8,13,12,14,14,16, 9,11,13,14,16,12,15,14))
print(summary(aov(result ~ a + s + Error(a:s), RB)))
print(summary(aov(result ~ a + Error(s), RB))) # Error(s)だけでよい
# 1要因に対応がある2要因分散分析(SPFp.qデザイン) データ数が同じ(p.107)
SPFp.q <- data.frame(
a = factor(c(rep(1,20), rep(2,20))),
b = factor(rep(c(rep(1,5), rep(2,5), rep(3,5), rep(4,5)),2)),
s = factor(c(rep(1:5, 4), rep(6:10, 4))),
result = c(3,3,1,3,5, 4,3,4,5,7, 6,6,6,4,8, 5,7,8,7,9,
3,5,2,4,6, 2,6,3,6,4, 3,2,3,6,5, 2,3,3,4,6))
print(summary(aov(result ~ a * b + Error(s:a + s:a:b), SPFp.q)))
print(summary(aov(result ~ a * b + Error(s), SPFp.q))) # 上に同じ
# 反復測定要因が2要因以上のときは結果が異なる。松井先生のコードどおりにやらなくてはならない。うーむ、奥が深い。
Jonathan Baronの解説をいつか読もう
うーむ、奥が深い
http://www.ats.ucla.edu/stat/R/seminars/Repeated_Measures/repeated_measures.htm
# 松井先生のサイトの例からコピペさせていただく。感謝
# 対応がある1要因分散分析(RBデザイン) (p.91)
RB <- data.frame(
a = factor(c(rep(1,8), rep(2,8), rep(3,8), rep(4,8))),
s = factor(rep(1:8, 4)),
result = c(9,7,8,8,12,11,8,13, 6,5,6,3,6,7,10,9,
10,13,8,13,12,14,14,16, 9,11,13,14,16,12,15,14))
print(summary(aov(result ~ a + s + Error(a:s), RB)))
print(summary(aov(result ~ a + Error(s), RB))) # Error(s)だけでよい
# 1要因に対応がある2要因分散分析(SPFp.qデザイン) データ数が同じ(p.107)
SPFp.q <- data.frame(
a = factor(c(rep(1,20), rep(2,20))),
b = factor(rep(c(rep(1,5), rep(2,5), rep(3,5), rep(4,5)),2)),
s = factor(c(rep(1:5, 4), rep(6:10, 4))),
result = c(3,3,1,3,5, 4,3,4,5,7, 6,6,6,4,8, 5,7,8,7,9,
3,5,2,4,6, 2,6,3,6,4, 3,2,3,6,5, 2,3,3,4,6))
print(summary(aov(result ~ a * b + Error(s:a + s:a:b), SPFp.q)))
print(summary(aov(result ~ a * b + Error(s), SPFp.q))) # 上に同じ
# 反復測定要因が2要因以上のときは結果が異なる。松井先生のコードどおりにやらなくてはならない。うーむ、奥が深い。
Jonathan Baronの解説をいつか読もう
LearnR
http://www.fort.usgs.gov/brdscience/LearnR.htm
どっかでやった集中講義らしい
r-help検索。tolstoyやfinziよりも見やすい
http://www.mail-archive.com/r-help@stat.math.ethz.ch/info.html
http://tolstoy.newcastle.edu.au/R/
http://finzi.psych.upenn.edu/search.html
r-site searchのFirefox アドイン。未導入
http://addictedtor.free.fr/rsitesearch/
UCLAのRトップ
UCLA R learning module
Rのコーディングシステム
Practical Regression and Anova using R
by Julian J. Faraway
http://cran.r-project.org/doc/contrib/Faraway-PRA.pdf
Using R for Data Analysis and Graphics: Introduction, Examples and Commentary
by John Maindonald
http://wwwmaths.anu.edu.au/~johnm/r/usingR-3.pdf
http://wwwmaths.anu.edu.au/~johnm/r/
R 入門
http://brieger.esalq.usp.br/CRAN/doc/contrib/manuals-jp/R-intro-170.jp.pdf
Using SAS Proc Mixed to Fit Multilevel Models, Hierarchical Models, and Individual Growth Models という論文の例をRでやってみた
http://www.ats.ucla.edu/stat/R/paperexamples/singer/default.htm
元の論文: http://gseweb.harvard.edu/~faculty/singer/Papers/sasprocmixed.pdf
http://www.fort.usgs.gov/brdscience/LearnR.htm
どっかでやった集中講義らしい
r-help検索。tolstoyやfinziよりも見やすい
http://www.mail-archive.com/r-help@stat.math.ethz.ch/info.html
http://tolstoy.newcastle.edu.au/R/
http://finzi.psych.upenn.edu/search.html
r-site searchのFirefox アドイン。未導入
http://addictedtor.free.fr/rsitesearch/
UCLAのRトップ
UCLA R learning module
Rのコーディングシステム
Practical Regression and Anova using R
by Julian J. Faraway
http://cran.r-project.org/doc/contrib/Faraway-PRA.pdf
Using R for Data Analysis and Graphics: Introduction, Examples and Commentary
by John Maindonald
http://wwwmaths.anu.edu.au/~johnm/r/usingR-3.pdf
http://wwwmaths.anu.edu.au/~johnm/r/
R 入門
http://brieger.esalq.usp.br/CRAN/doc/contrib/manuals-jp/R-intro-170.jp.pdf
Using SAS Proc Mixed to Fit Multilevel Models, Hierarchical Models, and Individual Growth Models という論文の例をRでやってみた
http://www.ats.ucla.edu/stat/R/paperexamples/singer/default.htm
元の論文: http://gseweb.harvard.edu/~faculty/singer/Papers/sasprocmixed.pdf
ここがわかりやすかったので、覚書。というかコピペ。官舎
dv <- c(8,6,7,7,6,4,5,5,6,3,5,3,3,6,2,3,4,2,2,3,6,7,5,4,6)
grp <- gl(5,5)
aovres <- aov(dv~grp)
summary(aovres)
model.tables(aovres)
# デフォルトの対比
options()$contrasts
contrasts(grp)
# ヘルマート対比で検定
coef.hel <- contr.helmert(5)
coef.hel
contrasts(grp) <- coef.hel
aovres.h <- aov(dv~grp)
summary(aovres.h, split=list(grp=c(1,2,3,4)))
# grp: C1 1 12.100 12.100 9.0299 0.0069971 ** grp1 vs. grp2
# grp: C2 1 12.033 12.033 8.9801 0.0071289 ** grp1&2 vs. grp3
# grp: C3 1 19.267 19.267 14.3781 0.0011441 ** grp1&2&3 vs. grp4
# grp: C4 1 4.840 4.840 3.6119 0.0718780 . grp1&2&3&4 vs. grp5
# 多項式対比で検定
coef.ply <- contr.poly(5)
coef.ply
contrasts(grp) <- coef.ply
aovres.p <- aov(dv~grp)
summary(aovres.p, split=list(grp=c(1,2,3,4)))
# 自分で対比を設定する
coef.s <- cbind(c(1,-1,1,0,-1), c(-1,-1,-1,4,-1))
coef.s
sum(coef.s)
contrasts(grp) <- coef.s
aovres.s <- aov(dv~grp)
summary(aovres.s, intercept=T, split=list(grp=c(1,2)))
# grp: C1 1 0.20 0.20 0.1493 0.7033289 # grp1&3 vs. grp2&5
# grp: C2 1 23.04 23.04 17.1940 0.0004994 *** # grp4 vs. grp1&2&3&5
dv <- c(8,6,7,7,6,4,5,5,6,3,5,3,3,6,2,3,4,2,2,3,6,7,5,4,6)
grp <- gl(5,5)
aovres <- aov(dv~grp)
summary(aovres)
model.tables(aovres)
# デフォルトの対比
options()$contrasts
contrasts(grp)
# ヘルマート対比で検定
coef.hel <- contr.helmert(5)
coef.hel
contrasts(grp) <- coef.hel
aovres.h <- aov(dv~grp)
summary(aovres.h, split=list(grp=c(1,2,3,4)))
# grp: C1 1 12.100 12.100 9.0299 0.0069971 ** grp1 vs. grp2
# grp: C2 1 12.033 12.033 8.9801 0.0071289 ** grp1&2 vs. grp3
# grp: C3 1 19.267 19.267 14.3781 0.0011441 ** grp1&2&3 vs. grp4
# grp: C4 1 4.840 4.840 3.6119 0.0718780 . grp1&2&3&4 vs. grp5
# 多項式対比で検定
coef.ply <- contr.poly(5)
coef.ply
contrasts(grp) <- coef.ply
aovres.p <- aov(dv~grp)
summary(aovres.p, split=list(grp=c(1,2,3,4)))
# 自分で対比を設定する
coef.s <- cbind(c(1,-1,1,0,-1), c(-1,-1,-1,4,-1))
coef.s
sum(coef.s)
contrasts(grp) <- coef.s
aovres.s <- aov(dv~grp)
summary(aovres.s, intercept=T, split=list(grp=c(1,2)))
# grp: C1 1 0.20 0.20 0.1493 0.7033289 # grp1&3 vs. grp2&5
# grp: C2 1 23.04 23.04 17.1940 0.0004994 *** # grp4 vs. grp1&2&3&5
dat <- data.frame(
ps=paste("p", 1:27, sep=""),
trt=gl(3, 9, labels=c("a1", "a2", "a3")), # "notreat", "spaced", "massed"
grp=rep(gl(3, 3, labels=c("b1", "b2", "b3")),3), #"washers", "checkers", "seekers"
value=c(5,4,5,4,5,3,3,5,4,5,6,3,4,6,3,12,10,13,10,12,16,11,10,12,13,12,10)
)
agdat <- aggregate(dat[4], list(dat[,2], dat[,3]), mean)
lmres <- lm(value~trt*grp, dat)
library(car)
Anovares <- Anova(lmres)
Anovares
smry <- summary(lmres)
mse <- smry$sigma^2 # 残差の平均平方和
by(dat, dat$trt, function(x) Anova(lm(value~grp, data=x), error=lmres))
by(dat, dat$grp, function(x) Anova(lm(value~trt, data=x), error=lmres))
vnm <- paste(agdat[,1], agdat[,2], sep="")
for (i in 1:9) assign(vnm[i], agdat[,3][i])
cv <- c(1,-1,1,-1) # 対比係数
n <- 3 # 各群の人数
nm <- (sum(cv*c(a1b1, a2b1, a1b3, a2b3)))^2 # (a1b1-a2b1+a1b3-a2b3) b1水準のa1, a2の差とb3水準のa1, a2の差の差をとって二乗しただけ
dn <- (sum((cv^2)*mse))/n
(fv <- nm/dn)
(tv <- sqrt(fv))
1-pf(fv, df1=1, df2=18) # 分子の自由度は1。分母の自由度は分散分析の誤差自由度
# F(1, 18) = 18.03, p < .001. b1群とb3群はa1 -> a2への反応に違いがあるということ (原文: The significant F indicates that Washers (b1) and Seekers (b3) respond differently at the two levels of Factor A (a1 and a2) .
(1-pt(tv, 18))*2
参考:
http://www.psychology.emory.edu/clinical/mcdowell/PSYCH560/factorw.htm
http://core.ecu.edu/psyc/wuenschk/docs30/Factorial-Computations.doc
http://www.sprweb.org/articles/Keselman98.pdf
ps=paste("p", 1:27, sep=""),
trt=gl(3, 9, labels=c("a1", "a2", "a3")), # "notreat", "spaced", "massed"
grp=rep(gl(3, 3, labels=c("b1", "b2", "b3")),3), #"washers", "checkers", "seekers"
value=c(5,4,5,4,5,3,3,5,4,5,6,3,4,6,3,12,10,13,10,12,16,11,10,12,13,12,10)
)
agdat <- aggregate(dat[4], list(dat[,2], dat[,3]), mean)
lmres <- lm(value~trt*grp, dat)
library(car)
Anovares <- Anova(lmres)
Anovares
smry <- summary(lmres)
mse <- smry$sigma^2 # 残差の平均平方和
by(dat, dat$trt, function(x) Anova(lm(value~grp, data=x), error=lmres))
by(dat, dat$grp, function(x) Anova(lm(value~trt, data=x), error=lmres))
vnm <- paste(agdat[,1], agdat[,2], sep="")
for (i in 1:9) assign(vnm[i], agdat[,3][i])
cv <- c(1,-1,1,-1) # 対比係数
n <- 3 # 各群の人数
nm <- (sum(cv*c(a1b1, a2b1, a1b3, a2b3)))^2 # (a1b1-a2b1+a1b3-a2b3) b1水準のa1, a2の差とb3水準のa1, a2の差の差をとって二乗しただけ
dn <- (sum((cv^2)*mse))/n
(fv <- nm/dn)
(tv <- sqrt(fv))
1-pf(fv, df1=1, df2=18) # 分子の自由度は1。分母の自由度は分散分析の誤差自由度
# F(1, 18) = 18.03, p < .001. b1群とb3群はa1 -> a2への反応に違いがあるということ (原文: The significant F indicates that Washers (b1) and Seekers (b3) respond differently at the two levels of Factor A (a1 and a2) .
(1-pt(tv, 18))*2
参考:
http://www.psychology.emory.edu/clinical/mcdowell/PSYCH560/factorw.htm
http://core.ecu.edu/psyc/wuenschk/docs30/Factorial-Computations.doc
http://www.sprweb.org/articles/Keselman98.pdf
