myscratchpad

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/mixed.htm 感謝。
dat <- data.frame(
cnd = factor(rep(c("wt", "atp", "trt"), each=4)),
pre=c(13, 10, 13, 4, 5, 8, 14, 12, 13, 9, 14, 8),
hlf=c(14, 11, 19, 12, 10, 15, 16, 21, 24, 22, 22, 18),
pst=c(20, 14, 21, 15, 21, 24, 23, 26, 30, 24, 28, 28),
fu=c(17, 15, 18, 14, 17, 22, 23, 26, 28, 22, 28, 27)
)

by(dat[,2:5], list(dat[,1]), mean)
by(dat[,2:5], list(dat[,1]), sd)

## とりあえずプロット
x <-aggregate(dat[2:5], list(dat[,1]), mean)
gx <- as.matrix(x[2:5])
rownames(gx) <- names(dat[,1])
barplot(gx, beside=T)
dev.off()

## 水準名とか
ncl <- ncol(dat)
anms <- levels(dat[,1])
  nanms <- sapply(anms, function(x) sum(x==dat[,1]))
  cmba <- combn(anms, 2)
bnms <- names(dat[2:ncl])
  mmbs <- mean(aggregate(dat[,2:ncl], list(dat[,1]), mean)[,2:ncl])
  cmbb <- combn(bnms, 2)

## 全体の分散分析
options(contrasts = c("contr.sum", "contr.sum"))
lmres <- lm(as.matrix(dat[,2:ncl]) ~ dat[,1] , data = dat)
bfact <- factor(bnms)
idat <- data.frame(bfact)
library(car)
Anovares <- Anova(lmres, idata = idat, idesign = ~ bfact)
(Ares <- Astat(Anovares)$unv)


## 要因の主効果の多重比較: 参加者間要因
val <- rowMeans(dat[,2:ncl])
a <- dat[,1]
factmc.a <- holm.mc(data.frame(a, val))
factmc.a

## 要因の主効果の多重比較: 参加者内要因
factmc.blist <- list()
ps <- vector()
library(car)
for (i in 1:ncol(cmbb)) {
  mcbfact <- factor(cmbb[,i])
  mcidata <- data.frame(mcbfact)
  mclmres <- lm(as.matrix(dat[,c(cmbb[,i])])~dat[,1], dat)
  mcAnovares <- Anova(mclmres, idata=mcidata, idesign=~mcbfact)
  mcAres <- Astat(mcAnovares)$unv
  mcmse <- mcAres[2,3]/mcAres[2,4]
  mcdf <- mcAres[2,4]
  den1 <- (mean(1/nanms)/nlevels(dat[,1]))*2
  diff <- mmbs[cmbb[1,i]] - mmbs[cmbb[2,i]]
  tv <- diff/sqrt(mcmse*den1)
  ps[i] <- pt(abs(tv), mcdf, lower.tail = F)*2
  resv <- c(mmbs[cmbb[1,i]], mmbs[cmbb[2,i]], diff, mcmse, mcdf, tv, pt(abs(tv),   mcdf, lower.tail = F)*2)
  factmc.blist [[i]] <- resv
}
factmc.b <- t(data.frame(factmc.blist ))
ps.holm <- p.adjust(ps, "holm")
factmc.b <- cbind(factmc.b, ps.holm)
colnames(factmc.b) <- c("m1", "m2", "diff", "mse", "df", "tv", "pv", "Pr(>|t|).holm")
rownames(factmc.b) <- paste(cmbb[1,], "_", cmbb[2,], sep="")
class(factmc.b) <- "anova"
factmc.b


## 単純主効果 ##
## 水準別誤差項で単純主効果 ##
## b (参加者内水準) におけるaの効果とその多重比較
seffect.a <- list()
seffect.a <- lapply(dat[,2:ncl], function(x) Anova(lm(x~a, dat)))
mcse.a <- list()
for (nm in bnms) {
  val <- dat[nm]
  x <- dat[,1]
  mcse.a[[nm]] <- holm.mc(data.frame(x, val))
}
names(seffect.a) <- paste(names(dat)[1], " effect at ", bnms, sep="")
names(mcse.a) <- paste(names(dat)[1], " effect at ", bnms, sep="")
seffect.a # 単純主効果
mcse.a # その多重比較


