mglmmTMB: Fitting Mixed Models Separately for Many Responses Using...

In nyiuab/NBZIMM: Negative Binomial and Zero-Inflated Mixed Models

Fitting Mixed Models Separately for Many Responses Using glmmTMB in the package glmmTMB

Description

This function fits mixed models separately for many responses using glmmTMB in the package glmmTMB.

Usage

  
mglmmTMB(y, formula, data, family = nbinom2(), 
        zi = ~0, disp = ~1,
        min.p = 0, verbose = TRUE)

Arguments

y	a matrix of responses; each column is a response and rows are samples.
formula	a one-sided formula of the form ~ x (i.e., the respose is omitted); the right side of ~ is the same as in glmmTMB.
data, family, zi, disp	the same as in glmmTMB.
min.p	a value in [0, 1). The responses with the proportion of non-zero values > min.p are analyzed.
verbose	logical. If TRUE, print out computational time.

Details

This function analyzes the responses in y by repeated calls to glmmTMB in the package glmmTMB.

Value

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;

Author(s)

Nengjun Yi, [e-mail]

See Also

glmmTMB, mms, mgam

Examples

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 = mglmmTMB(y = otu, formula = ~ Days + Age + Race + preg + offset(log(N))+ (1 | subject), 
        data = data, family = nbinom2, min.p = 0.2)

# zero-inflated negative binomial mixed model
f2 = mglmmTMB(y = otu, formula = ~ Days + Age + Race + preg + offset(log(N))+ (1 | subject), 
              data = data, family = nbinom2, zi = ~1,
              min.p = 0.2)

# display the results
f = f1
out = fixed(f)
out = out$cond
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)

nyiuab/NBZIMM documentation built on April 21, 2022, 7 a.m.