myscratchpad

source("https://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1276531669")
source("https://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1273264153")


## データはテクニカルブック 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
(mmat <- matrix(mean(dat), nr=4, dimnames=list(bnms, anms)))
(sdmat <- matrix(sd(dat), nr=4, dimnames=list(bnms, anms)))

## とりあえずプロット
barplot(mmat, beside=T)
dev.off()

## 水準とかの情報
lvnms <- c("a1b1", "a1b2", "a1b3", "a1b4", "a2b1", "a2b2", "a2b3", "a2b4")
ncl <- ncol(dat)
nlva <- 2 # 自分で入力
nlvb <- 4 # 自分で入力
anms <- c("a1", "a2") # 自分で入力
bnms <- c("b1", "b2", "b3", "b4") # 自分で入力
## データフレームの分割
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


## 全体の分散分析
lmres <- lm(as.matrix(dat[lvnms])~1, dat)
library(car)
afact <- factor(rep(anms, each=4))
bfact <- factor(rep(bnms, 2))
idat <- data.frame(afact, bfact)
Anovares <- Anova(lmres, idata=idat, idesign=~afact*bfact)
Ares <- Astat(Anovares)
Aunv <- Ares$unv
Aunv


## 要因の主効果の多重比較: 水準別誤差項
# a要因
tmp <- list()
for (nm in anms) {
  tmp[[nm]] <- rowMeans(adatlist[[nm]])
}
ameans <- data.frame(tmp)
factmc.a <- holm.mc(ameans, datw=T, psd=F, paired=T)
factmc.a

# b要因
tmp <- list()
for (nm in bnms) {
  tmp[[nm]] <- rowMeans(bdatlist[[nm]])
}
bmeans <- data.frame(tmp)
factmc.b <- holm.mc(bmeans, datw=T, psd=F, paired=T)
factmc.b

## 単純主効果 ##
## 水準別誤差項で単純主効果 ##
## bにおけるaの効果とその多重比較
afct <- factor(anms); iadat <- data.frame(afct)
seffect.alist <- list(); seffect.a <- list()
seffect.alist <- lapply(bdatlist, function(x) Anova(lm(as.matrix(x)~1, x), idata=iadat, idesign=~afct))
mcse.a <- list()
mcse.a <- lapply(bdatlist, function(x) holm.mc(x, datw=T, psd=F, paired=T))
seffect.a <- lapply(seffect.alist, function(x) Astat(x)$unv)
seffect.a # 単純主効果
mcse.a # 多重比較

## aにおけるbの効果とその多重比較
bfct <- factor(bnms); ibdat <- data.frame(bfct)
seffect.blist <- list(); seffect.b <- list()
seffect.blist <- lapply(adatlist, function(x) Anova(lm(as.matrix(x)~1, x), idata=ibdat, idesign=~bfct))
mcse.b <- list()
mcse.b <- lapply(adatlist, function(x) holm.mc(x, datw=T, psd=F, paired=T))
seffect.b <- lapply(seffect.blist, function(x) Astat(x)$unv)
seffect.b # 単純主効果
mcse.b # 多重比較


## 交互作用対比
### 処理用オブジェクト
icreslist <- list(); n <- 1
lisnm <- vector()
cmba <- combn(anms, 2)
cmbb <- combn(bnms, 2)
### 交互作用対比の分析
for (i in 1:ncol(cmba)) {
  x1 <- adatlist[[cmba[1,i]]]
  x2 <- adatlist[[cmba[2,i]]]
  for (j in 1:ncol(cmbb)) {
    y1 <- x1[cmbb[,j]]
    y2 <- x2[cmbb[,j]]
    dattmp <- data.frame(y1, y2)
    names(dattmp) <- c(paste(cmba[1,i], cmbb[,j], sep=""), paste(cmba[2,i], cmbb[,j], sep=""))
    afctic <- factor(rep(cmba[,i], each=2))
    bfctic <- factor(rep(cmbb[,j], 2))
    idatic <- data.frame(afctic, bfctic)
    Anovaresic <- Anova(lm(as.matrix(dattmp)~1, dattmp), idata=idatic, idesign=~afctic*bfctic)
    Aresic <- Astat(Anovaresic)$unv
    lisnm <- paste(paste(paste(cmba[1,i], cmbb[,j], sep=""), collapse="_"), "_vs._",paste(paste(cmba[2,i], cmbb[,j], sep=""), collapse="_"), sep="")
    icreslist[[lisnm]] <- Aresic; n <- n+1
}
}
#names(icreslist) <- lisnm
icreslist
tmp <- lapply(icreslist, function(x) x[4,])
tmp2 <- t(data.frame(tmp))
pvholm <- p.adjust(tmp2[,6], "holm")
icres <- cbind(tmp2, pvholm)
colnames(icres) <- c("ss", "nmdf", "errorss", "dendf", "fv", "pv", "Pr(>F).holm")
class(icres) <- "anova"
icres


