fido: Fitting and Analysis of Multinomial Logistic Normal Models

Description

Provides methods for fitting and inspection of Bayesian Multinomial Logistic Normal Models using MAP estimation and Laplace Approximation. Key functionality is implemented in C++ for scalability.

Author(s)

Maintainer: Justin Silverman Justin.Silverman@psu.edu

Other contributors:

• Kim Roche kimberly.roche@duke.edu [contributor]

• Michelle Nixon pistner@psu.edu [contributor]

See Also

Useful links:

• https://jsilve24.github.io/fido/

• Report bugs at https://github.com/jsilve24/fido/issues

Data from Gevers et al. (2014)

Description

OTU data and metadata for 1,359 samples in a Crohn's disease study

Usage

data(RISK_CCFA)

Format

An otu table, sample data table, and taxonomy table.

Details

Study is described here: https://pubmed.ncbi.nlm.nih.gov/24629344/. Data was obtained from https://github.com/twbattaglia/MicrobeDS.

References

Gevers D, et al. The treatment-naive microbiome in new-onset Crohn's disease. Cell Host Microbe. 2014 Mar 12;15(3):382-392. doi: 10.1016/j.chom.2014.02.005. PMID: 24629344; PMCID: PMC4059512.

Data from Gevers et al. (2014)

Description

OTU data and metadata for 1,359 samples in a Crohn's disease study

Usage

data(RISK_CCFA)

Format

A matrix otu table.

Details

Study is described here: https://pubmed.ncbi.nlm.nih.gov/24629344/. Data was obtained from https://github.com/twbattaglia/MicrobeDS.

References

Gevers D, et al. The treatment-naive microbiome in new-onset Crohn's disease. Cell Host Microbe. 2014 Mar 12;15(3):382-392. doi: 10.1016/j.chom.2014.02.005. PMID: 24629344; PMCID: PMC4059512.

Data from Gevers et al. (2014)

Description

OTU data and metadata for 1,359 samples in a Crohn's disease study

Usage

data(RISK_CCFA)

Format

A sample data table.

Details

Study is described here: https://pubmed.ncbi.nlm.nih.gov/24629344/. Data was obtained from https://github.com/twbattaglia/MicrobeDS.

References

Gevers D, et al. The treatment-naive microbiome in new-onset Crohn's disease. Cell Host Microbe. 2014 Mar 12;15(3):382-392. doi: 10.1016/j.chom.2014.02.005. PMID: 24629344; PMCID: PMC4059512.

Data from Gevers et al. (2014)

Description

OTU data and metadata for 1,359 samples in a Crohn's disease study

Usage

data(RISK_CCFA)

Format

A taxonomy table.

Details

Study is described here: https://pubmed.ncbi.nlm.nih.gov/24629344/. Data was obtained from https://github.com/twbattaglia/MicrobeDS.

References

Gevers D, et al. The treatment-naive microbiome in new-onset Crohn's disease. Cell Host Microbe. 2014 Mar 12;15(3):382-392. doi: 10.1016/j.chom.2014.02.005. PMID: 24629344; PMCID: PMC4059512.

Data from Silverman et al. (2019) bioRxiv

Description

Mock communities and calibration samples created for measuring and validating model of PCR bias. This data has been preprocessed as in the original manuscript.

Format

an matrix Y (counts for each community member)

References

Justin D. Silverman, Rachael J. Bloom, Sharon Jiang, Heather K. Durand, Sayan Mukherjee, Lawrence A. David. (2019) Measuring and Mitigating PCR Bias in Microbiome Data. bioRxiv 604025; doi: https://doi.org/10.1101/604025

Compute the ALR of a matrix

Description

Compute the ALR of a matrix

Usage

alr(x, d = NULL)

Arguments

Value

matrix

Compute the inverse ALR of a matrix

Description

Compute the inverse ALR of a matrix

Usage

alrInv(y, d = NULL)

Arguments

Value

matrix

Compute the ALR of an array

Description

Compute the ALR of an array

Usage

alrInv_array(y, d = dim(y)[coords] + 1, coords)

Arguments

Value

array

Compute the ALR of an array

Description

Compute the ALR of an array

Usage

alr_array(x, d = dim(x)[parts], parts)

Arguments

Value

array

Convert object of class orthusfit to a list

Description

Convert object of class orthusfit to a list

Usage

## S3 method for class 'orthusfit'
as.list(x, ...)

Arguments

Value

A list of the converted orthusfit object

Convert object of class pibblefit to a list

Description

Convert object of class pibblefit to a list

Usage

## S3 method for class 'pibblefit'
as.list(x, ...)

Arguments

Value

A list from the converted pibblefit object.

convert list to orthusfit

Description

convert list to orthusfit

Usage

as.orthusfit(object)

Arguments

Value

An orthusfit object

convert list to pibblefit

Description

convert list to pibblefit

Usage

as.pibblefit(object)

Arguments

Value

A pibblefit object

Interface to fit basset models

Description

Basset (A Lazy Learner) - non-linear regression models in fido

Usage

basset(
  Y = NULL,
  X,
  upsilon = NULL,
  Theta = NULL,
  Gamma = NULL,
  Xi = NULL,
  linear = NULL,
  init = NULL,
  pars = c("Eta", "Lambda", "Sigma"),
  newdata = NULL,
  ...
)

## S3 method for class 'bassetfit'
refit(m, pars = c("Eta", "Lambda", "Sigma"), ...)

Arguments

Details

the full model is given by:

Y_j \sim Multinomial(\pi_j)

\pi_j = \Phi^{-1}(\eta_j)

\eta \sim MN_{D-1 \times N}(\Lambda, \Sigma, I_N)

\Lambda \sim GP_{D-1 \times Q}(\Theta(X), \Sigma, \Gamma(X))

\Sigma \sim InvWish(\upsilon, \Xi)

Where \Gamma(X) is short hand for the Gram matrix of the Kernel function.

Alternatively can be used to fit an additive GP of the form:

Y_j \sim Multinomial(\pi_j)

\pi_j = \Phi^{-1}(\eta_j)

