Description Usage Arguments Details Value References Examples
Metropolis-within-Gibbs sampler for Binomial and Poisson Mixture Link models.
1 2 3 4 5 6 | 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 |
invlink |
The inverse link function for the mean. Default is |
report.period |
Report progress every |
save.psi |
Save draws for \bm{ψ}_1, …, \bm{ψ}_n. Default is
|
Beta.init |
Starting value for \bm{β}. |
Pi.init |
Starting value for \bm{π}. The length of |
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 |
proposal.VPi |
Covariance matrix for J-1 dimensional multivariate normal
proposal. If left at the default (NULL), we use |
proposal.Vkappa |
Covariance matrix for univariate normal proposal. If left at
the default (NULL), we use |
proposal.Vpsi |
Covariance matrix for univariate normal proposal. If left at
the default (NULL), we use |
hyper |
A list with hyperparameters corresponding to the prior from the
Details section. |
trials |
Number of success/failure trials. |
offset |
Offset term to add to regression function, as commonly used in count models. |
family |
Can be either |
fixed.kappa |
Keep κ fixed at |
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
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 | # ----- 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)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.