mixlink (version 0.1.5)

Metropolis-within-Gibbs Sampler: Metropolis-within-Gibbs Sampler

Description

Metropolis-within-Gibbs sampler for Binomial and Poisson Mixture Link models.

Usage

gibbs.mixlink(y, X, R, burn = 0, thin = 1, invlink = NULL,
  report.period = R + 1, save.psi = FALSE, Beta.init = NULL, Pi.init,
  kappa.init = NULL, psi.init = NULL, proposal.VBeta = NULL,
  proposal.VPi = NULL, proposal.Vkappa = NULL, proposal.Vpsi = NULL,
  hyper = NULL, trials = NULL, offset = rep(0, length(y)),
  family = NULL, fixed.kappa = FALSE)

Arguments

y

Observations.

X

Design matrix for regression.

R

Number of MCMC draws to take.

burn

Number of initial MCMC draws to discard.

thin

After burn-in period, save one of every thin draws.

invlink

The inverse link function for the mean. Default is NULL, which indicates to use plogis for Binomial and exp for Poisson.

report.period

Report progress every report.period draws.

save.psi

Save draws for \(\bm{\psi}_1, \ldots, \bm{\psi}_n\). Default is FALSE.

Beta.init

Starting value for \(\bm{\beta}\).

Pi.init

Starting value for \(\bm{\pi}\). The length of Pi.init determines \(J\) used in MCMC.

kappa.init

Starting value for \(\kappa\).

psi.init

Starting value for \(\bm{\psi}_1, \ldots, \bm{\psi}_n\)

proposal.VBeta

Covariance matrix for \(d\) dimensional multivariate normal proposal. If left at the default (NULL), we use diag(0.01^2, d).

proposal.VPi

Covariance matrix for \(J-1\) dimensional multivariate normal proposal. If left at the default (NULL), we use diag(0.01^2, J-1).

proposal.Vkappa

Covariance matrix for univariate normal proposal. If left at the default (NULL), we use 0.01^2.

proposal.Vpsi

Covariance matrix for univariate normal proposal. If left at the default (NULL), we use diag(0.01^2, J).

hyper

A list with hyperparameters corresponding to the prior from the Details section. var.Beta is \(V_{\beta}\) with default 1000. alpha.Pi is \(\bm{\alpha}_\pi\) with default rep(1,J). a.kappa and b.kappa correspond to \((a_\kappa, b_\kappa)\) which have default values a.kappa = 1 and b.kappa = 1/10.

trials

Number of success/failure trials.

offset

Offset term to add to regression function, as commonly used in count models.

family

Can be either binomial or poisson.

fixed.kappa

Keep \(\kappa\) fixed at kappa.init (Default FALSE).

Value

A list with the MCMC results:

Beta.hist

\(R \times d\) matrix of saved \(\bm{\beta}\) draws.

Pi.hist

\(R \times J\) matrix of saved \(\bm{\pi}\) draws.

kappa.hist

\(R \times 1\) vector of \(\kappa\) draws.

psi.hist

A list of \(J\) \(R \times n\) matrices. The \(j\)th matrix contains the MCMC draws for \(\bm{\psi}_{ij}\).

accept.Beta

Percentages of MCMC proposals for \(\bm{\beta}\) were accepted.

accept.Pi

Percentages of MCMC proposals for \(\bm{\pi}\) were accepted.

accept.kappa

Percentages of MCMC proposals for \(\kappa\) were accepted.

accept.psi

Percentages of MCMC proposals for \(\bm{\psi}\) were accepted.

elapsed.sec

Elapsed time for sampling, in seconds.

R

Number of total draws.

R.keep

Number of draws kept, after thinning and burn-in.

X.names

Names of columns of design matrix.

family

Name of family / outcome.

Can be accessed with the functions print, summary, and DIC.

Details

Priors for Bayesian Mixture Link model are

• \( \bm{\beta} \sim \textrm{N}(\bm{0}, V_{\beta} \bm{I}) \),

• \( \bm{\pi} \sim \textrm{Dirichlet}_J(\bm{\alpha}_\pi) \),

• \( \kappa \sim \textrm{Gamma}(a_\kappa, b_\kappa) \), parameterized with \(\textrm{E}(\kappa) = a_\kappa / b_\kappa\).

References

Andrew M. Raim, Nagaraj K. Neerchal, and Jorge G. Morel. An Extension of Generalized Linear Models to Finite Mixture Outcomes. arXiv preprint: 1612.03302

Examples

# NOT RUN {
# ----- Generate data -----
n <- 200
m <- rep(20, n)
x <- rnorm(n, 0, 1)
X <- model.matrix(~ x)
Beta.true <- c(-1, 1)
mean.true <- plogis(X %*% Beta.true)
kappa.true <- 1
Pi.true <- c(1,3)/4
d <- ncol(X)
J <- length(Pi.true)
y <- r.mixlink.binom(n, mean.true, Pi.true, kappa.true, m)

# ----- Run Metropolis-within-Gibbs sampler -----
hyper <- list(VBeta = diag(1000, d), alpha.Pi = rep(1, J),
	kappa.a = 1, kappa.b = 1/2)
gibbs.out <- gibbs.mixlink(y, X, R = 10, burn = 5, thin = 1,
	invlink = plogis, report.period = 100, Pi.init = c(1,9)/10,
	proposal.VBeta = solve(t(X) %*% X), proposal.VPi = diag(0.25^2, J-1),
	proposal.Vkappa = 0.5^2, proposal.Vpsi = diag(0.5^2, J),
	hyper = hyper, family = "binomial", trials = m)

print(gibbs.out)
DIC(gibbs.out)

# }