R/getPooledSE.R

In earlycapistran/consLettersUtils: Conservation Letters Utilities

# =============================================================================
# Pooled standard error from MICE
# earlycapistran@comunidad.unam.mx - March 2021
# =============================================================================

#' @title Pooled standard error from 'mice'
#' 
#' @description Calculates and pools standard errors from multiply imputed 
#' repeated analyses, generated by package 'mice', according to Rubin's Rules. T
#' his function assumes normally distributed errors. 
#'
#' @param mira A 'mira' object generated by package 'mice'
#' @return A numerical value (standard error)
#' @export 
#' @usage
#' getPooledSE(mira)
#' 
#' @importFrom mice is.mira
#' @importFrom broom augment

getPooledSE <- function(mira) {
  if (!is.mira(mira)) 
    stop("The object must have class 'mira'")
  # Calculate predictions from m imputed models
  predList <- lapply(mira$analyses, broom::augment)
  # Select vector of fitted values from each list item
  predList <- lapply(predList, `[[`, ".fitted")
  m <- length(predList)
  means <- sapply(predList, mean)
  meanBar <- mean(means)

  # Calculate within-models variance
  sse <- sapply(predList, sd)/sqrt(m) #Squared SE for each dataset
  sse <- sse^2
  # Pool across datasets
  wV <- mean(sse)
  
  # Calculate between-imputation variance (mean differences)
  bV <- (sum((means - meanBar)**2))/(m-1) 
  
  # Calculate total variance
  vT <- wV + bV + bV/(length(predList))
  
  # Calculate pooled standard error 
  seT <- sqrt(vT)
  
  result <- cbind("Pooled Standard Error" = seT)
  rownames(result) <- c("")
  print(result)
  return(seT)
}