hUM.post: Posterior sampling from a hierarchical...

In MADPop: MHC Allele-Based Differencing Between Populations

Posterior sampling from a hierarchical Unconstrained-Multinomial model

Description

MCMC sampling from a Dirichlet-Multinomial model using stan.

Usage

hUM.post(nsamples, X, popId, rhoId, full.stan.out = FALSE, ...)

Arguments

nsamples	Number of posterior samples
X	4-column or 5-column matrix of observations in the correct format. See UM.suff.
popId	Optional vector of population identifiers. See UM.suff.
rhoId	Populations for which posterior samples of the genotype probability vector rho are desired. Defaults to all populations. Set rhoId = NULL not to output these for any populations.
full.stan.out	Logical. Whether or not to return the full stan output. For monitoring convergence of the MCMC sampling.
...	Further arguments to be passed to the sampling function in rstan.

Details

The hierarchical Dirichlet-Multinomial model is given by

Y_k \mid \rho_k \sim_{\textrm{ind}} \textrm{Multinomial}(\rho_k, N_k),

\rho_k \sim_{\textrm{iid}} \textrm{Dirichlet}(\alpha).

where \alpha_0 = \sum_{i=1}^C \alpha_i and \bar \alpha = \alpha/\alpha_0. MCMC sampling is achieved with the rstan package, which is listed as a dependency for MADPop so as to expose rstan's sophisticated tuning mechanism and convergence diagnostics.

Value

A list with elements

• A: The unique allele names.

• G: The 4-column matrix Package libcurl was not found in the pkg-config search path.of unique genotype combinations.

• rho: A matrix with ncol(rho) == nrow(G), where each row is a draw from the posterior distribution of inheritance probabilities.

• sfit: If full.stan.out = TRUE, the fitted stan object.

Examples

# fit hierarchical model to fish215 data

# only output posterior samples for lake Simcoe
rhoId <- "Simcoe"
nsamples <- 500
hUM.fit <- hUM.post(nsamples = nsamples, X = fish215,
                    rhoId = rhoId,
                    chains = 1) # number of MCMC chains

# plot first 20 posterior probabilities in lake Simcoe
rho.post <- hUM.fit$rho[,1,]
boxplot(rho.post[,1:20], las = 2,
        xlab = "Genotype", ylab = "Posterior Probability",
        pch = ".", col = "grey")

MADPop documentation built on Oct. 13, 2023, 5:09 p.m.