fitZig: Computes the weighted fold-change estimates and t-statistics.

In HCBravoLab/metagenomeSeq: Statistical analysis for sparse high-throughput sequencing

Computes the weighted fold-change estimates and t-statistics.

Description

Wrapper to actually run the Expectation-maximization algorithm and estimate $f_count$ fits. Maximum-likelihood estimates are approximated using the EM algorithm where we treat mixture membership $delta_ij = 1$ if $y_ij$ is generated from the zero point mass as latent indicator variables. The density is defined as $f_zig(y_ij = pi_j(S_j)*f_0(y_ij) +(1-pi_j (S_j)) * f_count(y_ij; mu_i, sigma_i^2)$. The log-likelihood in this extended model is: $(1-delta_ij) log f_count(y;mu_i,sigma_i^2 )+delta_ij log pi_j(s_j)+(1-delta_ij) log (1-pi_j (s_j))$. The responsibilities are defined as $z_ij = pr(delta_ij=1 | data)$.

Usage

fitZig(
  obj,
  mod,
  zeroMod = NULL,
  useCSSoffset = TRUE,
  control = zigControl(),
  useMixedModel = FALSE,
  ...
)

Arguments

obj	A MRexperiment object with count data.
mod	The model for the count distribution.
zeroMod	The zero model, the model to account for the change in the number of OTUs observed as a linear effect of the depth of coverage.
useCSSoffset	Boolean, whether to include the default scaling parameters in the model or not.
control	The settings for fitZig.
useMixedModel	Estimate the correlation between duplicate features or replicates using duplicateCorrelation.
...	Additional parameters for duplicateCorrelation.

Value

A list of objects including:

• call - the call made to fitZig

• fit - 'MLArrayLM' Limma object of the weighted fit

• countResiduals - standardized residuals of the fit

• z - matrix of the posterior probabilities

• eb - output of eBayes, moderated t-statistics, moderated F-statistics, etc

• taxa - vector of the taxa names

• counts - the original count matrix input

• zeroMod - the zero model matrix

• zeroCoef - the zero model fitted results

• stillActive - convergence

• stillActiveNLL - nll at convergence

• dupcor - correlation of duplicates

See Also

cumNorm zigControl

Examples

# This is a simple demonstration 
data(lungData)
k = grep("Extraction.Control",pData(lungData)$SampleType)
lungTrim = lungData[,-k]
k = which(rowSums(MRcounts(lungTrim)>0)<30)
lungTrim = cumNorm(lungTrim)
lungTrim = lungTrim[-k,]
smokingStatus = pData(lungTrim)$SmokingStatus
mod = model.matrix(~smokingStatus)
# The maxit is not meant to be 1 - this is for demonstration/speed
settings = zigControl(maxit=1,verbose=FALSE)
fit = fitZig(obj = lungTrim,mod=mod,control=settings)

HCBravoLab/metagenomeSeq documentation built on Dec. 20, 2024, 4:06 p.m.