fact.mc <- function(dat) {
dat <- dat
result <- dat[,3]; a <- dat[,1]; b <- dat[,2]
library(car); options(contrasts=c("contr.sum", "contr.sum"))
Anovares <- Anova(lm(result~a*b))
errorss <- Astat(Anovares)$unv[4,1]
errordf <- Astat(Anovares)$unv[4,2]
mse <- errorss/errordf
tvlist.a <- list()
tvlist.b <- list()
pvs <- vector()
datnames <- names(dat)
meandf <- aggregate(dat[3], list(dat[,1], dat[,2]), mean)
meandf.a <- aggregate(meandf[3], list(meandf[,1]), mean)
meandf.b <- aggregate(meandf[3], list(meandf[,2]), mean)
lendf <- aggregate(dat[3], list(dat[,1], dat[,2]), length)
lvname <- levels(dat[,1])
cbn <- combn(lvname, 2)

for (i in 1:ncol(cbn)) {
x1 <- lendf[,3][lendf[1]==cbn[1,i]]
x2 <- lendf[,3][lendf[1]==cbn[2,i]]
meanv <- meandf.a[,2]
factv <- meandf.a[,1]
diff1 <- meanv[factv==cbn[1,i]]
diff2 <- meanv[factv==cbn[2,i]]
diff <- diff1-diff2
den1 <- sum(c(mean(1/x1), mean(1/x2)))/length(x1)
den2 <- sqrt(mse*den1)
tv <- diff/den2
pv <- pt(abs(tv), errordf, lower.tail = F)*2
pvs[i] <- pv
ltval <- diff-(-1*qt(0.025, errordf)*den2)
utval <- diff+(-1*qt(0.025, errordf)*den2)
tvlist.a[[i]] <- c(diff1, diff2, diff, ltval, utval, mse, errordf, tv, pv)
}
pvs <- p.adjust(pvs, "holm")
tmptdf <- cbind(t(data.frame(tvlist.a)), pvs)
rownames(tmptdf) <- paste(cbn[1,], "_", cbn[2,], sep="")
colnames(tmptdf) <- c("mean1", "mean2", "diff", "lower", "upper", "mse", "errordf", "tv", "pv", "Pr(>t).holm")
res.a <- tmptdf; class(res.a) <- "anova"


lvname <- levels(dat[,2])
cbn <- combn(lvname, 2)
for (i in 1:ncol(cbn)) {
x1 <- lendf[,3][lendf[2]==cbn[1,i]]
x2 <- lendf[,3][lendf[2]==cbn[2,i]]
meanv <- meandf.b[,2]
factv <- meandf.b[,1]
diff1 <- meanv[factv==cbn[1,i]]
diff2 <- meanv[factv==cbn[2,i]]
diff <- diff1-diff2
den1 <- sum(c(mean(1/x1), mean(1/x2)))/length(x1)
den2 <- sqrt(mse*den1)
tv <- diff/den2
pv <- pt(abs(tv), errordf, lower.tail = F)*2
pvs[i] <- pv
ltval <- diff-(-1*qt(0.025, errordf)*den2)
utval <- diff+(-1*qt(0.025, errordf)*den2)
tvlist.b[[i]] <- c(diff1, diff2, diff, ltval, utval, mse, errordf, tv, pv)
}
pvs <- p.adjust(pvs, "holm")
tmptdf <- cbind(t(data.frame(tvlist.b)), pvs)
rownames(tmptdf) <- paste(cbn[1,], "_", cbn[2,], sep="")
colnames(tmptdf) <- c("mean1", "mean2", "diff", "lower", "upper", "mse", "errordf", "tv", "pv", "Pr(>t).holm")
res.b <- tmptdf; class(res.b) <- "anova"


reslist <- list()
reslist[[1]] <- res.a
reslist[[2]] <- res.b
names(reslist) <- names(dat) [1:2]
options(contrasts = c("contr.treatment", "contr.poly"))
return(reslist)
}


seffect <- function (dat) {
    dat <- dat
    result<- dat[,3]; a <- dat[,1]; b <- dat[,2]
    library(car)
    Anovares <- Anova(lm(result~a*b))
    errorss <- Astat(Anovares)$unv[4,1]
    errordf <- Astat(Anovares)$unv[4,2]
    mse <- errorss/errordf
    errorss <- Astat(Anovares)$unv[4,1]
    errordf <- Astat(Anovares)$unv[4,2]
    (mse <- errorss/errordf)
    lvsa <- levels(dat[,1]); nlvsa <- nlevels(dat[,1])
    lvsb <- levels(dat[,2]); nlvsb <- nlevels(dat[,2])
    selist.a <- list(); selist.b <- list()
    for (i in 1:nlvsa) {
      x <- subset(dat, dat[1]==lvsa[i])
      x2 <- Astat(Anova(lm(x[,3]~x[,2], x)))$unv
      x3 <- c(x2[1,1], x2[1,2], errorss, errordf, mse, (x2[1,1]/x2[1,2])/mse, 1-pf((x2[1,1]/x2[1,2])/mse, x2[1,2], errordf))
      selist.a[[i]] <- x3
    }
    for (i in 1:nlvsb) {
      x <- subset(dat, dat[2]==lvsb[i])
      x2 <- Astat(Anova(lm(x[,3]~x[,1], x)))$unv
      x3 <- c(x2[1,1], x2[1,2], errorss, errordf, mse, (x2[1,1]/x2[1,2])/mse, 1-pf((x2[1,1]/x2[1,2])/mse, x2[1,2], errordf))
      selist.b[[i]] <- x3
    }
    tmpdf.a <- t(data.frame(selist.a))
    class(tmpdf.a) <- "anova"
    bname <- names(dat)[2]
    rownames(tmpdf.a) <- paste(bname, "_at_", lvsa)
    colnames(tmpdf.a) <- c("ss", "nmdf", "errorss", "errordf", "mse", "Fvalue", "Pr(>F)")
    mnef.a <- tmpdf.a
    tmpdf.b <- t(data.frame(selist.b))
    class(tmpdf.b) <- "anova"
    aname <- names(dat)[1]
    rownames(tmpdf.b) <- paste(aname, "_at_", lvsb)
    colnames(tmpdf.b) <- c("ss", "nmdf", "errorss", "errordf", "mse", "Fvalue", "Pr(>F)")
    mnef.b <- tmpdf.b
    dflist <- list()
    dflist[[1]] <- mnef.b; dflist[[2]] <- mnef.a; dflist[[3]] <- Astat(Anovares)
    names(dflist) <- c(paste(c(aname, bname), "effect", sep="."), "OmnibusAnova")
    return(dflist)
}