Metropolis-within-Gibbs-Sampler: Metropolis-within-Gibbs Sampler
In mixlink: Mixture Link Regression

Description Usage Arguments Details Value References Examples

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

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)

`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{ψ}_1, …, \bm{ψ}_n. Default is `FALSE`.
`Beta.init`	Starting value for \bm{β}.
`Pi.init`	Starting value for \bm{π}. The length of `Pi.init` determines J used in MCMC.
`kappa.init`	Starting value for κ.
`psi.init`	Starting value for \bm{ψ}_1, …, \bm{ψ}_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_{β} with default 1000. `alpha.Pi` is \bm{α}_π with default `rep(1,J)`. `a.kappa` and `b.kappa` correspond to (a_κ, b_κ) 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 κ fixed at `kappa.init` (Default FALSE).

Priors for Bayesian Mixture Link model are

\bm{β} \sim \textrm{N}(\bm{0}, V_{β} \bm{I}) ,
\bm{π} \sim \textrm{Dirichlet}_J(\bm{α}_π) ,
κ \sim \textrm{Gamma}(a_κ, b_κ) , parameterized with \textrm{E}(κ) = a_κ / b_κ.

A list with the MCMC results:

`Beta.hist`	R \times d matrix of saved \bm{β} draws.
`Pi.hist`	R \times J matrix of saved \bm{π} draws.
`kappa.hist`	R \times 1 vector of κ draws.
`psi.hist`	A list of J R \times n matrices. The jth matrix contains the MCMC draws for \bm{ψ}_{ij}.
`accept.Beta`	Percentages of MCMC proposals for \bm{β} were accepted.
`accept.Pi`	Percentages of MCMC proposals for \bm{π} were accepted.
`accept.kappa`	Percentages of MCMC proposals for κ were accepted.
`accept.psi`	Percentages of MCMC proposals for \bm{ψ} 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.

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

# ----- 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)