Metropolis-Hastings-Sampler: Metropolis-Hastings Sampler
In mixlink: Mixture Link Regression

Description Usage Arguments Details Value References Examples

Random-Walk Metropolis Hastings sampler for Binomial and Poisson Mixture Link models.

rwmetrop.mixlink.binomial(y, X, m, R, burn = 0, thin = 1,
  invlink = plogis, Beta.init = NULL, Pi.init, kappa.init = NULL,
  hyper = NULL, report.period = R + 1, use.laplace.approx = TRUE,
  proposal = NULL, param.grp = NULL, fixed.kappa = FALSE)

rwmetrop.mixlink.poisson(y, X, offset = rep(0, length(y)), R, burn = 0,
  thin = 1, invlink = exp, Beta.init = NULL, Pi.init, kappa.init = NULL,
  hyper = NULL, report.period = R + 1, use.laplace.approx = TRUE,
  proposal = NULL, param.grp = NULL, fixed.kappa = FALSE)

`y`	Observations.
`X`	Design matrix for regression.
`m`	Number of success/failure trials.
`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 `plogis` for Binomial and `exp` for Poisson.
`Beta.init`	Starting value for \bm{β}.
`Pi.init`	Starting value for \bm{π}.
`kappa.init`	Starting value for κ.
`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`.
`report.period`	Report progress every `report.period` draws.
`use.laplace.approx`	Maximize a Laplace approximation to the posterior, to find a starting value for MCMC.
`proposal`	A list with two elements. `var` is the covariance matrix for a d + (J-1) + 1 dimensional multivariate normal proposal, `scale` is a scalar multiplied with `var`. Defaults are `var = diag(d+(J-1)+1)` and `scale = 0.02`.
`param.grp`	A vector of integers of length d + (J-1) + 1, where `d = ncol(X)`, which indicates the grouping of parameters in MCMC. Parameters with common integers are sampled together. The first d correspond to \bm{β}, the next J-1 correspond to \bm{π}, and the last one corresponds to κ. At the default value (`param.grp = NULL`), \bm{β}, \bm{π}, and κ are sampled one at a time, each in their entirety.
`fixed.kappa`	Keep κ fixed at `kappa.init` (Default FALSE).
`offset`	Constant offset term to add to x^T β.

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:

`par.hist`	R \times [d + (J-1) + 1] matrix of saved MCMC draws before transformations are applied. Most users will not need this.
`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
`accept`	Percentages that MCMC proposals were accepted. Corresponds to `param.grp`
`elapsed.sec`	Elapsed time for sampling, in seconds.
`laplace.out`	Output of Laplace approximation.
`R.keep`	Number of draws kept, after thinning and burn-in.
`X.names`	Names of columns of design matrix.

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

## Not run: 
library(Matrix)

# ----- Generate data -----
n <- 200
x <- runif(n, 1, 3)
X <- model.matrix(~ x)
Beta.true <- c(0, 1)
mean.true <- exp(X %*% Beta.true)
kappa.true <- 0.95
Pi.true <- c(1,3)/4
d <- ncol(X)
J <- length(Pi.true)
y <- r.mixlink.pois(n, mean.true, Pi.true, kappa.true)

# ----- Run Metropolis-within-Gibbs sampler -----
hyper <- list(VBeta = diag(1000, d), alpha.Pi = rep(1, J),
	kappa.a = 1, kappa.b = 1/2)
proposal <- list(
var = bdiag(solve(t(X) %*% X), diag(J-1), 1),
scale = 0.5)
metrop.out <- rwmetrop.mixlink.poisson(y, X, R = 20000, burn = 1000,
	thin = 10, Pi.init = c(1,9)/10, hyper = hyper,
	report.period = 1000, use.laplace.approx = TRUE, proposal = proposal)

print(metrop.out)
DIC(metrop.out)

## End(Not run)