myscratchpad

各水準間の差得点を条件間で比較するだけ
入戸野先生のホームページの情報、pdfが大変分かりやすかった。感謝

# データの準備。データ解析テクニカルブックのp. 116より
dat <- data.frame(
s=factor(c("s1", "s2", "s3", "s4", "s5")),
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
sapply(dat[2:9], sum)
sapply(dat[2:9], mean)

options(contrasts = c("contr.sum", "contr.sum"))
lmres <- lm(cbind(a1b1,a1b2,a1b3,a1b4,a2b1,a2b2,a2b3,a2b4)~1, dat)
library(car)
afact <- factor(c(rep("a1", 4), rep("a2", 4)))
bfact <- factor(c(rep(c("b1","b2","b3","b4"),2)))
idat <- data.frame(afact, bfact)
Anovares <- Anova(lmres, idata=idat, idesign=~afact*bfact, type=3)
summary(Anovares)

# 交互作用対比の準備
data1 <- dat[2:5]
data2 <- dat[6:9]
commat <- combn(x=1:4, 2) # combn関数で組み合わせをつくる
combnN <- ncol(commat)
c1 <- commat[1,]
c2 <- commat[2,]

# t検定結果格納用のオブジェクトを用意する
tres <- data.frame(matrix(1, 6, 9))
colnames(tres) <- c("d1mean", "d2mean", "d1sd", "d2sd", "t.value", "df", "p.value", "CIlower", "CIupper")
tresam <- data.frame(matrix(1,6,2)) # 文字列用
colnames(tresam) <- c("alternative", "method")
rwnames <- letters[1:6] # 行名用
# t検定繰り返し
for (i in 1:combnN) {
d11 <- data1[c1[i]]
d12 <- data1[c2[i]]
d21 <- data2[c1[i]]
d22 <- data2[c2[i]]
d1 <- d11[,1]-d12[,1]
d2 <- d21[,1]-d22[,1]
tobject <- t.test(d1, d2, alternative="two.sided", paired=T)
tres.vct <- c(mean(d1), mean(d2), sd(d1), sd(d2), tobject[[1]], tobject[[2]], tobject[[3]], tobject[[4]][1], tobject[[4]][2])
tres[i,] <- tres.vct
tresam[i,] <- c(tobject[[7]], tobject[[8]])
cntname <- paste(names(d11), "-", names(d12), " vs. ", names(d21), "-", names(d22), sep="")
rwnames[i] <- cntname
}

rownames(tres) <- rwnames
tres <- data.frame(tres, tresam)
# t検定結果
tres
# holm法でp値の調整
(p.holm <- p.adjust(tres["p.value"][,1], "holm"))
(tresp <- data.frame(tres, p.holm=p.holm))

## 入戸野先生の解説は有意水準を=.05として、これをp値の順位で割って有意水準自体を調整している。
hprank <- c(3,5,6,1,4,2) # 差の小ささ。最も差が小さいのは3番目ということ
hp <- rep(0.05, 6)
hp/hprank # これとtres["p.value"][,1]を比較する
(tresp2 <- data.frame(tres, p.holm=p.holm, p.ntn=hp/hprank))
## なので、有意水準を.05とし、tres["p.value"][,1]を順位倍して.05を超えるかどうか見ても同じ。p.adjustでやってるのはこれ。ただし、第n+1ステップが第nステップの有意水準より小さくならないように調整する。小さくなった場合はn+1の水準を使う


## 別のやり方
dat2 <- dat[c(2,5,6,9)]
options()$contrasts

lmres2 <- lm(cbind(a1b1, a1b4, a2b1, a2b4)~1, dat2)
afact2 <- factor(c("a1", "a1", "a2", "a2"))
bfact2 <- factor(c("b1", "b4", "b1", "b4"))
idat2 <- data.frame(afact2, bfact2)
idat2
library(car)
Anovares2 <- Anova(lmres2, idata=idat2, idesign=~afact2*bfact2, type=3)
Ares2 <- Astat(Anovares2)
Ares2$unv
tres[3,]

## グラフにしてみる
pmat <- matrix(sapply(dat2, mean), 2)
dimnames(pmat) <- list(c("b1", "b4"), c("a1", "a2"))

palette(grey.colors(2))
barplot(pmat, beside=T)
legend("top", rownames(pmat), fill = 1:2)
win.graph()
barplot(t(pmat), beside=T)
legend("top", colnames(pmat), fill = 1:2)