## a (参加者間水準) におけるBの効果とその多重比較
seffect.b <- list(); seffect.b <- list()
mcse.b <- list()
for (nm in anms) {
  dattmp <- subset(dat, dat[,1]==nm)
  Anovares.sepa <- Anova(lm(as.matrix(dattmp[,2:ncl])~1, data=dattmp), idata=idat, idesign=~bfact)
  seffect.b[[nm]] <- Anovares.sepa
  mcse.b[[nm]] <- holm.mc(dattmp[,2:ncl], datw=T, psd=F, paired=T)
}
seffect.b <- lapply(seffect.b, function(x) Astat(x)$unv)
names(seffect.b) <- paste(names(dat)[2], " effect at ", anms, sep="")
names(mcse.b) <- paste(names(dat)[2], " effect at ", anms, sep="")
seffect.b # 単純主効果
mcse.b # 多重比較


## 交互作用対比
### 処理用オブジェクト
iclist <- list()
ps <- vector()
rnms <- vector()
cmba
cmbb
icn <- 1
### 交互作用対比の分析
for (i in 1:ncol(cmba)) {
  dattmpa <- subset(dat, dat[,1]==cmba[1,i]|dat[,1]==cmba[2,i])
  for (j in 1:ncol(cmbb)) {
    dattmp <- data.frame(dattmpa[1], dattmpa[cmbb[,j]])
    bfacttmp <- factor(cmbb[,j])
    idattmp <- data.frame(bfacttmp)
    lmrestmp <- lm(as.matrix(dattmp[,2:3])~dattmp[,1], dattmp)
    Anovarestmp <- Anova(lmrestmp, idata=idattmp, idesign=~bfacttmp)
    Arestmp <- Astat(Anovarestmp)$unv
    iclist[[icn]] <- Arestmp[3,]; ps[icn] <- Arestmp[3,6]; rnms[icn] <- paste(paste(cmba[,i], collapse=""), " vs. ", paste(cmbb[,j], collapse=""), sep="")
    icn <- icn+1
  }
}
icres <- t(data.frame(iclist))
ps.holm <- p.adjust(ps, "holm")
icres <- cbind(icres, ps)
colnames(icres) <- c("ss", "nmdf", "errorss", "dendf", "fv", "pv", "Pr(>F).holm")
rownames(icres) <- rnms
class(icres) <- "anova"
icres # それぞれのmseはicres[,3]/icres[,4]



## プールされた誤差項で分析してみる
Ares
### それぞれの誤差項をとりだしておく
mse.a <- Ares[1,3]/Ares[1,4] # 参加者間要因のmse: aの主効果の検定に使う
df.a <- Ares[1,4]
mse.ap <- sum(Ares[c(1,3), 3])/sum(Ares[c(1,3), 4]) # 参加者間要因を含むプールしたmse: aの単純主効果検定につかう
df.ap <- sum(Ares[c(1,3), 4]) # mse.aの自由度
mse.b <- Ares[2,3]/Ares[2,4] # 参加者内要因のプールしたmse
df.b <- Ares[2,4] # 参加者内要因の自由度
mse.ab <- Ares[3,3]/Ares[3,4] # 交互作用項のプールしたmse: bの単純主効果検定に使う
df.ab <- Ares[3,4] # mse.abの自由度

## 要因の主効果の多重比較: 参加者間要因
val <- rowMeans(dat[,2:ncl])
a <- dat[,1]
factmc.a.pl <- holm.mc(data.frame(a, val), mse=mse.a, mse.df=df.a)
factmc.a.pl

