×
[PR]上記の広告は3ヶ月以上新規記事投稿のないブログに表示されています。新しい記事を書く事で広告が消えます。
source("http://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1276622688") source("http://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1276531669") source("http://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1273264153") ## データ ## 参考: http://www.psychology.emory.edu/clinical/mcdowell/PSYCH560/psych560.htm 感謝。 dat <- data.frame( ps = factor(rep(c("ws", "ch", "sk"), each=3, time=3)), cnd = factor(rep(c("nt", "sp", "ms"), each=9)), val=c( c(5,4,5,4,5,3,3,5,4), c(5,6,3,4,6,3,12,10,13), c(10,12,16,11,10,12,13,12,10) )) tapply(dat[,3], list(dat[,1], dat[,2]), mean) tapply(dat[,3], list(dat[,1], dat[,2]), sd) ## 水準名とか anms <- levels(dat[,1]) bnms <- levels(dat[,2]) ## 等分散性の検定 grp <- factor(paste(dat[,1], dat[,2], sep="")) val <- dat[,3] library(car) levene.test(val, grp) ## とりあえずプロット x <- tapply(dat[,3], list(dat[,1], dat[,2]), mean) barplot(x, beside=T) dev.off() ## 全体の分散分析 lmres <- lm(val~ps*cnd, contrasts=list(ps=contr.sum, cnd=contr.sum), data=dat) library(car) Anovares <- Anova(lmres) Anovares errordf <- Anovares[4,2]; errorss <- Anovares[4,1]; mse <- errorss/errordf ## 要因の主効果の多重比較 fact.mc(dat) ## 単純主効果 seffect(dat)[1:2] ## 単純主効果の多重比較 mcse.a <- list() mcse.b <- list() ## b水準のa要因の効果 for (nm in bnms) { x <- subset(dat, dat[2]==nm)[,c(1,3)] mcse.a[[nm]] <- holm.mc(x, mse=mse, errordf) } names(mcse.a) <- paste(names(dat)[1], " effect at ", bnms, sep="") mcse.a ## a水準のb要因の効果 for (nm in anms) { x <- subset(dat, dat[1]==nm)[,c(2,3)] mcse.b[[nm]] <- holm.mc(x, mse=mse, errordf) } names(mcse.b) <- paste(names(dat)[2], " effect at ", anms, sep="") mcse.b ## 交互作用対比 ## データフレームの分割 cmba <- combn(anms,2) cmbb <- combn(bnms,2) n <- 1 nms <- vector() datlist <- list() for (i in 1: ncol(cmba)) { for (j in 1:ncol(cmbb)) { atf <- (dat[,1]==cmba[1,i]|dat[,1]==cmba[2,i]) btf <- (dat[,2]==cmbb[1,j]|dat[,2]==cmbb[2,j]) nms[n] <- paste(paste(cmba[1,i], "-", cmba[2,i], sep=""), " vs. ", paste(cmbb[1,j], "-", cmbb[2,j], sep=""), sep="") datlist[[n]] <- subset(dat, atf&btf); n <- n+1 } } names(datlist) <- nms ### 交互作用対比の分析 iclist <- lapply(datlist, function(x) Anova(lm(val~ps*cnd, x), error=lmres)) tmp <- lapply(iclist, function(x) unlist(x[3,])) tmp2 <- as.matrix(t(sapply(tmp, rbind))) ps <- tmp2[,4] ps.holm <- p.adjust(ps, "holm") icres <- cbind(tmp2[,1:2], errorss, errordf, mse, tmp2[,3:4], sqrt(tmp2[,3]), ps.holm) colnames(icres) <- c("ss", "nmdf", "errorss", "errordf", "mse", "fv", "pv", "tv", "Pr(>F).holm") class(icres) <- "anova" icres # こんな交互作用対比のやり方もありえんわなー
PR
dat <- data.frame(
a1b1 = rnorm(10),
a1b2 = rnorm(10),
a1b3 = rnorm(10, 100, 10),
a1b4 = rnorm(10, 10, 1),
a2b1 = rnorm(10),
a2b2 = rnorm(10),
a2b3 = rnorm(10, 100, 10),
a2b4 = rnorm(10, 10, 1),
a3b1 = rnorm(10),
a3b2 = rnorm(10),
a3b3 = rnorm(10, 100, 10),
a4b4 = rnorm(10, 10, 1)
)
dat
nlva <- 3 # 自分で入力
nlvb <- 4 # 自分で入力
anms <- c("a1", "a2", "a3") # 自分で入力
bnms <- c("b1", "b2", "b3", "b4") # 自分で入力
ncl <- ncol(dat)
clnm <- 1:ncl
clmt <- matrix(clnm, nlvb)
adatlist <- list()
bdatlist <- list()
for (i in 1:nlva) {
x <- dat[clmt[,i]] ; names(x) <- bnms
adatlist[[i]] <- x
}
names(adatlist) <- anms
for (i in 1:nlvb) {
x <- dat[clmt[i,]] ; names(x) <- anms
bdatlist[[i]] <- x
}
names(bdatlist) <- bnms
adatlist
bdatlist
吉田寿夫・森 敏昭 (1990). 心理学のためのデータ解析テクニカルブック 北大路書房
# 1要因対応のなし。p. 87
dat <- data.frame(
a = factor(rep(c("a1","a2","a3","a4"),each=8)),
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)
)
dat
# 1要因対応のあり。p. 92
dat <- data.frame(
a1=c(9,7,8,8,12,11,8,13),
a2=c(6,5,6,3,6,7,10,9),
a3=c(10,13,8,13,12,14,14,16),
a4=c(9,11,13,14,16,12,15,14)
)
dat
# 2要因参加者間。p95
dat <- data.frame(
a=factor(rep(c("a1", "a2"), each=20)),
b=factor(rep(c("b1","b2", "b3","b4"), each=5, time=2)),
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))
dat
# 2要因混合p. 107
dat <- data.frame(
a=factor(rep(c("a1", "a2"), each=5)),
b1=c(3, 3, 1, 3, 5, 3, 5, 2, 4, 6),
b2=c(4, 3, 4, 5, 7, 2, 6, 3, 6, 4),
b3=c(6, 6, 6, 4, 8, 3, 2, 3, 6, 5),
b4=c(5, 7, 8, 7, 9, 2, 3, 3, 4, 6)
)
dat
# 2要因対応あり。p117
dat <- data.frame(
a1b1=c(3, 3, 1, 3, 5),
a1b2=c(4, 3, 4, 5, 7),
a1b3=c(6, 6, 6, 4, 8),
a1b4=c(5, 7, 8, 7, 9),
a2b1=c(3, 5, 2, 4, 6),
a2b2=c(2, 6, 3, 6, 4),
a2b3=c(3, 2, 3, 6, 5),
a2b4=c(2, 3, 3, 4, 6)
)
dat
# 3要因対応なし p121
dat <- data.frame(
a = rep(c("a1", "a2"), each=24),
b = rep(c("b1", "b2", "b1", "b2"), each=12),
c = rep(c("c1", "c2", "c3", "c1", "c2", "c3","c1", "c2", "c3", "c1", "c2", "c3"), each=4),
val = c(2,6,5,7, 5,7,9,9, 9,10,13,14, 6,6,8,10, 3,6,5,8, 6,7,5,6,
1,2,5,2, 2,1,4,5, 1, 5, 3, 5, 5,3,4, 6, 5,6,6,7, 5,5,9,7)
)
dat
# 3要因混合 (1要因対応あり) p137
dat <- data.frame(
a = rep(c("a1", "a2"), each=8),
b = rep(c("b1", "b2", "b1", "b2"), each=4),
c1 = c(2,6,5,7,6,6,8,10,1,2,5,2,5,3,4,6),
c2 = c(5,7,9,9,3,6,5,8,2,1,4,5,5,6,6,7),
c3 = c(9,10,13,14,6,7,5,6,1,5,3,5,5,5,9,7)
)
# 3要因混合 (2要因対応あり) 。p144
dat <- data.frame(
a=factor(c(rep("a1", 4), rep("a2", 4))),
b1c1=c(2,6,5,7,1,2,5,2), b1c2=c(5,7,9,9,2,1,4,5), b1c3=c(9,10,13,14,1,5,3,5),
b2c1=c(6,6,8,10,5,3,4,6), b2c2=c(3,6,5,8,5,6,6,7), b2c3=c(6,7,5,6,5,5,9,7)
)
dat
# 3要因被験者内 p152
dat <- data.frame(
a1b1c1 = c(2,6,5,7),
a1b1c2 = c(5,7,9,9),
a1b1c3 = c(9,10,13,14),
a1b2c1 = c(6,6,8,10),
a1b2c2 = c(3,6,5,8),
a1b2c3 = c(6,7,5,6),
a2b1c1 = c(1,2,5,2),
a2b1c2 = c(2,1,4,5),
a2b1c3 = c(1,5,3,5),
a2b2c1 = c(5,3,4,6),
a2b2c2 = c(5,6,6,7),
a2b2c2 = c(5,5,9,7)
)
##南風原朝和 (2002). 心理統計学の基礎――統合的理解のために―― 有斐閣
## p. 9 表1-1、逸脱行動と変化量データ。
sex <- factor(rep(c("male", "female"), each=20))
sixth <- c(4,9,14,16,15,14,7,13,18,20,12,6,10,16,12,16,13,12,12,15,3,9,4,8,8,5,5,7,12,11,8,11,13,10,15,10,12,10,8,9)
eighth <- c(8,12,11,20,16,18,7,19,20,14,14,12,16,14,13,17,19,10,19,14,12,9,10,13,13,10,7,18,17,14,12,18,17,20,13,16,13,12,12,13)
cng <- c(4,3,-3,4,1,4,0,6,2,-6,2,6,6,-2,1,1,6,-2,7,-1,9,0,6,5,5,5,2,11,5,3,4,7,4,10,-2,6,1,2,4,4)
dat <- data.frame(sex, sixth, eighth, cng)
dat
# p226 。表8-1。50組の母子の協調性データ。
x1 <- c(12,12,7,17,14,9,10,13,15,12,12,15,11,14,17,17,16,15,15,10,12,9,12,12,19,11,14,15,15,15,16,15,12,10,11,12,15,13,15,12,12,12,13,17,13,11,14,16,12,12) # 母親価値
x2 <- c(2,2,2,3,2,2,3,3,3,1,3,3,2,2,4,2,4,3,4,2,2,1,2,2,4,2,3,2,3,3,2,3,2,2,3,1,2,3,2,2,2,3,3,3,2,3,2,4,2,2) # 通園年数
y <- c(6,11,11,13,13,10,10,15,11,11,16,14,10,13,12,15,16,14,14,8,13,12,12,11,16,9,12,13,13,14,12,15,8,12,11,6,12,15,9,13,9,11,14,12,13,9,11,14,16,8) # 協調性
dat <- data.frame(x1,x2, y)
dat
# p271 表9-1、ねたみ感情の測定データ
dat <- data.frame(
dmn = factor(rep(c("app", "acd", "aff"), each=15)),
att = factor(rep(rep(c("fr", "mid","hst"), each=5), 3)),
envy = c(2,4,3,2,1,7,4,6,5,8,4,3,2,5,2,2,2,1,2,2,3,4,2,1,2,1,3,3,2,1,3,1,1,0,1,4,7,5,6,4,3,2,1,2,1)
)
dat
# 表9-6 ブロック別の妬み感情データ
blc <- factor(1:15)
app <- c(1,2,2,2,2,3,3,4,4,4,5,5,6,7,8)
acd <- c(1,1,1,1,2,2,2,2,2,2,2,3,3,3,4)
aff <- c(0,1,1,1,1,1,2,2,3,3,4,4,5,6,7)
dat <- data.frame(blc,app,acd,aff)
# 縦長
blc <- rep(blc, 3)
fct <- factor(rep(c("app", "acd", "aff"),each=15))
val <- c(app, acd, aff)
dat <- data.frame(blc, fct, val)
# p312. 表9-8、共分散分析のためのデータ
dmn <- factor(rep(c("app", "acd", "aff"), each=15)) # 優越属性
cv <- c(4,9,6,2,1,14,10,13,11,15,8,7,3,12,5,6,10,2,8,5,13,15,9,4,7,1,12,14,11,3,9,4,2,1,6,12,15,13,14,11,10,7,5,8,3) # 共変量
envy <- c(2,4,3,2,1,7,4,6,5,8,4,3,2,5,2,2,2,1,2,2,3,4,2,1,2,1,3,3,2,1,3,1,1,0,1,4,7,5,6,4,3,2,1,2,1) # ねたみ感情
dat <- data.frame(dmn, cv, envy)
## http://www.psychology.emory.edu/clinical/mcdowell/PSYCH560/factorw.htm
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)
)
dat
# p318 相関行列
mat <- matrix(c(
1.000,0.033,0.315,0.456,0.266,0.607,0.228,0.419,
0.033,1.000,0.637,0.250,0.528,0.195,0.522,0.420,
0.315,0.637,1.000,0.333,0.880,0.237,0.750,0.328,
0.456,0.250,0.333,1.000,0.362,0.432,0.398,0.449,
0.266,0.528,0.880,0.362,1.000,0.252,0.738,0.269,
0.607,0.195,0.237,0.432,0.252,1.000,0.335,0.463,
0.228,0.522,0.750,0.398,0.738,0.335,1.000,0.238,
0.419,0.420,0.328,0.449,0.269,0.463,0.238,1.000),
nrow=8, dimnames=list(c("onw", "yok", "gai", "sin", "sha", "kyo", "sek", "sun"), c("onw", "yok", "gai", "sin", "sha", "kyo", "sek", "sun"))
)
# 表10-8 観測変数間の相関係数 p354
mat <- matrix(c(
1.000,0.160,0.302,0.461,0.299,0.152,0.134,0.182,0.251,0.372,0.157,0.203,
0.160,1.000,0.341,0.400,0.404,0.320,0.403,0.374,0.285,0.100,0.291,-0.014,
0.302,0.341,1.000,0.372,0.552,0.476,0.467,0.572,0.316,0.408,0.393,0.369,
0.461,0.400,0.372,1.000,0.302,0.225,0.256,0.255,0.164,0.236,0.229,0.224,
0.299,0.404,0.552,0.302,1.000,0.708,0.623,0.776,0.361,0.294,0.472,0.342,
0.152,0.320,0.476,0.225,0.708,1.000,0.324,0.769,0.295,0.206,0.351,0.202,
0.134,0.403,0.467,0.256,0.623,0.324,1.000,0.724,0.260,0.071,0.204,0.152,
0.182,0.374,0.572,0.255,0.776,0.769,0.724,1.000,0.284,0.142,0.320,0.189,
0.251,0.285,0.316,0.164,0.361,0.295,0.260,0.284,1.000,0.295,0.290,0.418,
0.372,0.100,0.408,0.236,0.294,0.206,0.071,0.142,0.295,1.000,0.468,0.351,
0.157,0.291,0.393,0.229,0.472,0.351,0.204,0.320,0.290,0.468,1.000,0.385,
0.203,-0.014,0.369,0.224,0.342,0.202,0.152,0.189,0.418,0.351,0.385,1.000),
nr=12,
dimnames = list(paste("y", 1:12, sep=""), paste("y", 1:12, sep=""))
)
mat
## 媒介分析。http://davidakenny.net/dtt/mediate.htm
library(foreign)
dat <- data.frame(read.spss("http://davidakenny.net/dtt/morse_et_al.sav"))
dat
## 信頼区間グラフ
## Masson, M. E. J., & Loftus, G. R. (2003). Using confidence intervals for graphically based data interpretation. Canadian Journal of Experimental Psychology, 57, 203-220.
# http://web.uvic.ca/psyc/masson/ML.pdf
dat <- data.frame(
incong=c(784,853,622,954,634,751,918,894),
cong=c(632,702,598,873,600,729,877,801),
neutral=c(651,689,606,855,595,740,893,822))
dat
# 1要因対応のなし。p. 87
dat <- data.frame(
a = factor(rep(c("a1","a2","a3","a4"),each=8)),
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)
)
dat
# 1要因対応のあり。p. 92
dat <- data.frame(
a1=c(9,7,8,8,12,11,8,13),
a2=c(6,5,6,3,6,7,10,9),
a3=c(10,13,8,13,12,14,14,16),
a4=c(9,11,13,14,16,12,15,14)
)
dat
# 2要因参加者間。p95
dat <- data.frame(
a=factor(rep(c("a1", "a2"), each=20)),
b=factor(rep(c("b1","b2", "b3","b4"), each=5, time=2)),
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))
dat
# 2要因混合p. 107
dat <- data.frame(
a=factor(rep(c("a1", "a2"), each=5)),
b1=c(3, 3, 1, 3, 5, 3, 5, 2, 4, 6),
b2=c(4, 3, 4, 5, 7, 2, 6, 3, 6, 4),
b3=c(6, 6, 6, 4, 8, 3, 2, 3, 6, 5),
b4=c(5, 7, 8, 7, 9, 2, 3, 3, 4, 6)
)
dat
# 2要因対応あり。p117
dat <- data.frame(
a1b1=c(3, 3, 1, 3, 5),
a1b2=c(4, 3, 4, 5, 7),
a1b3=c(6, 6, 6, 4, 8),
a1b4=c(5, 7, 8, 7, 9),
a2b1=c(3, 5, 2, 4, 6),
a2b2=c(2, 6, 3, 6, 4),
a2b3=c(3, 2, 3, 6, 5),
a2b4=c(2, 3, 3, 4, 6)
)
dat
# 3要因対応なし p121
dat <- data.frame(
a = rep(c("a1", "a2"), each=24),
b = rep(c("b1", "b2", "b1", "b2"), each=12),
c = rep(c("c1", "c2", "c3", "c1", "c2", "c3","c1", "c2", "c3", "c1", "c2", "c3"), each=4),
val = c(2,6,5,7, 5,7,9,9, 9,10,13,14, 6,6,8,10, 3,6,5,8, 6,7,5,6,
1,2,5,2, 2,1,4,5, 1, 5, 3, 5, 5,3,4, 6, 5,6,6,7, 5,5,9,7)
)
dat
# 3要因混合 (1要因対応あり) p137
dat <- data.frame(
a = rep(c("a1", "a2"), each=8),
b = rep(c("b1", "b2", "b1", "b2"), each=4),
c1 = c(2,6,5,7,6,6,8,10,1,2,5,2,5,3,4,6),
c2 = c(5,7,9,9,3,6,5,8,2,1,4,5,5,6,6,7),
c3 = c(9,10,13,14,6,7,5,6,1,5,3,5,5,5,9,7)
)
# 3要因混合 (2要因対応あり) 。p144
dat <- data.frame(
a=factor(c(rep("a1", 4), rep("a2", 4))),
b1c1=c(2,6,5,7,1,2,5,2), b1c2=c(5,7,9,9,2,1,4,5), b1c3=c(9,10,13,14,1,5,3,5),
b2c1=c(6,6,8,10,5,3,4,6), b2c2=c(3,6,5,8,5,6,6,7), b2c3=c(6,7,5,6,5,5,9,7)
)
dat
# 3要因被験者内 p152
dat <- data.frame(
a1b1c1 = c(2,6,5,7),
a1b1c2 = c(5,7,9,9),
a1b1c3 = c(9,10,13,14),
a1b2c1 = c(6,6,8,10),
a1b2c2 = c(3,6,5,8),
a1b2c3 = c(6,7,5,6),
a2b1c1 = c(1,2,5,2),
a2b1c2 = c(2,1,4,5),
a2b1c3 = c(1,5,3,5),
a2b2c1 = c(5,3,4,6),
a2b2c2 = c(5,6,6,7),
a2b2c2 = c(5,5,9,7)
)
##南風原朝和 (2002). 心理統計学の基礎――統合的理解のために―― 有斐閣
## p. 9 表1-1、逸脱行動と変化量データ。
sex <- factor(rep(c("male", "female"), each=20))
sixth <- c(4,9,14,16,15,14,7,13,18,20,12,6,10,16,12,16,13,12,12,15,3,9,4,8,8,5,5,7,12,11,8,11,13,10,15,10,12,10,8,9)
eighth <- c(8,12,11,20,16,18,7,19,20,14,14,12,16,14,13,17,19,10,19,14,12,9,10,13,13,10,7,18,17,14,12,18,17,20,13,16,13,12,12,13)
cng <- c(4,3,-3,4,1,4,0,6,2,-6,2,6,6,-2,1,1,6,-2,7,-1,9,0,6,5,5,5,2,11,5,3,4,7,4,10,-2,6,1,2,4,4)
dat <- data.frame(sex, sixth, eighth, cng)
dat
# p226 。表8-1。50組の母子の協調性データ。
x1 <- c(12,12,7,17,14,9,10,13,15,12,12,15,11,14,17,17,16,15,15,10,12,9,12,12,19,11,14,15,15,15,16,15,12,10,11,12,15,13,15,12,12,12,13,17,13,11,14,16,12,12) # 母親価値
x2 <- c(2,2,2,3,2,2,3,3,3,1,3,3,2,2,4,2,4,3,4,2,2,1,2,2,4,2,3,2,3,3,2,3,2,2,3,1,2,3,2,2,2,3,3,3,2,3,2,4,2,2) # 通園年数
y <- c(6,11,11,13,13,10,10,15,11,11,16,14,10,13,12,15,16,14,14,8,13,12,12,11,16,9,12,13,13,14,12,15,8,12,11,6,12,15,9,13,9,11,14,12,13,9,11,14,16,8) # 協調性
dat <- data.frame(x1,x2, y)
dat
# p271 表9-1、ねたみ感情の測定データ
dat <- data.frame(
dmn = factor(rep(c("app", "acd", "aff"), each=15)),
att = factor(rep(rep(c("fr", "mid","hst"), each=5), 3)),
envy = c(2,4,3,2,1,7,4,6,5,8,4,3,2,5,2,2,2,1,2,2,3,4,2,1,2,1,3,3,2,1,3,1,1,0,1,4,7,5,6,4,3,2,1,2,1)
)
dat
# 表9-6 ブロック別の妬み感情データ
blc <- factor(1:15)
app <- c(1,2,2,2,2,3,3,4,4,4,5,5,6,7,8)
acd <- c(1,1,1,1,2,2,2,2,2,2,2,3,3,3,4)
aff <- c(0,1,1,1,1,1,2,2,3,3,4,4,5,6,7)
dat <- data.frame(blc,app,acd,aff)
# 縦長
blc <- rep(blc, 3)
fct <- factor(rep(c("app", "acd", "aff"),each=15))
val <- c(app, acd, aff)
dat <- data.frame(blc, fct, val)
# p312. 表9-8、共分散分析のためのデータ
dmn <- factor(rep(c("app", "acd", "aff"), each=15)) # 優越属性
cv <- c(4,9,6,2,1,14,10,13,11,15,8,7,3,12,5,6,10,2,8,5,13,15,9,4,7,1,12,14,11,3,9,4,2,1,6,12,15,13,14,11,10,7,5,8,3) # 共変量
envy <- c(2,4,3,2,1,7,4,6,5,8,4,3,2,5,2,2,2,1,2,2,3,4,2,1,2,1,3,3,2,1,3,1,1,0,1,4,7,5,6,4,3,2,1,2,1) # ねたみ感情
dat <- data.frame(dmn, cv, envy)
## http://www.psychology.emory.edu/clinical/mcdowell/PSYCH560/factorw.htm
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)
)
dat
# p318 相関行列
mat <- matrix(c(
1.000,0.033,0.315,0.456,0.266,0.607,0.228,0.419,
0.033,1.000,0.637,0.250,0.528,0.195,0.522,0.420,
0.315,0.637,1.000,0.333,0.880,0.237,0.750,0.328,
0.456,0.250,0.333,1.000,0.362,0.432,0.398,0.449,
0.266,0.528,0.880,0.362,1.000,0.252,0.738,0.269,
0.607,0.195,0.237,0.432,0.252,1.000,0.335,0.463,
0.228,0.522,0.750,0.398,0.738,0.335,1.000,0.238,
0.419,0.420,0.328,0.449,0.269,0.463,0.238,1.000),
nrow=8, dimnames=list(c("onw", "yok", "gai", "sin", "sha", "kyo", "sek", "sun"), c("onw", "yok", "gai", "sin", "sha", "kyo", "sek", "sun"))
)
# 表10-8 観測変数間の相関係数 p354
mat <- matrix(c(
1.000,0.160,0.302,0.461,0.299,0.152,0.134,0.182,0.251,0.372,0.157,0.203,
0.160,1.000,0.341,0.400,0.404,0.320,0.403,0.374,0.285,0.100,0.291,-0.014,
0.302,0.341,1.000,0.372,0.552,0.476,0.467,0.572,0.316,0.408,0.393,0.369,
0.461,0.400,0.372,1.000,0.302,0.225,0.256,0.255,0.164,0.236,0.229,0.224,
0.299,0.404,0.552,0.302,1.000,0.708,0.623,0.776,0.361,0.294,0.472,0.342,
0.152,0.320,0.476,0.225,0.708,1.000,0.324,0.769,0.295,0.206,0.351,0.202,
0.134,0.403,0.467,0.256,0.623,0.324,1.000,0.724,0.260,0.071,0.204,0.152,
0.182,0.374,0.572,0.255,0.776,0.769,0.724,1.000,0.284,0.142,0.320,0.189,
0.251,0.285,0.316,0.164,0.361,0.295,0.260,0.284,1.000,0.295,0.290,0.418,
0.372,0.100,0.408,0.236,0.294,0.206,0.071,0.142,0.295,1.000,0.468,0.351,
0.157,0.291,0.393,0.229,0.472,0.351,0.204,0.320,0.290,0.468,1.000,0.385,
0.203,-0.014,0.369,0.224,0.342,0.202,0.152,0.189,0.418,0.351,0.385,1.000),
nr=12,
dimnames = list(paste("y", 1:12, sep=""), paste("y", 1:12, sep=""))
)
mat
## 媒介分析。http://davidakenny.net/dtt/mediate.htm
library(foreign)
dat <- data.frame(read.spss("http://davidakenny.net/dtt/morse_et_al.sav"))
dat
## 信頼区間グラフ
## Masson, M. E. J., & Loftus, G. R. (2003). Using confidence intervals for graphically based data interpretation. Canadian Journal of Experimental Psychology, 57, 203-220.
# http://web.uvic.ca/psyc/masson/ML.pdf
dat <- data.frame(
incong=c(784,853,622,954,634,751,918,894),
cong=c(632,702,598,873,600,729,877,801),
neutral=c(651,689,606,855,595,740,893,822))
dat
有意水準調整型の多重比較の関数 holm.mc
source("http://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1276531669")
## 使い方の覚書
# 与えるデータは以下のように1列目に因子変数、2列目に数値変数。横長のデータはdatw=TRUEと指定する。横長の場合は変数名が因子名で、データは数値のみ
# デフォルトでプールされた標準偏差で多重比較する
# mse=で数値を入れれば、これを元にプールされた標準偏差を使って検定する
# 調整はデフォルトでholm。p.adjust.methodにあるものは選べる
# paired=Tのときは個々にt検定を繰り返すだけ
# psd=F, paired=Fで個々の対応なし検定繰り返し
# デフォルトでWelchの検定。var.equal=Tにすると普通のやつ
dat <- data.frame(a = factor(c(rep("a1",8), rep("a2",8), rep("a3",8), rep("a4",8))), 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))
dat # 縦長
dat2 <- data.frame(a1=c(9,7,8,8,12,11,8,13), a2=c(6,5,6,3,6,7,10,9), a3=c(10,13,8,13,12,14,14,16), a4=c(9,11,13,14,16,12,15,14))
dat2 #横長
res <- holm.mc(dat) # デフォルト
print.anova(res)
res1 <- holm.mc(dat, mse=5) # 適当にmseを指定する。内部ではsqrt(mse) を計算する。mseを指定した場合はその自由度をmse.df=で指定しないと警告が出る
res1.2 <- holm.mc(dat, mse=5, mse.df=30)
print.anova(res1.2)
res2 <- holm.mc(dat, psd=F, paired=T)
print.anova(res2)
res3 <- holm.mc(dat, psd=F, paired=F)
print.anova(res3)
res4 <- holm.mc(dat, psd=F, paired=F, var.equal=T)
print.anova(res4)
res5 <- holm.mc(dat2, datw=T, paired=T, psd=F) # 横長、対応あり、プールしない。対応ありのときは対象となる2群でプールするので全体では基本的にプールしない (たぶん)
print.anova(res5)
## うーむ、どうも引数の指定がオシャレじゃないな…
source("http://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1276531669")
## 使い方の覚書
# 与えるデータは以下のように1列目に因子変数、2列目に数値変数。横長のデータはdatw=TRUEと指定する。横長の場合は変数名が因子名で、データは数値のみ
# デフォルトでプールされた標準偏差で多重比較する
# mse=で数値を入れれば、これを元にプールされた標準偏差を使って検定する
# 調整はデフォルトでholm。p.adjust.methodにあるものは選べる
# paired=Tのときは個々にt検定を繰り返すだけ
# psd=F, paired=Fで個々の対応なし検定繰り返し
# デフォルトでWelchの検定。var.equal=Tにすると普通のやつ
dat <- data.frame(a = factor(c(rep("a1",8), rep("a2",8), rep("a3",8), rep("a4",8))), 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))
dat # 縦長
dat2 <- data.frame(a1=c(9,7,8,8,12,11,8,13), a2=c(6,5,6,3,6,7,10,9), a3=c(10,13,8,13,12,14,14,16), a4=c(9,11,13,14,16,12,15,14))
dat2 #横長
res <- holm.mc(dat) # デフォルト
print.anova(res)
res1 <- holm.mc(dat, mse=5) # 適当にmseを指定する。内部ではsqrt(mse) を計算する。mseを指定した場合はその自由度をmse.df=で指定しないと警告が出る
res1.2 <- holm.mc(dat, mse=5, mse.df=30)
print.anova(res1.2)
res2 <- holm.mc(dat, psd=F, paired=T)
print.anova(res2)
res3 <- holm.mc(dat, psd=F, paired=F)
print.anova(res3)
res4 <- holm.mc(dat, psd=F, paired=F, var.equal=T)
print.anova(res4)
res5 <- holm.mc(dat2, datw=T, paired=T, psd=F) # 横長、対応あり、プールしない。対応ありのときは対象となる2群でプールするので全体では基本的にプールしない (たぶん)
print.anova(res5)
## うーむ、どうも引数の指定がオシャレじゃないな…
Holm, Scheffe, TukeyHSD, Dunnet あたりがメジャー。あと、ペリの方法が検定力が高いらしい。
基本的に、群間の独立性を仮定している。すなわち、対応なしの検定である。
ただし、Holm等、有意水準調整型の場合は対応ありでも一応使える
dat <- data.frame(a = factor(c(rep("a1",8), rep("a2",8), rep("a3",8), rep("a4",8))), 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))
# 平均と合計
aggregate(dat[2], list(dat[,1]), mean)
aggregate(dat[2], list(dat[,1]), sum)
# lmで分析
lmres <- lm(result~a, dat)
library(car)
Anova(lmres)
## 色々な多重比較
# Holm法、プールされていない標準偏差
## こっちと同じ
source("http://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1275792703")
val <- dat$result
idv <- dat$a
cbn <- combn(levels(idv), 2)
tresdat <- tobj(t.test(rnorm(10), rnorm(10), paired=F)); tresdat[1,] <- NA # 結果格納用データフレーム
for (i in 1:ncol(cbn)) {
tres <- t.test(val[which(idv==cbn[1,i])], val[which(idv==cbn[2,i])], paired=F)
tob <- tobj(tres)
tresdat[i,] <- tob
}
rownames(tresdat) <- apply(cbn, 2, function(x) paste(x, collapse=","))
tresdat
## Holm法でp値の調整をする
pholms <- p.adjust(tresdat[,"p.vl"], "holm")
(tresdat <- data.frame(tresdat, pholms))
colnames(tresdat) <- c(colnames(tresdat[-ncol(tresdat)]), "Pr(>|t|).Holm")
print.anova(tresdat) ## 有意の星をつける。文字型の変数は01にされる
# Holm法、プールされた標準偏差
## pairwise.t.test
source("http://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1275849015")
pairwise.t.test(dat$result, dat$a, p.adj="holm")
ptest(dat$result, dat$a)
## 青木先生の関数。感謝
source("http://aoki2.si.gunma-u.ac.jp/R/src/Bonferroni.R", encoding="euc-jp")
Bonferroni(dat$result, dat$a, method="Holm")
print.anova(Bonferroni(dat$result, dat$a, method="Holm")$result2)
x <- Bonferroni(dat$result, dat$a, method="Holm")$result2[,2]
round(p.adjust(x, "holm"), 5)
# scheffe
## 青木先生のコードより。感謝
source("http://aoki2.si.gunma-u.ac.jp/R/src/scheffe.R", encoding="euc-jp")
ns <- tapply(dat$result, dat$a, length)
ms <- tapply(dat$result, dat$a, mean)
us <- tapply(dat$result, dat$a, var)
scheffe(ns, ms, us, 1,2)
scheffe(ns, ms, us, 1,3)
scheffe(ns, ms, us, 1,4)
scheffe(ns, ms, us, 2,3)
scheffe(ns, ms, us, 2,4)
scheffe(ns, ms, us, 3,4)
## こういう関数も見つけた。
http://www.biw.kuleuven.be/vakken/statisticsbyr/ANOVAbyRr/multiplecompJIMRC.htm
source("http://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1276112608")
aovres <- aov(result~a, dat)
summary(aovres)
nis <- tapply(dat$result, dat$a, length); nis
ms <- tapply(dat$result, dat$a, mean); ms
mse.df <- aovres$df.residual; mse.df
mse <- sum(aovres$residuals^2)/mse.df; mse
scheffeCI(ms, nis, mse.df, mse, conf=.95)
## ちょっと改造してF値とp値を出すようにした
scres <- scheffeCI2(ms, nis, mse.df, mse, conf=.95)
print.anova(scres)
# Tukey
library(multcomp)
summary(glht(lmres,linfct=mcp(a="Tukey"))) # これはなんかおかしい
summary(glht(aov(lmres),linfct=mcp(a="Tukey"))) # 教科書と合うのはこっち
TukeyHSD(aov(lmres))
## 小塩先生のspssでの分析例をやってみた。感謝
cond <- factor(rep(1:3, each=7))
result <- c(4,1,3,2,2,4,3,6,8,5,9,8,7,7,4,3,4,6,5,5,5)
aovres <- aov(result~cond)
TukeyHSD(aovres)
library(multcomp)
summary(glht(aovres, linfct=mcp(cond="Tukey")))
## TukeyHSDとglhtで微妙に結果が違うのはなんでだろーなー
## Dunnet 対照群は一番上の水準 ("a1")
summary(glht(aov(lmres),linfct=mcp(a= "Dunnett")))
## Williams なんだかよくわからない。水準間に順序が仮定できるときに使うらしい。
summary(glht(aov(lmres),linfct=mcp(a= "Williams")))
## そのうちやろう
# 2要因以上の分散分析で交互作用が出たときはプールされた分散をつかい統計量を算出する
# glhtの使い方。特にTukey
# 各群のサンプルサイズが違う場合の多重比較。
参考
http://www.ibaraki-kodomo.com/toukei/posthoc.html
http://www.gen-info.osaka-u.ac.jp/testdocs/tomocom/tazyu.html
http://home.hiroshima-u.ac.jp/keiroh/maeda/statsarekore/posthoc.html
基本的に、群間の独立性を仮定している。すなわち、対応なしの検定である。
ただし、Holm等、有意水準調整型の場合は対応ありでも一応使える
dat <- data.frame(a = factor(c(rep("a1",8), rep("a2",8), rep("a3",8), rep("a4",8))), 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))
# 平均と合計
aggregate(dat[2], list(dat[,1]), mean)
aggregate(dat[2], list(dat[,1]), sum)
# lmで分析
lmres <- lm(result~a, dat)
library(car)
Anova(lmres)
## 色々な多重比較
# Holm法、プールされていない標準偏差
## こっちと同じ
source("http://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1275792703")
val <- dat$result
idv <- dat$a
cbn <- combn(levels(idv), 2)
tresdat <- tobj(t.test(rnorm(10), rnorm(10), paired=F)); tresdat[1,] <- NA # 結果格納用データフレーム
for (i in 1:ncol(cbn)) {
tres <- t.test(val[which(idv==cbn[1,i])], val[which(idv==cbn[2,i])], paired=F)
tob <- tobj(tres)
tresdat[i,] <- tob
}
rownames(tresdat) <- apply(cbn, 2, function(x) paste(x, collapse=","))
tresdat
## Holm法でp値の調整をする
pholms <- p.adjust(tresdat[,"p.vl"], "holm")
(tresdat <- data.frame(tresdat, pholms))
colnames(tresdat) <- c(colnames(tresdat[-ncol(tresdat)]), "Pr(>|t|).Holm")
print.anova(tresdat) ## 有意の星をつける。文字型の変数は01にされる
# Holm法、プールされた標準偏差
## pairwise.t.test
source("http://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1275849015")
pairwise.t.test(dat$result, dat$a, p.adj="holm")
ptest(dat$result, dat$a)
## 青木先生の関数。感謝
source("http://aoki2.si.gunma-u.ac.jp/R/src/Bonferroni.R", encoding="euc-jp")
Bonferroni(dat$result, dat$a, method="Holm")
print.anova(Bonferroni(dat$result, dat$a, method="Holm")$result2)
x <- Bonferroni(dat$result, dat$a, method="Holm")$result2[,2]
round(p.adjust(x, "holm"), 5)
# scheffe
## 青木先生のコードより。感謝
source("http://aoki2.si.gunma-u.ac.jp/R/src/scheffe.R", encoding="euc-jp")
ns <- tapply(dat$result, dat$a, length)
ms <- tapply(dat$result, dat$a, mean)
us <- tapply(dat$result, dat$a, var)
scheffe(ns, ms, us, 1,2)
scheffe(ns, ms, us, 1,3)
scheffe(ns, ms, us, 1,4)
scheffe(ns, ms, us, 2,3)
scheffe(ns, ms, us, 2,4)
scheffe(ns, ms, us, 3,4)
## こういう関数も見つけた。
http://www.biw.kuleuven.be/vakken/statisticsbyr/ANOVAbyRr/multiplecompJIMRC.htm
source("http://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1276112608")
aovres <- aov(result~a, dat)
summary(aovres)
nis <- tapply(dat$result, dat$a, length); nis
ms <- tapply(dat$result, dat$a, mean); ms
mse.df <- aovres$df.residual; mse.df
mse <- sum(aovres$residuals^2)/mse.df; mse
scheffeCI(ms, nis, mse.df, mse, conf=.95)
## ちょっと改造してF値とp値を出すようにした
scres <- scheffeCI2(ms, nis, mse.df, mse, conf=.95)
print.anova(scres)
# Tukey
library(multcomp)
summary(glht(lmres,linfct=mcp(a="Tukey"))) # これはなんかおかしい
summary(glht(aov(lmres),linfct=mcp(a="Tukey"))) # 教科書と合うのはこっち
TukeyHSD(aov(lmres))
## 小塩先生のspssでの分析例をやってみた。感謝
cond <- factor(rep(1:3, each=7))
result <- c(4,1,3,2,2,4,3,6,8,5,9,8,7,7,4,3,4,6,5,5,5)
aovres <- aov(result~cond)
TukeyHSD(aovres)
library(multcomp)
summary(glht(aovres, linfct=mcp(cond="Tukey")))
## TukeyHSDとglhtで微妙に結果が違うのはなんでだろーなー
## Dunnet 対照群は一番上の水準 ("a1")
summary(glht(aov(lmres),linfct=mcp(a= "Dunnett")))
## Williams なんだかよくわからない。水準間に順序が仮定できるときに使うらしい。
summary(glht(aov(lmres),linfct=mcp(a= "Williams")))
## そのうちやろう
# 2要因以上の分散分析で交互作用が出たときはプールされた分散をつかい統計量を算出する
# glhtの使い方。特にTukey
# 各群のサンプルサイズが違う場合の多重比較。
参考
http://www.ibaraki-kodomo.com/toukei/posthoc.html
http://www.gen-info.osaka-u.ac.jp/testdocs/tomocom/tazyu.html
http://home.hiroshima-u.ac.jp/keiroh/maeda/statsarekore/posthoc.html
