mcmcpvalue <- function(X) {
   ##From Bates pers comm. September 14, 2006 2:59:23 PM EDT
   ## elementary version that creates an empirical p-value for the
   ## hypothesis that the columns of samp have mean zero versus a
   ## general multivariate distribution with elliptical contours.
# samp <- rnorm(10000, m=3)
   ## differences from the mean standardized by the observed
   ## variance-covariance factor
if(class(X) == "numeric") samp <- as.matrix(X) else{
  samp <- t( X ) }
std <- backsolve(chol(var(samp)),
                    cbind(0, t(samp)) - colMeans(samp),
                    transpose = TRUE)
   sqdist <- colSums(std * std)
   sum(sqdist[-1] > sqdist[1])/nrow(samp)
}

aov.Bayes <- function(mod, mod.mc=NULL, n=100, mcSamp=TRUE, ...) {
  ## after ideas by Bates and others Sept. 2005
## mod.mc is mcmcsamp output
if(any(attr(mod @ X, "contrasts") == "contr.treatment") == TRUE) {
  print("Danger! Danger! Contrasts not orthogonal! This will alter the conclusions from the Bates.Bayes empirical P-values") }

  #mc.out <- mixed.mod.out.fixed(mod, N=n)
  sources <- attr(mod @ X, "assign")
  aov.lmer <- anova(mod)
  mod.reml <- summary(mod)
  df1 <- dim(mod @ X)[1] - dim(mod @ X)[2]

  dfhatTr <- dim(mod @ X)[1] - sum(aov.lmer[,1])
  MS.Resid <- mod.reml @ sigma ^ 2
  F <- aov.lmer$"Mean Sq"/MS.Resid
  P.sum <-  1-pf(F, aov.lmer$Df, dfhatTr)
  aov.out <- data.frame(aov.lmer, P.sum=round(P.sum,4)) 
  nextrows <- dim(aov.out)[1]+1
  aov.out[nextrows,] <- matrix( c(dfhatTr, NA, MS.Resid, NA, NA), nrow=1)
  rownames(aov.out)[nextrows] <- "Residuals (Tr(hat))"
  ## mc.out is the output of mixed.mod.out
  ## sources is a list which includes the groups of coefficients
  ## that comprise a "factor"
if(mcSamp==TRUE){
  mod.mc <- mcmcsamp(mod, n=n)
hpds <- lme4::HPDinterval(mod.mc)}
  ## Bayes <- sapply(sources, function(x) if(length(x)==1) {
  ## mcmcpvalue1(mcsam[,x])} else mcmcpvalue(mcsam[,x]) )
  Bayes <- sapply(unique(sources[-1]), function(x)  {
    rows <- which(sources==x)
    mcmcpvalue(mod.mc@fixef[rows,])}
                  )
  #aov.out$"P(H|D)" <- round(c(Bayes, NA), 4)
  aov.out$"P(H|D)" <- round(c(Bayes, NA), floor(log10(n)))

  list(aov.tab=aov.out, coeff = hpds)
}