R/getPooledR2.R

In earlycapistran/consLettersUtils: Conservation Letters Utilities

# =============================================================================
# R^2 for multiply imputed repeated
# earlycapistran@comunidad.unam.mx - July 2020
# =============================================================================

#' Utility to calculate R^2  values by pooling values from 
#' multiply imputed repeated analyses. 
#' Based on van Buuren (https://rdrr.io/cran/mice/src/R/pool.r.squared.R)
#' 
#' @param mira An object of class 'mira' generated by 'mice'
#' @return R^2 value
#' @export
#' @usage
#' getPooledR2(mira)
#' 
#' @importFrom mice is.mira
#' @importFrom mice pool.scalar
#' 

getPooledR2 <- function(mira) {
  if (!is.mira(mira)) 
    stop("The object must have class 'mira'")
  if ((m <- length(mira$analyses)) < 2) 
    stop("At least two imputations are needed for pooling.\n")
  
  # Calculate quasi-R2 values for each of the multiply imputed repeated analyses
  analyses <- mira$analyses
  n <- length(analyses)
  r2<-list()
  for(i in seq_len(n))
  {
    r2[[i]] <- getMiceR2(analyses[[i]])
  } 
  
  # Set up array 'res2' to store R2 values, Fisher z-transformations of R2 values and its variance.
  res2 <- matrix(NA, nrow = n, ncol = 3, dimnames = list(seq_len(n), c("R^2", "Fisher trans F^2", "se()")))
  # Fill arrays
  for (i in seq_len(n)) {
    fit <- analyses[[i]]
    r.squared <- r2[[i]]
    res2[i, 1] <- sqrt(r.squared)
    res2[i, 2] <- 0.5 * log((res2[i, 1] + 1)/(1 - res2[i, 1]))
    res2[i, 3] <- 1/(length(residuals(fit) - 3))
  }
  
  # Compute within, between and total variances following Rubin's rules  with function pool.scalar().
  fit <- pool.scalar(res2[, 2], res2[, 3])
  
  # Make table with results.
  qbar <- fit$qbar
  table <- array(((exp(2 * qbar) - 1)/(1 + exp(2 * qbar)))^2,
                 dim = c(1, 3))
  
  dimnames(table) <- list("R^2", c("est", "lo 95", "hi 95"))
  
  table[, 2] <- ((exp(2 * (qbar - 1.96 * sqrt(fit$t))) - 1)/(1 + exp(2 * (qbar - 1.96 * sqrt(fit$t)))))^2
  table[, 3] <- ((exp(2 * (qbar + 1.96 * sqrt(fit$t))) - 1)/(1 + exp(2 * (qbar + 1.96 * sqrt(fit$t)))))^2
  return(table)
}