R/lmer_analysis_two_way_mixed_effects_model.R

In seonjoo/slmisc: Miscelaneous functions to conduct LeeLab projects

#' function: Two-way linear mixed effects model analysis:
#' first version of lmer analysis
#'
#' @param fits List. Fitted TWO-WAY linear mixed effect model.
#' @param xvar string. explanatory variable to be analyzed.
#' @return List of anova, lsmeans and lsdiff (group differences) and figure of lsmeans vs time points
#' @examples TBA
#' @author Jongwoo Choi, \email{jc4816@columbia.edu}
#' @references TBA
#' @keywords twoway anova mixedeffect lsmean
#' @import dplyr
#' @import lmerTest
#' @import ggplot2
#' @import knitr
#' @export

lmer_analysis <- function(fits, xvar, interaction=TRUE, ways=NULL){

  # fits: List. Fitted linear mixed effect model
  # xvar: string. explanatory variable interests
  # Perform Anova, find lsmeans and group differences and return the list of results

  anova <- anova(fits[[xvar]])  # anova
  lsmean <- data.frame(ls_means(fits[[xvar]]))  # find lsmean
  lsdiff <- difflsmeans(fits[[xvar]])  # find group differences

  if(interaction==TRUE){
    if(ways=='2'){
      lsmean=lsmean[grep(':',lsmean$term),]  # Only for two-way mixed effects model
      varnames=strsplit(lsmean$term[1],':')[[1]] # find two main effects variable name
      lsmean[,varnames[1]] <- unlist(lapply(lsmean$levels, function(x)strsplit(x,':')[[1]][1]))
      lsmean[,varnames[2]] <- unlist(lapply(lsmean$levels, function(x)strsplit(x,':')[[1]][2]))

      var1 <- lsmean[,varnames[1]]  # first main effect variable
      var2 <- lsmean[,varnames[2]]  # second main effect variable

      # plot
      fig <- ggplot(lsmean, aes(x=var2, y=Estimate, colour=var1, shape=var1, group=var1)) +
        geom_errorbar(aes(ymin=lower, ymax=upper), size=1, width=.2) +
        geom_line() + geom_point(size=2, show.legend=F)
    }
    if(ways=='3'){
      lsmean=lsmean[31:46,]
      varnames=strsplit(lsmean$term[1],':')[[1]]
      lsmean[,varnames[1]] <- unlist(lapply(lsmean$levels, function(x)strsplit(x,':')[[1]][1]))
      lsmean[,varnames[2]] <- unlist(lapply(lsmean$levels, function(x)strsplit(x,':')[[1]][2]))
      lsmean[,varnames[3]] <- unlist(lapply(lsmean$levels, function(x)strsplit(x,':')[[1]][3]))

      var1 <- lsmean[,varnames[1]]  # first main effect variable
      var2 <- lsmean[,varnames[2]]  # second main effect variable
      var3 <- lsmean[,varnames[3]]  # third main effect variable

      var1d <- var1[var2=='Yes']; var1dn <- var1[var2=='No']
      var2d <- var2[var2=='Yes']; var2dn <- var2[var2=='No']
      var3d <- var3[var2=='Yes']; var3dn <- var3[var2=='No']

      g1 <- ggplot(lsmean[lsmean$demented=='Yes',], aes(x=var3d, y=Estimate, colour=var1d, shape=var1d, group=var1d)) +
        geom_errorbar(aes(ymin=lower, ymax=upper), size=1, width=.2) +
        geom_line() + geom_point(size=2, show.legend=F)

      g2 <- ggplot(lsmean[lsmean$demented=='No',], aes(x=var3dn, y=Estimate, colour=var1dn, shape=var1dn, group=var1dn)) +
        geom_errorbar(aes(ymin=lower, ymax=upper), size=1, width=.2) +
        geom_line() + geom_point(size=2, show.legend=F)

      fig <- list(g1, g2)
    }}

  if(interaction==FALSE){
    time_points <- c('eval0','eval1','eval2','eval3')
    lsmean <- lsmean[time_points,]
    #plot
    fig <- ggplot(lsmean, aes(x=levels, y=Estimate)) +
      geom_errorbar(aes(ymin=lower, ymax=upper), size=.5, width=.1) +
      + geom_line() geom_point(size=1, show.legend = F)
  }

  return_list <- list(anova, lsmean, lsdiff, fig)
  names(return_list) <- c('anova', 'lsmean', 'lsdiff', 'plot')
  return(return_list)
}