R/clmm_analysis.R

In seonjoo/slmisc: Miscelaneous functions to conduct LeeLab projects

#' function: ordinal logistic regression, Cumulative Link "Mixed" Models
#' first version of clmm_analysis
#'
#' @param fits List. Fitted ordinal logistic regression from clmm package.
#' @param xvar string. explanatory variable to be analyzed.
#' @return List of emmeans and lsdiff (group differences) and figure of lsmeans vs time points
#' @examples TBA
#' @author Jongwoo Choi, \email{jc4816@columbia.edu}
#' @references TBA
#' @keywords ordinal logistic regression cumulative link mixed models anova lsmean emmean random effect
#' @import emmeans
#' @import ggplot2
#' @import knitr
#' @export

clmm_analysis <- function(fits, xvar){

  # Here we will create an object of emmeans ouput called marginal.
  # The formula in the emmeans function indicates that pairwise comparisons should be conducted for the xvar
  marginal = eval(parse(text=paste('emmeans(fits[[xvar]], pairwise~', xvar, "*eval, mode='linear.predictor')",  sep='')))  # mode can be changed
  # See details here:https://cran.r-project.org/web/packages/emmeans/vignettes/models.html#O

  emdiff <- confint(marginal$contrasts) # for group differences
  emdiff.intersts <- as.data.frame(emdiff[c(1,14,23,28),]) # interest at each time point. # Need to be generalized

  estimate.interests <- exp(emdiff.intersts$estimate)
  # var2 <-as.character(emdiff.intersts[,'contrast'])
  time_point <- c('0','1','2','3')

  #ci <- emdiff.intersts$estimate + qnorm(0.95)* sqrt(emdiff.intersts$SE) %o% c(-1, 1)
  # Can we go with asymptotic confidence interval since we have large sample size?

  # plot
  fig <- ggplot(emdiff.intersts, aes(x=time_point, y=estimate.interests)) +
    geom_errorbar(aes(ymin=exp(asymp.LCL), ymax=exp(asymp.UCL)), size=1, width=.2) +
    geom_point(size=3) +
    geom_hline(yintercept=exp(0), col='red') +
    labs(x='Time Point', y='Odds Ratio of EMMean differences')  # can be changed
  #+ ggtitle(paste('EMMeans,',xvar, 'vs eval')) # can be changed

  return_list <- list(emdiff.intersts, fig)
  names(return_list) <- c('emdiff', 'plot')
  return(return_list)
}