\eta \sim MN_{D-1 \times N}(\Lambda, \Sigma, I_N)

\Lambda = \Lambda_1 + ... + \Lambda_p + B X

\Lambda_1 \sim GP_{D-1 \times Q}(\Theta_1(X), \Sigma, \Gamma_1(X))

...

\Lambda_p \sim GP_{D-1 \times Q}(\Theta_p(X), \Sigma, \Gamma_p(X))

B \sim MN(\Theta_B, \Sigma, \Gamma_B)

\Sigma \sim InvWish(\upsilon, \Xi)

Where \Gamma(X) is short hand for the Gram matrix of the Kernel function.

Default behavior is to use MAP estimate for uncollaping the LTP model if laplace approximation is not preformed.

Value

an object of class bassetfit

Check vector/matrix/data.frame for expected dimensions or throw error

Description

Check vector/matrix/data.frame for expected dimensions or throw error

Usage

check_dims(x, d, par)

Arguments

Value

nothing if no error, otherwise throws error

Examples

y <- c(1,3,4)
check_dims(y, 3, "y")

Compute the CLR of an array

Description

Compute the CLR of an array

Usage

clr_array(x, parts)

Arguments

Value

array

Return regression coefficients of orthus object

Description

Orthus: Returned as array of dimension (D-1+P) x Q x iter (if in ALR or ILR) otherwise (D+P) x Q x iter.

Usage

## S3 method for class 'orthusfit'
coef(object, ...)

Arguments

Details

Other arguments:

• use_names if column and row names were passed for Y and X in call to pibble, should these names be applied to output array.

Value

Array of dimension (D-1) x Q x iter

Return regression coefficients of pibblefit object

Description

Pibble: Returned as array of dimension (D-1) x Q x iter (if in ALR or ILR) otherwise DxQxiter (if in proportions or clr).

Usage

## S3 method for class 'pibblefit'
coef(object, ...)

Arguments

Details

Other arguments:

• 'use_names' if column and row names were passed for Y and X in call to pibble, should these names be applied to output array.

Value

Array of dimension (D-1) x Q x iter

Solve Bayesian Multivariate Conjugate Linear Model

Description

See details for model. Notation: N is number of samples, D is the dimension of the response, Q is number of covariates.

Usage

conjugateLinearModel(Y, X, Theta, Gamma, Xi, upsilon, n_samples = 2000L)

Arguments

Details

Y \sim MN_{D \times N}(\Lambda \mathbf{X}, \Sigma, I_N)

\Lambda \sim MN_{D \times Q}(\Theta, \Sigma, \Gamma)

\Sigma \sim InvWish(\upsilon, \Xi)

This function provides a means of sampling from the posterior distribution of Lambda and Sigma.

Value

List with components

1. Lambda Array of dimension D x Q x n_samples (posterior samples)

2. Sigma Array of dimension D x D x n_samples (posterior samples)

Examples

sim <- pibble_sim()
eta.hat <- t(alr(t(sim$Y+0.65)))
fit <- conjugateLinearModel(eta.hat, sim$X, sim$Theta, sim$Gamma, 
                            sim$Xi, sim$upsilon, n_samples=2000)

Convert orthus covariance matricies between representations

Description

Convert orthus covariance matricies between representations

Usage

oilrvar2ilrvar(Sigma, s, V1, V2)

oilrvar2clrvar(Sigma, s, V)

oclrvar2ilrvar(Sigma, s, V)

oalrvar2clrvar(Sigma, s, d1)

oclrvar2alrvar(Sigma, s, d2)

oalrvar2alrvar(Sigma, s, d1, d2)

oalrvar2ilrvar(Sigma, s, d1, V2)

oilrvar2alrvar(Sigma, s, V1, d2)

Arguments

Value

matrix

Create a default ILR base

Description

Create a default ILR base

Usage

create_default_ilr_base(D)

Arguments

Value

A matrix

Transform Fit fido Parameters to other representations

Description

These are a collection of convenience functions for transforming fido fit objects to a number of different representations including ILR bases, CLR coordinates, ALR coordinates, and proportions.

Usage

to_proportions(m)

to_alr(m, d)

to_ilr(m, V = NULL)

to_clr(m)

## S3 method for class 'pibblefit'
to_proportions(m)

## S3 method for class 'orthusfit'
to_proportions(m)

## S3 method for class 'pibblefit'
to_alr(m, d)

## S3 method for class 'orthusfit'
to_alr(m, d)

## S3 method for class 'pibblefit'
to_ilr(m, V = NULL)

## S3 method for class 'orthusfit'
to_ilr(m, V = NULL)

## S3 method for class 'pibblefit'
to_clr(m)

## S3 method for class 'orthusfit'
to_clr(m)

Arguments

Details

For orthus, transforms only appleid to log-ratio parameters

Note: that there is a degeneracy of representations for a covariance matrix represented in terms of proportions. As such the function to_proportions does not attempt to transform parameters Sigma or prior Xi and instead just removes them from the pibblefit object returned.

Value

object

Gather Multidimensional Array to Tidy Tibble

Description

Gather Multidimensional Array to Tidy Tibble

Usage

gather_array(a, value, ..., .id = NULL)

Arguments

Value

data.frame

See Also

spread_array

Examples

a <- array(1:100, dim =c(10, 5, 2))
gather_array(a, sequence, A, B, C)

Multivariate RBF Kernel

Description

Designed to be partially specified. (see examples)

Usage

SE(X, sigma = 1, rho = median(as.matrix(dist(t(X)))), jitter = 1e-10)

LINEAR(X, sigma = 1, c = rep(0, nrow(X)))

Arguments

Details

Gram matrix G is given by

SE (squared exponential):

