anovaMI <- function(obj1, obj2) {
	library(mice)
	stopifnot(is(obj1, "mira"))
	stopifnot(is(obj2, "mira"))
	out1 <- obj1$analyses
	out2 <- obj2$analyses
	df1 <- attributes(logLik(out1[[1]]))$df
	df2 <- attributes(logLik(out2[[1]]))$df
	chis <- mapply(function(x, y) abs(logLik(x) - logLik(y)), x=out1, y=out2)
	result <- lmrrPooledChi(chis, abs(df1 - df2))
	print(round(result, digits = 3))
	invisible(result)
}

lmrrPooledChi <- function(chis, df) {
	# From Li, Meng, Raghunathan, & Rubin (1991)
	if(is.matrix(chis)) {
		ifelse(ncol(chis) == 1 | nrow(chis) == 1, chis <- as.vector(chis), stop("Please put a vector of chi-square values"))
	}
	m <- length(chis)
	dbar <- mean(chis)
	sqrtd <- sqrt(chis)
	xbarsqrtd <- mean(sqrtd)
	# Equation 2.2
	r <- (1 + 1/m) * (sum((sqrtd - xbarsqrtd)^2)/(m - 1))
	# Equation 2.1
	D <- (dbar/df - ((m + 1) * r /(m - 1)))/(1 + r)
	if(D < 0) D <- 0
	# Equation 2.16 and 2.17
	aw <- df^(-(3/m)) * (m - 1) * (1 + (1/r))^2
	p <- 1 - pf(D, df, aw)
	result <- c(D, df, aw, p)
	names(result) <- c("F", "df1", "df2", "p.F")
	return(result)
}
#Examples:
#lmrrPooledChi(c(89.864, 81.116,71.500,49.022,61.986,64.422,55.256,57.890,79.416,63.944), 2)

ranefMI <- function(obj) {
  library(mice)
  stopifnot(is(obj, "mira"))
  out <- obj$analyses
  sumout <- lapply(out, summary)
  vc <- lapply(sumout, "[[", "varcor")
  vc <- sapply(vc, rbind)
  result <- Reduce("+", vc) / length(vc)
  result2 <- result
  attr(result2, "stddev") <- NULL 
  attr(result2, "correlation") <- NULL
  attr(result, "stddev") <- sqrt(diag(result2))
  attr(result, "correlation") <- cov2cor(result2)
  return(result)
}

sigmaMI <- function(obj) {
  library(mice)
  stopifnot(is(obj, "mira"))
  out <- obj$analyses
  sumout <- lapply(out, summary)
  return(mean(sapply(sumout, "[[", "sigma")^2))
}