myscratchpad

source("https://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1276622688")
source("https://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1276531669")
source("https://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  # こんな交互作用対比のやり方もありえんわなー