G = \sigma^2 * exp(-[(X-c)'(X-c)]/(s*\rho^2))

LINEAR:

G = \sigma^2*(X-c)'(X-c)

Value

Gram Matrix (N x N) (e.g., the Kernel evaluated at each pair of points)

Transform Lambda into IQLR (Inter-Quantile Log-Ratio)

Description

Takes idea from Wu et al. (citation below) and calculates IQLR for Lambda, potentially useful if you believe there is an invariant group of categories (e.g., taxa / genes) that are not changing (in absolute abundance) between samples. IQLR is defined as

IQLR_x = log(x_i/g(IQVF))

for i in 1,...,D. IQVF are the CLR coordinates whose variance is within the inter-quantile range (defined by probs argument to this function). A different IQVF is fit for each posteior sample as the IQVFs are calculted based on posterior estimates for Lambda. The variance of a CLR coordinate is defined as the norm of each row of Lambda[,focus.cov] (i.e., the covariation in Eta, explained by those covariates). This definition of variance allows uses to exclude variation from technical / trivial sources in calculation of IQVF/IQLR.

Usage

lambda_to_iqlr(m, focus.cov = NULL, probs = c(0.25, 0.75))

Arguments

Details

Primarily intended for doing differential expression analysis under assumption that only small group of categories (e.g., taxa / genes) are changing

Value

array of dimension (D, Q, iter) where D is number of taxa, Q is number of covariates, and iter is number of posterior samples.

References

Jia R. Wu, Jean M. Macklaim, Briana L. Genge, Gregory B. Gloor (2017) Finding the center: corrections for asymmetry in high-throughput sequencing datasets. arxiv:1704.01841v1

Examples

sim <- pibble_sim()
fit <- pibble(sim$Y, sim$X)
# Use first two covariates to define iqlr, just show first 5 samples
lambda_to_iqlr(fit, 1:2)[,,1:5] 

Log of Multivarate Gamma Function - Gamma_p(a)

Description

Log of Multivarate Gamma Function - Gamma_p(a)

Usage

lmvgamma(a, p)

Arguments

Value

Numeric

References

https://en.wikipedia.org/wiki/Multivariate_gamma_function

Derivative of Log of Multivariate Gamma Function - Gamma_p(a)

Description

Derivative of Log of Multivariate Gamma Function - Gamma_p(a)

Usage

lmvgamma_deriv(a, p)

Arguments

Value

Numeric

References

https://en.wikipedia.org/wiki/Multivariate_gamma_function

Calculations for the Collapsed Pibble Model

Description

Functions providing access to the Log Likelihood, Gradient, and Hessian of the collapsed pibble model. Note: These are convenience functions but are not as optimized as direct coding of the PibbleCollapsed C++ class due to a lack of Memoization. By contrast function optimPibbleCollapsed is much more optimized and massively cuts down on repeated calculations. A more efficient Rcpp module based implementation of these functions may following if the future. For model details see optimPibbleCollapsed documentation

Usage

loglikPibbleCollapsed(Y, upsilon, ThetaX, KInv, AInv, eta, sylv = FALSE)

gradPibbleCollapsed(Y, upsilon, ThetaX, KInv, AInv, eta, sylv = FALSE)

hessPibbleCollapsed(Y, upsilon, ThetaX, KInv, AInv, eta, sylv = FALSE)

Arguments

Value

see below

• loglikPibbleCollapsed - double

• gradPibbleCollapsed - vector

• hessPibbleCollapsed- matrix

Examples

D <- 10
Q <- 2
N <- 30

# Simulate Data
Sigma <- diag(sample(1:8, D-1, replace=TRUE))
Sigma[2, 3] <- Sigma[3,2] <- -1
Gamma <- diag(sqrt(rnorm(Q)^2))
Theta <- matrix(0, D-1, Q)
Phi <-  Theta + t(chol(Sigma))%*%matrix(rnorm(Q*(D-1)), nrow=D-1)%*%chol(Gamma)
X <- matrix(rnorm(N*(Q-1)), Q-1, N)
X <- rbind(1, X)
Eta <- Phi%*%X + t(chol(Sigma))%*%matrix(rnorm(N*(D-1)), nrow=D-1)
Pi <- t(alrInv(t(Eta)))
Y <- matrix(0, D, N)
for (i in 1:N) Y[,i] <- rmultinom(1, sample(5000:10000), prob = Pi[,i])

# Priors
upsilon <- D+10
Xi <- Sigma*(upsilon-D)

# Precompute
KInv <- solve(Xi)
AInv <- solve(diag(N)+ t(X)%*%Gamma%*%X)
ThetaX <- Theta%*%X


loglikPibbleCollapsed(Y, upsilon, ThetaX, KInv, AInv, Eta)
gradPibbleCollapsed(Y, upsilon, ThetaX, KInv, AInv, Eta)[1:5]
hessPibbleCollapsed(Y, upsilon, ThetaX, KInv, AInv, Eta)[1:5,1:5]

Data from Silverman et al. (2018) Microbiome

Description

High Resolution (hourly and daily) sampling of 4 in vitro artificial gut models with many technical replicates to identify technical variation.

Usage

data(mallard)

Format

A list containing "otu_table", "sample_data", "tax_table", and "refseq".

Details

This data is at the sequence variant level. Data at the family level processed as in Silverman et al. 2018 is given in mallard_family

References

Silverman et al. "Dynamic linear models guide design and analysis of microbiota studies within artificial human guts". Microbiome 2018 6:202

Data from Silverman et al. (2018) Microbiome

Description

High Resolution (hourly and daily) sampling of 4 in vitro artificial gut models with many technical replicates to identify technical variation.

Usage

data(mallard_family)

Format

A list containing "otu_table", "sample_data", "tax_table", and "refseq".

Details

This data is at the family level and processed as in Silverman et al. 2018. Data at the sequence variant level without preprocessing is given in mallard

References

Silverman et al. "Dynamic linear models guide design and analysis of microbiota studies within artificial human guts". Microbiome 2018 6:202

Data from Silverman et al. (2019) bioRxiv

Description

Mock communities and calibration samples created for measuring and validating model of PCR bias. This data has been preprocessed as in the original manuscript.

Format

a data.frame metadata associated with the counts matrix 'Y'

References

Justin D. Silverman, Rachael J. Bloom, Sharon Jiang, Heather K. Durand, Sayan Mukherjee, Lawrence A. David. (2019) Measuring and Mitigating PCR Bias in Microbiome Data. bioRxiv 604025; doi: https://doi.org/10.1101/604025

Closure operator

Description

Closure operator

Usage

miniclo(x)

Arguments

Value

x with row entries divided by sum of row (converts vectors to row matricies)

Examples

x <- matrix(runif(30), 10, 3)
x <- miniclo(x)

Closure Operation applied to array on margin

Description

Array version of miniclo.

Usage

miniclo_array(x, parts)

Arguments

Value

array

Examples

x <- array(1:100, dim=c(10, 5, 2))
miniclo_array(x, parts=2)

mongrel

Description

This function is deprecated, please use pibble instead.

Usage

mongrel(
  Y = NULL,
  X = NULL,
  upsilon = NULL,
  Theta = NULL,
  Gamma = NULL,
  Xi = NULL,
  init = NULL,
  pars = c("Eta", "Lambda", "Sigma"),
  ...
)

Arguments

Value

An object of class pibblefit

Generic method for applying names to an object

Description

Intended to be called internally by package

Usage

name(m, ...)

Arguments

Value

object of same class but with names applied to dimensions

S3 for orthusfit apply names to orthusfit object

Description

To avoid confusion, assigned default names to multinomial categories (c1 etc...) and zdimensions (z1 etc...)

Usage

## S3 method for class 'orthusfit'
name(m, ...)

Arguments

Value

object of class orthusfit

S3 for pibblefit apply names to pibblefit object

Description

S3 for pibblefit apply names to pibblefit object

Usage

## S3 method for class 'pibblefit'
name(m, ...)

Arguments

Value

object of class pibblefit

Generic method for getting and setting dimension names of fit object

Description

Generic method for getting and setting dimension names of fit object

Usage

## S3 method for class 'pibblefit'
names_covariates(m)

## S3 method for class 'pibblefit'
names_samples(m)

## S3 method for class 'pibblefit'
names_categories(m)

## S3 method for class 'pibblefit'
names_coords(m)

## S3 replacement method for class 'pibblefit'
names_covariates(m) <- value

## S3 replacement method for class 'pibblefit'
names_samples(m) <- value

## S3 replacement method for class 'pibblefit'
names_categories(m) <- value

names_covariates(m)

names_samples(m)

names_categories(m)

names_coords(m)

names_covariates(m) <- value

names_samples(m) <- value

names_categories(m) <- value

Arguments

Details

names_coords is different than names_categories. names_categories provides access to the basic names of each multinomial category. In contrast, names_coords provides access to the names of the coordinates in which an object is represented. These coordinate names are based on the category names. For example, category names may be, (OTU1, ..., OTUD) where as coordinate names could be (log(OTU1/OTUD), etc...) if object is in default coordinate system.

Value

A vector of names

Generic method for accessing model fit dimensions

Description

Generic method for accessing model fit dimensions

Usage

## S3 method for class 'pibblefit'
ncategories(m)

## S3 method for class 'pibblefit'
nsamples(m)

## S3 method for class 'pibblefit'
ncovariates(m)

## S3 method for class 'pibblefit'
niter(m)

## S3 method for class 'orthusfit'
ncategories(m)

## S3 method for class 'orthusfit'
nsamples(m)

## S3 method for class 'orthusfit'
ncovariates(m)

## S3 method for class 'orthusfit'
niter(m)

ncategories(m)

nsamples(m)

ncovariates(m)

niter(m)

Arguments

Details

An alternative approach to accessing these dimensions is to access them directly from the pibblefit object using list indexing. * ncategories is equivalent to m$D * nsamples is equivalent to m$N * ncovariates is equivalent to m$Q

Value

integer

Function to Optimize the Collapsed Pibble Model

Description

See details for model. Should likely be followed by function uncollapsePibble. Notation: N is number of samples, D is number of multinomial categories, and Q is number of covariates.

Usage

optimPibbleCollapsed(
  Y,
  upsilon,
  ThetaX,
  KInv,
  AInv,
  init,
  n_samples = 2000L,
  calcGradHess = TRUE,
  b1 = 0.9,
  b2 = 0.99,
  step_size = 0.003,
  epsilon = 1e-06,
  eps_f = 1e-10,
  eps_g = 1e-04,
  max_iter = 10000L,
  verbose = FALSE,
  verbose_rate = 10L,
  decomp_method = "cholesky",
  optim_method = "lbfgs",
  eigvalthresh = 0,
  jitter = 0,
  multDirichletBoot = -1,
  useSylv = TRUE,
  ncores = -1L,
  seed = -1L
)

Arguments

Details

Notation: Let Z_j denote the J-th row of a matrix Z. Model:

Y_j \sim Multinomial(\pi_j)

\pi_j = \Phi^{-1}(\eta_j)

\eta \sim T_{D-1, N}(\upsilon, \Theta X, K, A)

Where A = I_N + X \Gamma X', K is a (D-1)x(D-1) covariance matrix, \Gamma is a Q x Q covariance matrix, and \Phi^{-1} is ALRInv_D transform.

Gradient and Hessian calculations are fast as they are computed using closed form solutions. That said, the Hessian matrix can be quite large [N*(D-1) x N*(D-1)] and storage may be an issue.

Note: Warnings about large negative eigenvalues can either signal that the optimizer did not reach an optima or (more commonly in my experience) that the prior / degrees of freedom for the covariance (given by parameters upsilon and KInv) were too specific and at odds with the observed data. If you get this warning try the following.

1. Try restarting the optimization using a different initial guess for eta

2. Try decreasing (or even increasing )step_size (by increments of 0.001 or 0.002) and increasing max_iter parameters in optimizer. Also can try increasing b1 to 0.99 and decreasing eps_f by a few orders of magnitude

3. Try relaxing prior assumptions regarding covariance matrix. (e.g., may want to consider decreasing parameter upsilon closer to a minimum value of D)

4. Try adding small amount of jitter (e.g., set jitter=1e-5) to address potential floating point errors.

Value

List containing (all with respect to found optima)

1. LogLik - Log Likelihood of collapsed model (up to proportionality constant)

2. Gradient - (if calcGradHess=true)

3. Hessian - (if calcGradHess=true) of the POSITIVE LOG POSTERIOR

4. Pars - Parameter value of eta at optima

5. Samples - (D-1) x N x n_samples array containing posterior samples of eta based on Laplace approximation (if n_samples>0)

6. Timer - Vector of Execution Times

7. logInvNegHessDet - the log determinant of the covariacne of the Laplace approximation, useful for calculating marginal likelihood

8. logMarginalLikelihood - A calculation of the log marginal likelihood based on the laplace approximation

References

S. Ruder (2016) An overview of gradient descent optimization algorithms. arXiv 1609.04747

JD Silverman K Roche, ZC Holmes, LA David, S Mukherjee. Bayesian Multinomial Logistic Normal Models through Marginally Latent Matrix-T Processes. 2022, Journal of Machine Learning

See Also

uncollapsePibble

Examples

sim <- pibble_sim()

# Fit model for eta
fit <- optimPibbleCollapsed(sim$Y, sim$upsilon, sim$Theta%*%sim$X, sim$KInv, 
                             sim$AInv, random_pibble_init(sim$Y))  

Interface to fit orthus models

Description

This function is largely a more user friendly wrapper around optimPibbleCollapsed and uncollapsePibble for fitting orthus models. See details for model specification. Notation: N is number of samples, P is the number of dimensions of observations in the second dataset, D is number of multinomial categories, Q is number of covariates, iter is the number of samples of eta (e.g., the parameter n_samples in the function optimPibbleCollapsed)

Usage

orthus(
  Y = NULL,
  Z = NULL,
  X = NULL,
  upsilon = NULL,
  Theta = NULL,
  Gamma = NULL,
  Xi = NULL,
  init = NULL,
  pars = c("Eta", "Lambda", "Sigma"),
  ...
)

Arguments

Details

the full model is given by:

Y_j \sim Multinomial(\pi_j)

\pi_j = \Phi^{-1}(\eta_j)

cbind(\eta, Z) \sim MN_{D-1+P \times N}(\Lambda X, \Sigma, I_N)

\Lambda \sim MN_{D-1+P \times Q}(\Theta, \Sigma, \Gamma)

\Sigma \sim InvWish(\upsilon, \Xi)

Where \Gamma is a Q x Q covariance matrix, and \Phi^{-1} is ALRInv_D transform. That is, the orthus model models the latent multinomial log-ratios (Eta) and the observations of the second dataset jointly as a linear model. This allows Sigma to also describe the covariation between the two datasets.

Default behavior is to use MAP estimate for uncollaping the LTP model if laplace approximation is not preformed.

Value

an object of class pibblefit

References

JD Silverman K Roche, ZC Holmes, LA David, S Mukherjee. Bayesian Multinomial Logistic Normal Models through Marginally Latent Matrix-T Processes. 2019, arXiv e-prints, arXiv:1903.11695

See Also

fido_transforms provide convenience methods for transforming the representation of pibblefit objects (e.g., conversion to proportions, alr, clr, or ilr coordinates.)

access_dims provides convenience methods for accessing dimensions of pibblefit object

Examples

sim <- orthus_sim()
fit <- orthus(sim$Y, sim$Z, sim$X)

Log-Ratio transforms for orthus objects

Description

Log-Ratio transforms for orthus objects

Usage

oglr(x, s, V)

oglrInv(x, s, V)

oalr(x, s, d = NULL)

oalrInv(y, s, d = NULL)

oilr(x, s, V = NULL)

oilrInv(y, s, V = NULL)

oclr(x, s)

oclrInv(x, s)

Arguments

Value

A matrix

Simulate simple orthus dataset and priors (for testing)

Description

Simulate simple orthus dataset and priors (for testing)

Usage

orthus_sim(
  D = 10,
  P = 10,
  N = 30,
  Q = 2,
  use_names = TRUE,
  true_priors = FALSE
)

Arguments

Value

list

Examples

sim <- orthus_sim()

Convert orthus samples of Eta Lambda and Sigma to tidy format

Description

Combines them all into a single tibble, see example for formatting and column headers. Primarily designed to be used by summary.orthusfit.

Usage

orthus_tidy_samples(m, use_names = FALSE, as_factor = FALSE)

Arguments

Value

tibble

Examples

sim <- orthus_sim()
fit <- orthus(sim$Y, sim$Z, sim$X)
fit_tidy <- orthus_tidy_samples(fit, use_names=TRUE)
head(fit_tidy)

Create orthusfit object

Description

Create orthusfit object

Usage

orthusfit(
  D,
  N,
  Q,
  P,
  coord_system,
  iter = NULL,
  alr_base = NULL,
  ilr_base = NULL,
  Eta = NULL,
  Lambda = NULL,
  Sigma = NULL,
  Sigma_default = NULL,
  Z = NULL,
  Y = NULL,
  X = NULL,
  upsilon = NULL,
  Theta = NULL,
  Xi = NULL,
  Xi_default = NULL,
  Gamma = NULL,
  init = NULL,
  names_categories = NULL,
  names_samples = NULL,
  names_Zdimensions = NULL,
  names_covariates = NULL
)

Arguments

Value

object of class orthusfit

See Also

pibble

Data from Silverman et al. (2019) bioRxiv

Description

Mock communities and calibration samples created for measuring and validating model of PCR bias. This data has been preprocessed as in the original manuscript.

Usage

data(pcrbias_mock)

Format

an matrix Y (counts for each community member) and a data.frame metadata

References

Justin D. Silverman, Rachael J. Bloom, Sharon Jiang, Heather K. Durand, Sayan Mukherjee, Lawrence A. David. (2019) Measuring and Mitigating PCR Bias in Microbiome Data. bioRxiv 604025; doi: https://doi.org/10.1101/604025

Interface to fit pibble models

Description

This function is largely a more user friendly wrapper around optimPibbleCollapsed and uncollapsePibble. See details for model specification. Notation: N is number of samples, D is number of multinomial categories, Q is number of covariates, iter is the number of samples of eta (e.g., the parameter n_samples in the function optimPibbleCollapsed)

Usage

pibble(
  Y = NULL,
  X = NULL,
  upsilon = NULL,
  Theta = NULL,
  Gamma = NULL,
  Xi = NULL,
  init = NULL,
  pars = c("Eta", "Lambda", "Sigma"),
  newdata = NULL,
  ...
)

## S3 method for class 'pibblefit'
refit(m, pars = c("Eta", "Lambda", "Sigma"), ...)

Arguments

Details

the full model is given by:

Y_j \sim Multinomial(\pi_j)

\pi_j = \Phi^{-1}(\eta_j)

\eta \sim MN_{D-1 \times N}(\Lambda X, \Sigma, I_N)

\Lambda \sim MN_{D-1 \times Q}(\Theta, \Sigma, \Gamma)

\Sigma \sim InvWish(\upsilon, \Xi)

Where \Gamma is a Q x Q covariance matrix, and \Phi^{-1} is ALRInv_D transform.

Default behavior is to use MAP estimate for uncollaping the LTP model if laplace approximation is not preformed.

Value

an object of class pibblefit

References

JD Silverman K Roche, ZC Holmes, LA David, S Mukherjee. Bayesian Multinomial Logistic Normal Models through Marginally Latent Matrix-T Processes. 2019, arXiv e-prints, arXiv:1903.11695

See Also

fido_transforms provide convenience methods for transforming the representation of pibblefit objects (e.g., conversion to proportions, alr, clr, or ilr coordinates.)

access_dims provides convenience methods for accessing dimensions of pibblefit object

Generic functions including summary, print, coef, as.list, predict, name, and sample_prior name_dims

Plotting functions provided by plot and ppc (posterior predictive checks)

Examples

sim <- pibble_sim()
fit <- pibble(sim$Y, sim$X)

Simulate simple pibble dataset and priors (for testing)

Description

Simulate simple pibble dataset and priors (for testing)

Usage

pibble_sim(D = 10, N = 30, Q = 2, use_names = TRUE, true_priors = FALSE)

Arguments

Value

list

Examples

sim <- pibble_sim()

Convert pibble samples of Eta Lambda and Sigma to tidy format

Description

Combines them all into a single tibble, see example for formatting and column headers. Primarily designed to be used by summary.pibblefit.

Usage

pibble_tidy_samples(m, use_names = FALSE, as_factor = FALSE)

Arguments

Value

tibble

Examples

sim <- pibble_sim()
fit <- pibble(sim$Y, sim$X)
fit_tidy <- pibble_tidy_samples(fit, use_names=TRUE)
head(fit_tidy)

Create pibblefit object

Description

Create pibblefit object

Usage

pibblefit(
  D,
  N,
  Q,
  coord_system,
  iter = NULL,
  alr_base = NULL,
  ilr_base = NULL,
  Eta = NULL,
  Lambda = NULL,
  Sigma = NULL,
  Sigma_default = NULL,
  Y = NULL,
  X = NULL,
  upsilon = NULL,
  Theta = NULL,
  Xi = NULL,
  Xi_default = NULL,
  Gamma = NULL,
  init = NULL,
  names_categories = NULL,
  names_samples = NULL,
  names_covariates = NULL
)

Arguments

Value

object of class pibblefit

See Also

pibble

Plot Summaries of Posterior Distribution of pibblefit Parameters

Description

Plot Summaries of Posterior Distribution of pibblefit Parameters

Usage

## S3 method for class 'pibblefit'
plot(x, ...)

Arguments

Details

Other arguments:

• 'par' parameter to plot (options: Lambda, Eta, and Sigma) (default="Lambda")

• 'focus.cov' vector of covariates to include in plot (plots all if NULL)

• 'focus.coord' vector of coordinates to include in plot (plots all if NULL)

• 'focus.sample' vector of samples to include in plot (plots all if NULL)

• 'use_names' if TRUE, uses dimension names found in data as plot labels rather than using dimension integer indices.

Value

ggplot object

Examples

sim <- pibble_sim(N=10, D=4, Q=3)
fit <- pibble(sim$Y, sim$X)
plot(fit, par="Lambda")
plot(fit, par="Sigma")

Generic method for visualizing posterior predictive checks

Description

Generic method for visualizing posterior predictive checks

Usage

ppc(m, ...)

Arguments

Value

A plot

Visualization of Posterior Predictive Check of fit model

Description

Visualization of Posterior Predictive Check of fit model

Usage

## S3 method for class 'pibblefit'
ppc(m, ...)

Arguments

Details

ppc.pibblefit accepts the following additional arguments:

• "type" type of plot (options "lines", "points", "bounds")

• "iter" number of samples from posterior predictive distribution to plot (currently must be <= m$iter) if type=="lines" default is 50, if type=="ribbon" default is to use all available iterations.

• "from_scratch" should predictions of Y come from fitted Eta or from predictions of Eta from posterior of Lambda? (default: false)

Value

ggplot object

Examples

sim <- pibble_sim()
fit <- pibble(sim$Y, sim$X)
ppc(fit)

Generic Method to Plot Posterior Predictive Summaries

Description

Generic Method to Plot Posterior Predictive Summaries

Usage

## S3 method for class 'pibblefit'
ppc_summary(m, from_scratch = FALSE, ...)

ppc_summary(m, ...)

Arguments

Value

vector

Predict using basset

Description

Predict using basset

Usage

## S3 method for class 'bassetfit'
predict(
  object,
  newdata = NULL,
  response = "Lambda",
  size = NULL,
  use_names = TRUE,
  summary = FALSE,
  iter = NULL,
  from_scratch = FALSE,
  ...
)

Arguments

Details

currently only implemented for pibblefit objects in coord_system "default" "alr", or "ilr".

Value

(if summary==FALSE) array D x N x iter; (if summary==TRUE) tibble with calculated posterior summaries

Predict response from new data

Description

Predict response from new data

Usage

## S3 method for class 'pibblefit'
predict(
  object,
  newdata = NULL,
  response = "LambdaX",
  size = NULL,
  use_names = TRUE,
  summary = FALSE,
  iter = NULL,
  from_scratch = FALSE,
  ...
)

Arguments

Details

currently only implemented for pibblefit objects in coord_system "default" "alr", or "ilr".

Value

(if summary==FALSE) array D x N x iter; (if summary==TRUE) tibble with calculated posterior summaries

Examples

sim <- pibble_sim()
fit <- pibble(sim$Y, sim$X)
predict(fit)[,,1:2] # just show 2 samples

Print dimensions and coordinate system information for orthusfit object.

Description

Print dimensions and coordinate system information for orthusfit object.

Usage

## S3 method for class 'orthusfit'
print(x, summary = FALSE, ...)

Arguments

Value

No direct return, prints out summary

See Also

summary.orthusfit summarizes posterior intervals

Examples

sim <- orthus_sim()
fit <- orthus(sim$Y, sim$Z, sim$X)
print(fit)

Print dimensions and coordinate system information for pibblefit object.

Description

Print dimensions and coordinate system information for pibblefit object.

Usage

## S3 method for class 'pibblefit'
print(x, summary = FALSE, ...)

Arguments

Value

No direct return, prints out summary

See Also

summary.pibblefit summarizes posterior intervals

Examples

sim <- pibble_sim()
fit <- pibble(sim$Y, sim$X)
print(fit)

Generic Method to Calculate R2 for Fitted Model

Description

Generic Method to Calculate R2 for Fitted Model

Usage

r2(m, ...)

## S3 method for class 'pibblefit'
r2(m, covariates = NULL, ...)

## S3 method for class 'bassetfit'
r2(m, covariates = NULL, components = NULL, ...)

Arguments

Details

Calculates Posterior over Linear Model R2 as:

1-\frac{SS_{res}}{SS_{tot}}

where SS is defined in terms of trace of variances

Method of calculating R2 is multivariate version of the Bayesian R2 proposed by Gelman, Goodrich, Gabry, and Vehtari, 2019

Calculates Posterior over Basset Model R2 as:

1-\frac{SS_{res}}{SS_{tot}}

Method of calculating R2 is multivariate version of the Bayesian R2 proposed by Gelman, Goodrich, Gabry, and Vehtari, 2019

Value

vector

Provide random initialization for pibble model

Description

Randomly initializes based on ALR transform of counts plus random pseudocounts uniformily distributed between 0 and 1.

Usage

random_pibble_init(Y)

Arguments

Details

Notation: N is number of samples and D is number of multinomial categories

Value

(D-1) x N matrix

Examples

Y <- matrix(sample(1:100, 100), 10, 10)
random_pibble_init(Y)

Generic method for fitting model from passed model fit object

Description

Generic method for fitting model from passed model fit object

Usage

refit(m, ...)

Arguments

Value

object of the same class as m

Generic method for ensuring object contains required elements

Description

Intended to be called internally by package

Usage

req(m, r)

Arguments

Value

throws error if required element is not present

require elements to be non-null in orthusfit or throw error

Description

require elements to be non-null in orthusfit or throw error

Usage

## S3 method for class 'orthusfit'
req(m, r)

Arguments

Value

None, throws an error if NULL

require elements to be non-null in pibblefit or throw error

Description

require elements to be non-null in pibblefit or throw error

Usage

## S3 method for class 'pibblefit'
req(m, r)

Arguments

Value

Nothing, throws an error if NULL

Generic method for sampling from prior distribution of object

Description

Generic method for sampling from prior distribution of object

Usage

sample_prior(m, n_samples = 2000L, ...)

Arguments

Value

object of the same class

Sample from the prior distribution of pibblefit object

Description

Note this can be used to sample from prior and then predict can be called to get counts or LambdaX (predict.pibblefit)

Usage

## S3 method for class 'pibblefit'
sample_prior(
  m,
  n_samples = 2000L,
  pars = c("Eta", "Lambda", "Sigma"),
  use_names = TRUE,
  ...
)

Arguments

Details

Could be greatly speed up in the future if needed by sampling directly from cholesky form of inverse wishart (currently implemented as header in this library - see MatDist.h).

Value

A pibblefit object

Examples

# Sample prior of already fitted  pibblefit object
sim <- pibble_sim()
attach(sim)
fit <- pibble(Y, X)
head(sample_prior(fit))

# Sample prior as part of model fitting
m <- pibblefit(N=as.integer(sim$N), D=as.integer(sim$D), Q=as.integer(sim$Q), 
                iter=2000L, upsilon=upsilon, 
                Xi=Xi, Gamma=Gamma, Theta=Theta, X=X, 
                coord_system="alr", alr_base=D)
m <- sample_prior(m)
plot(m) # plot prior distribution (defaults to parameter Lambda) 

Holds information on coordinates system to later be reapplied

Description

store_coord stores coordinate information for pibblefit object and can be reapplied with function reapply_coord. Some coordinate systems are not useful for computation and this makes it simple keep returned object from computations in the same coordinate system as the input. (Likely most useful inside of a package)

Usage

store_coord(m)

reapply_coord(m, l)

Arguments

Value

store_coord list with important information to identify c coordinate system of pibblefit object. reapply_coord pibblefit object in coordinate system previously stored.

Shortcut for summarize variable with quantiles and mean

Description

Shortcut for summarize variable with quantiles and mean

Usage

summarise_posterior(data, var, ...)

Arguments

Details

Notation: pX refers to the X% quantile

Value

data.frame

Examples

d <- data.frame("a"=sample(1:10, 50, TRUE),
                "b"=rnorm(50))

# Summarize posterior for b over grouping of a and also calcuate
# minmum of b (in addition to normal statistics returned)
d <- dplyr::group_by(d, a)
summarise_posterior(d, b, mean.b = mean(b), min=min(b))

Summarise orthusfit object and print posterior quantiles

Description

Default calculates median, mean, 50% and 95% credible interval

Usage

## S3 method for class 'orthusfit'
summary(
  object,
  pars = NULL,
  use_names = TRUE,
  as_factor = FALSE,
  gather_prob = FALSE,
  ...
)

Arguments

Value

A list

Summarise pibblefit object and print posterior quantiles

Description

Default calculates median, mean, 50% and 95% credible interval

Usage

## S3 method for class 'pibblefit'
summary(
  object,
  pars = NULL,
  use_names = TRUE,
  as_factor = FALSE,
  gather_prob = FALSE,
  ...
)

Arguments

Value

A list

Uncollapse output from optimPibbleCollapsed to full pibble Model

Description

See details for model. Should likely be called following optimPibbleCollapsed. Notation: N is number of samples, D is number of multinomial categories, Q is number of covariates, iter is the number of samples of eta (e.g., the parameter n_samples in the function optimPibbleCollapsed)

Usage

uncollapsePibble(
  eta,
  X,
  Theta,
  Gamma,
  Xi,
  upsilon,
  seed,
  ret_mean = FALSE,
  ncores = -1L
)

Arguments

Details

Notation: Let Z_j denote the J-th row of a matrix Z. While the collapsed model is given by:

Y_j \sim \text{Multinomial}(\pi_j)

\pi_j = \Phi^{-1}(\eta_j)

\eta \sim T_{D-1, N}(\upsilon, \Theta X, K, A)

Where A = I_N + X \Gamma X', K=\Xi is a (D-1)x(D-1) covariance matrix, \Gamma is a Q x Q covariance matrix, and \Phi^{-1} is ALRInv_D transform.

The uncollapsed model (Full pibble model) is given by:

Y_j \sim \text{Multinomial}(\pi_j)

\pi_j = \Phi^{-1}(\eta_j)

\eta \sim MN_{D-1 \times N}(\Lambda X, \Sigma, I_N)

\Lambda \sim MN_{D-1 x Q}(\Theta, \Sigma, \Gamma)

\Sigma \sim InvWish(\upsilon, \Xi)

This function provides a means of sampling from the posterior distribution of Lambda and Sigma given posterior samples of Eta from the collapsed model.

Value

List with components

1. Lambda Array of dimension (D-1) x Q x iter (posterior samples)

2. Sigma Array of dimension (D-1) x (D-1) x iter (posterior samples)

3. The number of cores used

4. Timer

References

JD Silverman K Roche, ZC Holmes, LA David, S Mukherjee. Bayesian Multinomial Logistic Normal Models through Marginally Latent Matrix-T Processes. 2019, arXiv e-prints, arXiv:1903.11695

See Also

optimPibbleCollapsed

Examples

sim <- pibble_sim()

# Fit model for eta
fit <- optimPibbleCollapsed(sim$Y, sim$upsilon, sim$Theta%*%sim$X, sim$KInv, 
                             sim$AInv, random_pibble_init(sim$Y))  

# Finally obtain samples from Lambda and Sigma
fit2 <- uncollapsePibble(fit$Samples, sim$X, sim$Theta, 
                                   sim$Gamma, sim$Xi, sim$upsilon, 
                                   seed=2849)

Uncollapse output from optimPibbleCollapsed to full pibble Model when Sigma is known

Description

See details for model. Should likely be called following optimPibbleCollapsed. Notation: N is number of samples, D is number of multinomial categories, Q is number of covariates, iter is the number of samples of eta (e.g., the parameter n_samples in the function optimPibbleCollapsed)

Usage

uncollapsePibble_sigmaKnown(
  eta,
  X,
  Theta,
  Gamma,
  GammaComb,
  Xi,
  sigma,
  upsilon,
  seed,
  ret_mean = FALSE,
  linear = FALSE,
  ncores = -1L
)

Arguments

Details

Notation: Let Z_j denote the J-th row of a matrix Z. While the collapsed model is given by:

Y_j \sim Multinomial(\pi_j)

\pi_j = \Phi^{-1}(\eta_j)

\eta \sim T_{D-1, N}(\upsilon, \Theta X, K, A)

Where A = I_N + X \Gamma X', K = \Xi is a (D-1)x(D-1) covariance matrix, \Gamma is a Q x Q covariance matrix, and \Phi^{-1} is ALRInv_D transform.

The uncollapsed model (Full pibble model) is given by:

Y_j \sim Multinomial(\pi_j)

\pi_j = \Phi^{-1}(\eta_j)

\eta \sim MN_{D-1 \times N}(\Lambda X, \Sigma, I_N)

\Lambda \sim MN_{D-1 \times Q}(\Theta, \Sigma, \Gamma)

\Sigma \sim InvWish(\upsilon, \Xi)

This function provides a means of sampling from the posterior distribution of Lambda and Sigma given posterior samples of Eta from the collapsed model.

Value

List with components

1. Lambda Array of dimension (D-1) x Q x iter (posterior samples)

2. Sigma Array of dimension (D-1) x (D-1) x iter (posterior samples)

3. The number of cores used

4. Timer

References

JD Silverman K Roche, ZC Holmes, LA David, S Mukherjee. Bayesian Multinomial Logistic Normal Models through Marginally Latent Matrix-T Processes. 2019, arXiv e-prints, arXiv:1903.11695

See Also

optimPibbleCollapsed

Generic method for verifying new objects

Description

Intended to be called internally by package or object creator

Usage

verify(m, ...)

Arguments

Value

throws error if verify test fails

Simple verification of passed bassetfit object

Description

Simple verification of passed bassetfit object

Usage

## S3 method for class 'bassetfit'
verify(m, ...)

Arguments

Value

throws error if any verification tests fail

Simple verification of passed orthusfit object

Description

Simple verification of passed orthusfit object

Usage

## S3 method for class 'orthusfit'
verify(m, ...)

Arguments

Value

throws error if any verification tests fail

Simple verification of passed pibblefit object

Description

Simple verification of passed pibblefit object

Usage

## S3 method for class 'pibblefit'
verify(m, ...)

Arguments

Value

throws error if any verification tests fail