R/caterpillar_mean_plot.R

In dacarras/r4sda: Collection of functions to produce descriptives analysis of secondary data and large scale assessments studies

#' caterpillar_mean_plot() it plots the rank of random effects across clusters, with ther point estimates and their uncertainty (1.96*SD), around it points estimates
#'
#' @param mlm_model it requires an mlm_model, generated by lme4. If this is a null, the plot display the unconditioned random effect, plus its intercept. If the model is not a null model, then the plot display the expected random effects given the model plus its intercept or grand mean.
#'
#' @return a plot, a `ggplot2` object.
#'
#' @details Given the return object is a common `ggplot2` user can further customize the figure, overwritting or defining parameter for the plot such as colours, axis lenght, and coord_fixed(ratio=.xx)
#'
#' @examples
#'
#' lme4::lmer(y ~ 1 + (1|id_j), data = clustered_data, REML = FALSE) %>%
#' caterpillar_mean_plot(.)
#'
#' @export
caterpillar_mean_plot <- function(mlm_model){
  ## required libraries
  require(ggplot2)
  require(lme4)
  require(dplyr)
  require(tibble)
  require(broom.mixed)

  ## extract intercept
  intercept <- broom.mixed::tidy(mlm_model) %>%
    dplyr::filter(effect == 'fixed') %>%
    dplyr::select(estimate) %>%
    .$estimate %>%
    as.numeric()

  ## extract random effects
  rand_intercept <- tibble::as_tibble(lme4::ranef(mlm_model)) %>%
    dplyr::filter(term == '(Intercept)') %>%
    mutate(id_j = as.numeric(as.character(grp))) %>%
    mutate(id_j = seq(1:nrow(.))) %>%
    mutate(id_j = as.numeric(id_j)) %>%
    mutate(u_j = intercept + condval) %>%
    mutate(ll = u_j-condsd*1.96) %>%
    mutate(ul = u_j+condsd*1.96) %>%
    dplyr::select(id_j, u_j,ll, ul)

  ## get icc
  var_com <- broom.mixed::tidy(mlm_model) %>%
    dplyr::filter(effect == 'ran_pars') %>%
    dplyr::mutate(var = estimate^2) %>%
    dplyr::mutate(per = var/sum(var))

  icc <- as.numeric(dplyr::select(var_com, per)[1,])

  ## plot function
  p <- ggplot(rand_intercept, aes(x=reorder(id_j,u_j), y=u_j, ymin=ll, ymax=ul))+
    geom_pointrange(colour='grey20', alpha = .25, size = .2)+
    geom_hline(yintercept = intercept, linetype=2, size = .25, colour = "grey50") +
    theme(plot.title = element_text(hjust = 0.5)) +
    theme(axis.text.y = element_text(size=6, colour = "grey50")) +
    theme(axis.text.x = element_blank()) +
    theme(axis.ticks = element_blank())+
    theme(panel.background = element_rect(fill = "white", colour = "grey50"))+
    theme(panel.grid.major.y = element_line(size = .05, colour = "grey50")) +
    theme(panel.grid.minor.y = element_line(size = .05, colour = "grey50")) +
    xlab('cluster ranking') +
    ylab('latent means') +
    annotate('text',
             # positions the ICC at the first fifth of the rank the clusters
             x = nrow(rand_intercept)/5,
             # positions the ICC at the centered grand mean
             y = intercept,
             label = paste0('ICC = ', format(round(icc, 2), nsmall = 2)))
  ## return plot
  return(p)
}