stopifnot(require(knitr)) opts_chunk$set( comment=NA, message = FALSE, warning = FALSE, eval = params$EVAL, dev = "png", dpi = 500, fig.asp = 0.618, fig.width = 14, out.width = "100%", fig.align = "center" )
library("NBZIMM")
This vignette focuses on demonstrations of fitting generalized additive models for multiple longitudinal microbiome taxa with NBZIMM.
The NBZIMM provides a mgam
functions for setting up and fitting generalized additive models by working with the function gam
from R package mgcv.
In this vignette we'll use the mgam
function in the
NBZIMM package to demonstrate the analysis of a public available dataset from Romero et al (2014).
This example used a published dataset from Romero et al.(2014). It was a case-control longitudinal study collected 22 pregnant women who delivered at term and 32 non-pregnant women, each measured at multiple time points. The data consist of two data components: OTU and SampleData. OTU contains microbiome count data for 900 samples with 143 taxa; SampleData contains data of sample variables for all the samples with 9 variables, including subject ID, pregnant status, total sequencing read, age, race, etc.
Here we import the dataset from the package NBZIMM and check the clinical variables first.
data(Romero) # see help("Romero") names(Romero)
Get the necessary variables to fit the models
otu = Romero$OTU; dim(otu) sam = Romero$SampleData; dim(sam) colnames(sam) N = sam[, "Total.Read.Counts"] Days = sam$GA_Days Age = sam$Age Race = sam$Race preg = sam$pregnant; table(preg) subject = sam[, "Subect_ID"]; table(subject)
f = mgam(y=otu, formula= ~ offset(log(N)) + Days + Age + Race + preg + s(subject, bs="re"), family=nb(), min.p=0.2) out = fixed(f) para = out$parametric_terms para = para[para[,2]!="(Intercept)", ] par(mfrow = c(1, 1), cex.axis = 1, mar = c(2, 12, 4, 4)) plot.fixed(res=para[,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 = para[para$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=para, 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.