R/obsImpResiPlot.R

In earlycapistran/consLettersUtils: Conservation Letters Utilities

# =============================================================================
# Plot of residuals vs fitted values for imputed and observed models
# earlycapistran@comunidad.unam.mx - March 2021
# =============================================================================

#' @title Plot residuals vs fitted for imputed and observed models
#' 
#' @description Provides a plot with residuals and fitted values from all 
#' multiply imputed models (grey points) and from the observed data model 
#' (black points).
#' 
#' Based on method proposed by: 
#' Rizopoulos, D., Verbeke, G., & Molenberghs, G. (2010). Multiple-Imputation-
#' Based Residuals and Diagnostic Plots for Joint Models of Longitudinal and 
#' Survival Outcomes. Biometrics, 66(1), 20-29. 
#' https://doi.org/10.1111/j.1541-0420.2009.01273.x
#'  
#' @param mira An object of class 'mira' generated by the 'with.mids' function 
#' in ' mice'
#' @param nls A fitted model object created with observed data only (e.g., an
#' nls model object created with observed dataset that has missing values).
#' @return A plot of residual vs. fited values
#' @export
#' @usage
#' obsImpResiPlot(mira, nls)
#' 
#' @importFrom mice is.mira
#' @import graphics

obsImpResiPlot <- function(mira, nls) {
  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")

  # Extract residual and predicted values from each of the 'm' 
  # complete-data models ------------------------------------------------------
  miResi <- sapply(mira$analyses, residuals)  
  miPred <- sapply(mira$analyses, predict)  

  # Extract residual and fitted values from the observed data medel -----------
  obsPred <- fitted(nls)
  obsResi <- residuals(nls)

  # Calculate ranges for x and y axes
  xlim <- range(c(miPred,obsPred))
  ylim <- range(c(miResi,obsResi))

  # Plot fitted vs. residual values for all imputed data models
  plot(x = miPred, y = miResi,
       col = "grey",
       abline(h=0),
       xlim=xlim, ylim=ylim,
       xlab = "Fitted values", ylab = "Residuals",
       main = "Residuals versus fitted values")

  par(new=TRUE)
  # Add points for fitted vs. residual values for observed model ---------------
  points(x = obsPred, y = obsResi, pch = 19)
}