consLettersUtils/R/analyseMiraResi.R

In earlycapistran/miceNls: Conservation Letters Utilities

# =============================================================================
# Residual analysis for multiply imputed repeated analyses
# earlycapistran@comunidad.unam.mx - July 2020
# =============================================================================

#' Provides separate residual analyses for each complete-data model generated
#' by with 'with.mids' function in 'mice'.
#' 
#' @param mira An object of class 'mira' generated by the 'with.mids' function
#' in 'mice'
#' @return Residual plots (normality, residuals vs. fitted values, lag plot)
#' Shapiro-Wilk normality test, Levene Test for homogeneity of variance, 
#' Runs test for randomness for each of the 'm' repeated complete-data models
#' 
#' @export
#' @usage
#' analyseMiraResi(mira)
#' 
#' @importFrom mice is.mira
#' @importFrom stats shapiro.test
#' @importFrom stats residuals
#' @importFrom stats predict
#' @importFrom stats qqnorm
#' @importFrom stats qqline
#' @importFrom car leveneTest
#' @import graphics

# To run this function, you must have 'car' installed
analyseMiraResi  <- 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")
  analyses <- mira$analyses
  n <- length(analyses)
  for(i in seq_len(n))
  {
    # Prepare residual and predicted values from each complete data model -----
    resi <- residuals(analyses[[i]]) 
    resDf <- as.data.frame(resi)
    pred <- predict(analyses[[i]]) 
    mean <- mean(resDf$resi)
    # Subset positive and negative values 
    resDf$sign <- as.factor(ifelse(resDf < 0, "negative", "positive"))

    # Make plots --------------------------------------------------------------
    # Plot normality
    par(mfrow = c(2, 2))
    qqnorm(resDf$resi, 
           main = paste("Normal Q-Q Plot", "(Model", i, ")"))
    qqline(resDf$resi)
    # Plot residual vs. fitted values
    res.plot <- plot(x = pred, y = resi,
                     xlab = "Fitted values", ylab = "Residuals",
                     main = paste("Fitted vs. residual", "(Model", i, ")"),
                     abline(h=0))
    
    # Run tests ---------------------------------------------------------------
    norm <- stats::shapiro.test(resDf$resi)
    levene <- car::leveneTest(resDf$resi ~ sign,
                              data = resDf)
    result <- list(normality = norm, mean = mean, levene = levene)
    names(result) <- c("Residual Normality Test", "Residual Mean", 
                       "Levene's Test")
    cat("Model", i,"\n")
    print(result)
  }
}