R/nerve/MixedDiagnosticsRE.R

In premodelling/mixed: Modelling Multilevel & Longitudinal Problems

# Title     : MixedDiagnosticsRE.R
# Objective : Diagnostics w.r.t. random effects (RE)
# Created by: greyhypotheses
# Created on: 06/04/2022



GaussianAssumptionRE <- function (patients) {

  # ... does each term have a normally distributed set of values?
  patients %>%
    gather(key = 'Term', value = 'Sample') %>%
    ggplot(aes(sample = Sample)) +
    facet_wrap(~Term, scales = 'free', labeller = as_labeller(c('intercept' = 'Intercept', 'painScoreStd' = 'std(Pain Score)')) ) +
    stat_qq() +
    stat_qq_line() +
    theme_minimal() +
    theme(panel.spacing = unit(x = 3, units = 'line'),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_line(size = 0.05)) +
    xlab(label = '\nTheoretical Quantiles\n') +
    ylab(label = '\nSample Quantiles\n')

}


CorrelationAssumptionRE <- function (patients) {

  # ... does the correlation assumption hold?
  graph <- ggplot(data = patients, mapping = aes(x = painScoreStd, y = intercept)) +
    geom_point(alpha = 0.65) +
    geom_abline(intercept = 0, slope = 1) +
    theme_minimal() +
    theme(panel.grid.minor = element_blank(),
          panel.grid.major = element_line(size = 0.05)) +
    xlab(label = '\nRandom Coefficients\nw.r.t. std(pain score), patient\n') +
    ylab(label = '\nRandom Intercepts\nw.r.t. std(pain score), patient\n')

}

MeasuresRE <- function (model) {

  # "This function calculates the estimated variances, standard deviations, and
  #  correlations between the random-effects terms in a linear mixed-effects
  #  model.  The within-group error variance and standard deviation are also calculated."
  measures <- nlme::VarCorr(x = model)
  data.frame(measures$patient)


  # correlation
  correlations <- attr(x = measures$patient, which = 'correlation')
  correlation <- correlations[upper.tri(correlations)]


  # standard deviations
  deviations <- attr(x = measures$patient, which = 'stddev')


  # covariance
  # covariance <- correlation * prod(deviations)
  covariance <- measures$patient['(Intercept)', 'painScoreStd']


  return(list(variances = data.frame(measures$patient), correlation = correlation,
              deviations = deviations, covariance = covariance))

}