## 要因の主効果の多重比較: 参加者内要因
factmc.b.pllist <- list()
ps <- vector()
library(car)
for (i in 1:ncol(cmbb)) {
  mcmse <- mse.ab
  mcdf <- df.ab
  den1 <- (mean(1/nanms)/nlevels(dat[,1]))*2
  diff <- mmbs[cmbb[1,i]] - mmbs[cmbb[2,i]]
  tv <- diff/sqrt(mcmse*den1)
  ps[i] <- pt(abs(tv), mcdf, lower.tail = F)*2
  resv <- c(mmbs[cmbb[1,i]], mmbs[cmbb[2,i]], diff, mcmse, mcdf, tv, pt(abs(tv),   mcdf, lower.tail = F)*2)
  factmc.b.pllist [[i]] <- resv
}
factmc.b.pl <- t(data.frame(factmc.b.pllist ))
ps.holm <- p.adjust(ps, "holm")
factmc.b.pl <- cbind(factmc.b.pl, ps.holm)
colnames(factmc.b.pl) <- c("m1", "m2", "diff", "mse", "df", "tv", "pv", "Pr(>|t|).holm")
rownames(factmc.b.pl) <- paste(cmbb[1,], "_", cmbb[2,], sep="")
class(factmc.b.pl) <- "anova"
factmc.b.pl



## プールされた誤差項で単純主効果 ##
## b (参加者内水準) におけるAの単純効果
seffect.a.pllist <- list()
for (nm in bnms) {
  Anovares.se <- Anova(lm(dat[,nm]~a, dat))
  Ares.se <- Astat(Anovares.se)$unv
  seffect.a.pllist[[nm]] <- c(unlist(Ares.se[1,1:2]), mse.ap, df.ap, ((Ares.se[1,1]/Ares.se[1,2])/mse.ap), 1-pf(((Ares.se[1,1]/Ares.se[1,2])/mse.ap), Ares.se[1,2], df.ap))
 }
seffect.a.pl <-t(data.frame(seffect.a.pllist))
colnames(seffect.a.pl) <- c("ss", "nmdf", "mse", "dendf", "fv", "Pr(>F)")
rownames(seffect.a.pl) <- paste(names(dat)[1], " effect at ", bnms, sep="")
class(seffect.a.pl) <- "anova"
seffect.a.pl
## 多重比較は勘弁してやろう

## a (参加者間水準) におけるBの単純効果
seffect.b.pllist <- list()
for (nm in anms) {
  dattmp <- subset(dat, dat[,1]==nm)
  Anovares.se <- Anova(lm(as.matrix(dattmp [,2:ncl])~1, data=dattmp)
, idata=idat, idesign=~bfact, type=2)
  Ares.se1 <- Astat(Anovares.se)$unv
  if (rownames(Ares.se1)[1]=="(Intercept)") {Ares.se <- Ares.se1[2:nrow(Ares.se1),]} else {Ares.se <- Ares.se1}

  seffect.b.pllist[[nm]] <- c(unlist(Ares.se[1:2]), mse.b, df.b, ((Ares.se[1]/Ares.se[2])/mse.b), 1-pf(((Ares.se[1]/Ares.se[2])/mse.b), Ares.se[2], df.b))
  }
seffect.b.pl <-t(data.frame(seffect.b.pllist))
colnames(seffect.b.pl ) <- c("ss", "nmdf", "mse", "dendf", "fv", "Pr(>F)")
rownames(seffect.b.pl) <- paste(names(dat)[2], " 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*df.ab
pv.holm <- p.adjust(1-(pf(pl.mss/mse.ab, icres[,2], df.ab)), "holm")
icres.pl <- cbind(icres.pltp, pl.errorss, df.ab, pl.mss/mse.ab, 1-(pf(pl.mss/mse.ab, icres[,2], df.ab)), pv.holm)
colnames(icres.pl) <- c("ss", "nmdf", "errorss", "dendf", "fv", "pv", "Pr(>F).holm")
class(icres.pl) <- "anova"
icres.pl