holm.mc <- function (dat, mse=0, mse.df= 0, psd=TRUE, datw=FALSE, paired=FALSE, alternative=c("two.sided", "less", "grater"), var.equal = FALSE, mtd = c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY","fdr", "none")) {
  if (datw==TRUE) {dat <- stack(dat)[2:1]}
  x <- dat[,2]
  g <- dat[,1]
  
        xbar <- tapply(x, g, mean, na.rm = TRUE)
        s <- tapply(x, g, sd, na.rm = TRUE)
        n <- tapply(!is.na(x), g, sum)
        degf <- n - 1
        alternative <- match.arg(alternative)
        mtd <- match.arg(mtd)
        total.degf <- sum(degf)
        cbn <- combn(levels(g), 2)
        tvallist <- list()
        dfmse <- mse
        if (psd) {
        if (mse==0) {pooled.sd <- sqrt(sum(s^2 * degf)/total.degf)} else {pooled.sd <- sqrt(mse); total.degf<- mse.df}
        plsd <- paste("pooled=", round(pooled.sd,4), sep=""); mse <- pooled.sd^2
        for (i in 1:ncol(cbn)) {
              m1 <- xbar[cbn[1,i]]; m2 <- xbar[cbn[2,i]]
              dif <- m1-m2
              se.dif <- pooled.sd * sqrt(1/n[cbn[1,i]] + 1/n[cbn[2,i]])
              t.val <- dif/se.dif
                if (alternative == "two.sided") 
                    pv <- 2 * pt(-abs(t.val), total.degf)
                else pv <- pt(t.val, total.degf, lower.tail = (alternative == 
                    "less"))
               ltval <- dif-(-1*qt(0.025, total.degf)*se.dif)
               utval <- dif+(-1*qt(0.025, total.degf)*se.dif)
               
           tvallist[[i]] <- c(m1, m2, dif, ltval, utval, pooled.sd^2, total.degf, t.val, pv)
         }
         
         }
         
         if (!psd) {
        plsd <- "non-pooled"
        for (i in 1:ncol(cbn)) {
          v1 <- na.omit(x[g==cbn[1,i]]); v2 <- na.omit(x[g==cbn[2,i]])
          m1 <- mean(v1); m2 <- mean(v2)
          dif <- m1-m2
          tres <- t.test(v1, v2, paired=paired, alternative=alternative, var.equal=var.equal)
          t.val <- tres$statistic
          pv <- tres$p.value
          ltval <- tres$conf.int[1]
          utval <- tres$conf.int[2]
          total.degf <- tres$parameter
          tvallist[[i]] <- c(m1, m2, dif, ltval, utval, mse, total.degf, t.val, pv)
             }}
         names(tvallist) <- apply(cbn, 2, function(x) paste(x, collapse="_"))
         tvallist <- t(data.frame(tvallist))
         pvh <- p.adjust(as.numeric(tvallist[,9]), method=mtd)
         tvallist <- cbind(tvallist, pvh)
         if (mse==0) {msea <- "mse=no"} else {msea <- paste("mse=", round(mse, 4), ", and dfmse=", round(dfmse, 4), sep="")}
         attr(tvallist, "heading") <- paste(plsd, msea, paste("paired=", paired, sep=""), alternative, paste("var.equal=", var.equal, sep=""), mtd, sep=",")
         class(tvallist) <- "anova"
         colnames(tvallist) <- c("mean1", "mean2", "diff", "lower", "upper", "mse", "errordf", "tv", "pv", paste("Pr(>t).", mtd, sep=""))
         return(tvallist)
         }