Alternative version of mitools:::MIcombine.default

AlternativeMIcombine.R

## Implement multivariate version of Rubin's formula
## eq 14 in https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4029775/
## Page 97 of https://www.amazon.com/Multiple-Imputation-Nonresponse-Surveys-Donald/dp/0471655740
AlternativeMIcombine.default <- function (results, variances, call = sys.call(), df.complete = Inf, ...) {
    m <- length(results)
    oldcall <- attr(results, "call")
    if (missing(variances)) {
        variances <- suppressWarnings(lapply(results, vcov))
        results <- lapply(results, coef)
    }
    vbar <- variances[[1]]
    cbar <- results[[1]]
    for (i in 2:m) {
        cbar <- cbar + results[[i]]
        vbar <- vbar + variances[[i]]
    }
    cbar <- cbar/m
    vbar <- vbar/m
    evar <- var(do.call("rbind", results))

    ## This part may need change?
    ## evar/vbar is element wise division, not a matrix operation.
    r <- (1 + 1/m) * evar/vbar

    df <- (m - 1) * (1 + 1/r)^2
    if (is.matrix(df))
        df <- diag(df)
    if (is.finite(df.complete)) {
        dfobs <- ((df.complete + 1)/(df.complete + 3)) * df.complete *
            vbar/(vbar + evar)
        if (is.matrix(dfobs))
            dfobs <- diag(dfobs)
        df <- 1/(1/dfobs + 1/df)
    }
    if (is.matrix(r))
        r <- diag(r)

    ##
    ## Original is different from Rubin 1987
    ## varest <- vbar + evar * (m + 1)/m
    ## Version consistent with Rubin 1987
    ## Ubar + (1 + 1/m) * tr(B Ubar^-1) / p * Ubar
    varest <- vbar + (m + 1)/m * sum(diag(evar %*% solve(vbar))) / length(cbar) * vbar

    ## Construct results
    rval <- list(coefficients = cbar, variance = varest,
                 call = c(oldcall, call), nimp = m, df = df,
                 missinfo = (r + 2/(df + 3))/(r + 1))
    class(rval) <- "MIresult"
    rval
}