Skip to content

Instantly share code, notes, and snippets.

@eclarke eclarke/glmm.R
Created Dec 15, 2016

Embed
What would you like to do?
# 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))
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.