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

View source: R/mgam.R

mgamR Documentation

Fitting Mixed Models Separately for Many Responses Using gam in the package mgcv

Description

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

Usage

  
mgam(y, formula, data, family=nb(), paraPen=NULL, 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 gam.

data, family, paraPen

the same as in gam.

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 gam in the package mgcv.

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, nyi@uab.edu

See Also

gam, mms, mglmmTMB

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
Age = sam$Age
Race = sam$Race
preg = sam$pregnant; table(preg)
subject = sam[, "Subect_ID"]; table(subject)

# analyze all taxa with a given nonzero proportion

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)

nyiuab/NBZIMM documentation built on June 14, 2024, 1:45 a.m.