mms | R Documentation |
glmm.nb
, lme
, glmm.zinb
, or lme.zig
This function fits mixed models separately for many responses using glmm.nb
, lme
, glmm.zinb
, or lme.zig
.
mms(y, fixed, random, data, method = c("nb", "lme", "zinb", "zig"), correlation, zi_fixed = ~1, zi_random = NULL, niter = 30, epsilon = 1e-05, min.p = 0, verbose = TRUE)
y |
a matrix of responses; each column is a response and rows are samples. |
fixed |
a one-sided formula of the form |
data, random, correlation, zi_fixed, zi_random, niter, epsilon |
These arguments are the same as in |
method |
character specification of the method:
|
min.p |
a value in [0, 1). The responses with the proportion of non-zero values > min.p are analyzed. |
verbose |
logical. If |
This function analyzes the responses in y
by repeated calls to glmm.nb
, lme
, glmm.zinb
, or lme.zig
.
A list including fit
, responses
, variables
, and call
:
fit |
fitted models for all the responses; |
responses |
names of all the responses; |
variables |
names of all covariates; |
Nengjun Yi, nyi@uab.edu
lme
, glmm.nb
, glmm.zinb
, lme.zig
library(NBZIMM) data(Romero) names(Romero) otu = Romero$OTU; dim(otu) sam = Romero$SampleData; dim(sam) colnames(sam) N = sam[, "Total.Read.Counts"] Days = sam$GA_Days; Days = scale(Days) Age = sam$Age; Age = scale(Age) Race = sam$Race preg = sam$pregnant; table(preg) subject = sam[, "Subect_ID"]; table(subject) # analyze all taxa with a given nonzero proportion data = data.frame(Days=Days, Age=Age, Race=Race, preg=preg, N=N, subject=subject) # negative binomial mixed model f1 = mms(y = otu, fixed = ~ Days + Age + Race + preg + offset(log(N)), random = ~ 1 | subject, data = data, min.p = 0.2, method = "nb") # zero-inflated negative binomial mixed model f2 = mms(y = otu, fixed = ~ Days + Age + Race + preg + offset(log(N)), random = ~ 1 | subject, zi_fixed = ~1, data = data, min.p = 0.2, method = "zinb") # zero-inflated gaussian mixed model f3 = mms(y = log2(otu+1), fixed = ~ Days + Age + Race + preg + offset(log(N)), random = ~ 1 | subject, zi_fixed = ~1, data = data, min.p = 0.2, method = "zig") # display the results f = f1 out = fixed(f) out = out$dist out = out[out[,2]!="(Intercept)", ] par(mfrow = c(1, 1), cex.axis = 1, mar = c(2, 12, 4, 4)) plot.fixed(res=out[,c("Estimate","Std.Error","padj")], threshold=0.001, gap=500, col.pts=c("black", "grey"), cex.axis=0.6, cex.var=0.7) preg.coefs = out[out$variables=="preg",] par(mfrow = c(1, 1), cex.axis = 1, mar = c(2, 10, 4, 4)) plot.fixed(res=preg.coefs[,c("Estimate","Std.Error","padj")], threshold=0.05, gap=300, main="Covariate: pregnant", cex.axis=0.6, cex.var=0.7) g = heat.p(df=out, p.breaks = c(0.001, 0.01, 0.05), colors = c("black", "darkgrey", "grey", "lightgrey"), zigzag=c(T,F), abbrv=c(T,F), margin=c(1,6), y.size=8, legend=T)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.