## プールされた誤差項で分析してみる。多重比較はやめておこう
Ares
Aunv
### それぞれの誤差項をとりだしておく
(mse.a <- (Aunv[2,3]+Aunv[4,3])/(Aunv[2,4]+Aunv[4,4]))
(mse.adf <- (Aunv[2,4]+Aunv[4,4]))
(mse.b <- (Aunv[3,3]+Aunv[4,3])/(Aunv[3,4]+Aunv[4,4]))
(mse.bdf <- (Aunv[3,4]+Aunv[4,4]))
(mse.ab <- (Aunv[4,3])/(Aunv[4,4]))
(mse.abdf <- (Aunv[4,4]))

## プールされた誤差項で単純主効果 ##
### 水準別誤差項の結果を利用する
# bにおけるAの単純効果
tmp <- list()
tmp <- lapply(seffect.a, function(x) Astat(x)$unv[2,])
tmp2 <- t(data.frame(tmp))
pvholm <- p.adjust(1-pf((tmp2[,1]/tmp2[,2])/mse.a, tmp2[,2], mse.adf), "holm")
seffect.a.pl <- cbind(tmp2[,1:2], mse.a*mse.adf, mse.adf, mse.a, (tmp2[,1]/tmp2[,2])/mse.a, 1-pf((tmp2[,1]/tmp2[,2])/mse.a, tmp2[,2], mse.adf), pvholm)
colnames(seffect.a.pl) <- c("ss", "nmdf", "errorss", "dendf", "mse", "fv", "pv", "Pr(>F).holm")
rownames(seffect.a.pl) <- paste("a effect at ", bnms, sep="")
class(seffect.a.pl) <- "anova"
seffect.a.pl

# Bの効果 at a水準
tmp <- list()
tmp <- lapply(seffect.b, function(x) Astat(x)$unv[2,])
tmp2 <- t(data.frame(tmp))
pvholm <- p.adjust(1-pf((tmp2[,1]/tmp2[,2])/mse.b, tmp2[,2], mse.bdf), "holm")
seffect.b.pl <- cbind(tmp2[,1:2], mse.b*mse.bdf, mse.bdf, mse.b, (tmp2[,1]/tmp2[,2])/mse.b, 1-pf((tmp2[,1]/tmp2[,2])/mse.b, tmp2[,2], mse.bdf), pvholm)
colnames(seffect.b.pl) <- c("ss", "nmdf", "errorss", "dendf", "mse", "fv", "pv", "Pr(>F).holm")
rownames(seffect.b.pl) <- paste("b effect at ", anms, sep="")
class(seffect.b.pl) <- "anova"
seffect.b.pl


## プールされた誤差項で交互作用対比…交互作用項の誤差項を使えばいいのか？
icres.pltp <- icres[,1:2]
pl.mss <- icres[,1]/icres[,2]
pl.errorss <- mse.ab*mse.abdf
pv.holm <- p.adjust(1-(pf(pl.mss/mse.ab, icres[,2], mse.abdf)), "holm")
icres.pl <- cbind(icres.pltp, pl.errorss, mse.abdf, pl.mss/mse.ab, 1-(pf(pl.mss/mse.ab, icres[,2], mse.abdf)), pv.holm)
colnames(icres.pl) <- c("ss", "nmdf", "errorss", "dendf", "fv", "pv", "Pr(>F).holm")
class(icres.pl) <- "anova"
icres.pl