glmm.R

# This is the testing function that builds the models
  testMixedEffects <- function(test.col) {
    # We return this default in case of error/no fit
    default <- list(null.model=NA, model=NA)
    
    counts <- .mat[, colnames(.mat)==test.col]
    if (sum(counts > 0)/length(counts) < sparsity) {
      message(sprintf("%s: too sparse, skipping", test.col))
      return(default) 
    }
    
    # Build positive and negative counts for taxa under test
    test.data <- .mat[, c(1:5)]
    test.data$neg.cts <- test.data[["total"]] - counts
    test.data$pos.cts <- counts
    
    # The model of interest contains an interaction between StudyGroup and 
    # SampleType; the null model does not. We are stricter about forcing 
    # convergence for this model than the null model. Errors return NA.
    model1 <- plyr::try_default(
      lme4::glmer(
      formula=cbind(pos.cts, neg.cts) ~ 
        StudyGroup * SampleType + 
        (1|OriginatingSampleID) + (1|SampleID),
      data=test.data,
      family=binomial,
      control=glmerControl(
        optimizer="bobyqa", 
        check.conv.singular = "warning",
        check.conv.grad = "warning",
        check.conv.hess = "warning")
    ), NA) 

    # The null model
    model0 <- plyr::try_default(
      lme4::glmer(
      formula=cbind(pos.cts, neg.cts) ~ 
        StudyGroup + SampleType + 
        (1|OriginatingSampleID) + (1|SampleID),
      data=test.data,
      family=binomial,
      control=glmerControl(
        optimizer="bobyqa")
    ), NA)    
        
    # Return NA if the primary model broke
    if (is.na(model1)) {
      message(sprintf("Nonconvergent: %s", test.col))
      return(default)
    } else {
      # browser()
      return(list(null.model=model0, model=model1))
    }
  }