This function runs MCMC for a generalized linear mixed model with possibly several response variables and possibly normal mixtures in the distributions of random effects.
GLMM_MCMC(y, dist = "gaussian", id, x, z, random.intercept,
prior.alpha, init.alpha, init2.alpha,
scale.b, prior.b, init.b, init2.b,
prior.eps, init.eps, init2.eps,
nMCMC = c(burn = 10, keep = 10, thin = 1, info = 10),
tuneMCMC = list(alpha = 1, b = 1),
store = c(b = FALSE), PED = TRUE, keep.chains = TRUE,
dens.zero = 1e-300, parallel = FALSE, cltype, silent = FALSE)# S3 method for GLMM_MCMC
print(x, ...)
# S3 method for GLMM_MCMClist
print(x, ...)
An object of class GLMM_MCMClist (if PED argument is
TRUE) or GLMM_MCMC (if PED argument is
FALSE).
vector, matrix or data frame with responses. If y is vector then
there is only one response in the model. If y is matrix or data frame then
each column gives values of one response. Missing values are allowed.
If there are several responses specified then continuous responses must be put in the first columns and discrete responses in the subsequent columns.
character (vector) which determines distribution (and a link function) for each response variable. Possible values are: “gaussian” for gaussian (normal) distribution (with identity link), “binomial(logit)” for binomial (0/1) distribution with a logit link. “poisson(log)” for Poisson distribution with a log link. Single value is recycled if necessary.
vector which determines longitudinally or otherwise dependent observations. If not given then it is assumed that there are no clusters and all observations of one response are independent.
matrix or a list of matrices with covariates (intercept not included) for fixed effects.
If there is more than one response, this must always be a list. Note that intercept in included
in all models. Use a character value “empty” as a component of the list x
if there are no covariates for a particular response.
matrix or a list of matrices with covariates (intercept not included) for random effects.
If there is more than one response, this must always be a list. Note that random intercept
is specified using the argument random.intercept.
REMARK: For a particular response, matrices x and z
may not have the same columns. That is, matrix x includes
covariates which are not involved among random effects and matrix
z includes covariates which are involved among random effects
(and implicitely among fixed effects as well).
logical (vector) which determines for which responses random intercept should be included.
a list which specifies prior distribution for fixed
effects (not the means of random effects). The prior distribution is
normal and the user can specify the mean and variances.
The list prior.alpha can have the components listed below.
a vector with prior means, defaults to zeros.
a vector with prior variances, defaults to 10000 for all components.
a numeric vector with initial values of fixed effects
(not the means of random effects) for the first chain. A sensible value is determined using the
maximum-likelihood fits (using lmer functions)
and does not have to be given by the user.
a numeric vector with initial values of fixed effects for the second chain.
a list specifying how to scale the random effects during
the MCMC. A sensible value is determined using the
maximum-likelihood fits (using lmer functions)
and does not have to be given by the user.
If the user wishes to influence the shift and scale constants, these
are given as components of the list scale.b. The components
are named:
see vignette PBCseq.pdf for details
see vignette PBCseq.pdf for details
a list which specifies prior distribution for (shifted and scaled) random effects. The prior is in principle a normal mixture (being a simple normal distribution if we restrict the number of mixture components to be equal to one).
The list prior.b can have the components listed below. Their
meaning is analogous to the components of the same name of the
argument prior of function NMixMCMC (see
therein for details).
a character string which specifies the assumed prior distribution for random effects. It can be either “normal” (multivaruate normal - default) or “MVT” (multivariate Student t distribution).
a character string which specifies the type of the prior for \(K\) (the number of mixture components).
a character string which specifies the type of the prior for mixture means and mixture variances.
maximal number of mixture components.
see vignette PBCseq.pdf for details
see vignette PBCseq.pdf for details
see vignette PBCseq.pdf for details
see vignette PBCseq.pdf for details
see vignette PBCseq.pdf for details
see vignette PBCseq.pdf for details
see vignette PBCseq.pdf for details
see vignette PBCseq.pdf for details
shape parameter of the prior distribution for the degrees of freedom if the random effects are assumed to follow the MVT distribution
rate parameter of the prior distribution for the degrees of freedom if the random effects are assumed to follow the MVT distribution
a list with initial values of the first chain for parameters related to the distribution of random effects and random effects themselves. Sensible initial values are determined by the function itself and do not have to be given by the user.
a list with initial values of the second chain for parameters related to the distribution of random effects and random effects themselves.
a list specifying prior distributions for
error terms for continuous responses. The list prior.eps can
have the components listed below. For all components, a sensible
value leading to weakly informative prior distribution can be
determined by the function.
see vignette PBCseq.pdf for details
see vignette PBCseq.pdf for details
see vignette PBCseq.pdf for details
a list with initial values of the first chain for parameters related to the
distribution of error terms of continuous responses. The list
init.eps can have the components listed below. For all
components, a sensible value can be determined by the function.
a numeric vector with the initial values for residual standard deviations for each continuous response.
a numeric vector with the initial values for the inverted components of the hyperparameter gamma for each continuous response.
a list with initial values of the second chain for parameters related to the distribution of error terms of continuous responses.
numeric vector of length 4 giving parameters of the MCMC simulation. Its components may be named (ordering is then unimportant) as:
length of the burn-in (after discarding the thinned values), can be equal to zero as well.
length of the kept chains (after discarding the thinned values), must be positive.
thinning interval, must be positive.
interval in which the progress information is printed on the screen.
In total \((M_{burn} + M_{keep}) \cdot M_{thin}\) MCMC scans are performed.
a list with tuning scale parameters for proposal distribution of fixed and random effects. It is used only when there are some discrete response profiles. The components of the list have the following meaning:
scale parameters by which we multiply the proposal covariance matrix when updating the fixed effects pertaining to the discrete response profiles. There is one scale parameter for each DISCRETE profile. A single value is recycled if necessary.
a scale parameter by which we multiply the proposal covariance matrix when updating the random effects. It is used only when there are some discrete response profiles in the model.
logical vector indicating whether the chains of parameters should be stored. Its components may be named (ordering is then unimportant) as:
if TRUE then the sampled values of random effects
are stored. Defaults to FALSE.
a logical value which indicates whether the penalized expected deviance (see Plummer, 2008 for more details) is to be computed (which requires two parallel chains).
logical. If FALSE, only summary statistics
are returned in the resulting object. This might be useful in the
model searching step to save some memory.
a small value used instead of zero when computing deviance related quantities.
a logical value which indicates whether parallel
computation (based on a package parallel)
should be used when running two chains for the purpose of PED
calculation.
optional argument applicable if parallel is
TRUE. If cltype is given, it is passed as the
type argument into the call to makeCluster.
a logical value indicating whether the information on the MCMC progress is to be supressed.
additional arguments passed to the default print method.
Object of class GLMM_MCMC can have the following
components (some of them may be missing according to the context
of the model):
index of the last iteration performed.
used value of the argument nMCMC.
a character vector of length R corresponding to the
dist argument.
a two component vector giving the number of continuous responses and the number of discrete responses.
a numeric vector of length R giving the number of non-intercept alpha parameters for each response.
a numeric vector of length R giving the number of non-intercept random effects for each response.
a logical vector of length R which indicates inclusion of fixed intercept for each response.
a logical vector of length R which indicates inclusion of random intercept for each response.
length of the vector of fixed effects.
dimension of the distribution of random effects.
a list containing the used value of the
argument prior.alpha.
a list containing the used value of the
argument prior.b.
a list containing the used value of the
argument prior.eps.
a numeric vector with the used value of the
argument init.alpha.
a list containing the used value of the
argument init.b.
a list containing the used value of the
argument init.eps.
a numeric vector with the first stored (after burn-in) value of fixed effects \(\alpha\).
a numeric vector with the last sampled value
of fixed effects \(\alpha\). It can be used as argument
init.alpha to restart MCMC.
a list with the first stored (after burn-in) values of parameters
related to the distribution of random effects. It has components
named b, K, w, mu, Sigma, Li, Q,
gammaInv, r.
a list with the last sampled values of parameters
related to the distribution of random effects. It has components
named b, K, w, mu, Sigma, Li, Q,
gammaInv, r. It can be used as argument
init.b to restart MCMC.
a list with the first stored (after burn-in) values of parameters
related to the distribution of residuals of continuous responses. It
has components named sigma, gammaInv.
a list with the last sampled values of parameters
related to the distribution of residuals of continuous responses. It
has components named sigma, gammaInv. It can be used as argument
init.eps to restart MCMC.
acceptance proportion from the Metropolis-Hastings algorithm for fixed effects (separately for each response type). Note that the acceptance proportion is equal to one for continuous responses since the Gibbs algorithm is used there.
acceptance proportion from the Metropolis-Hastings algorithm for random effects (separately for each cluster). Note that the acceptance proportion is equal to one for models with continuous responses only since the Gibbs algorithm is used there.
a list containing the used value of the argument
scale.b.
a data.frame with posterior summary
statistics for the deviance (approximated using the Laplacian
approximation) and conditional (given random effects) devience.
a data.frame with posterior summary statistics for fixed effects.
a matrix with posterior summary statistics for means of random effects.
a matrix with posterior summary statistics for standard deviations of random effects and correlations of each pair of random effects.
a matrix with posterior summary statistics for standard deviations of the error terms in the (mixed) models of continuous responses.
a matrix which is present in the output object
if the number of mixture components in the distribution of random
effects is fixed and equal to \(K\). In that case,
poster.comp.prob_u is a matrix with \(K\) columns and \(I\)
rows (\(I\) is the number of subjects defining the longitudinal
profiles or correlated observations) with estimated posterior component probabilities
-- posterior means of the components of the underlying 0/1
allocation vector.
WARNING: By default, the labels of components are based on artificial identifiability constraints based on ordering of the mixture means in the first margin. Very often, such identifiability constraint is not satisfactory!
a matrix which is present in the output object
if the number of mixture components in the distribution of random
effects is fixed and equal to \(K\). In that case,
poster.comp.prob_b is a matrix with \(K\) columns and \(I\)
rows (\(I\) is the number of subjects defining the longitudinal
profiles or correlated observations)
with estimated posterior component probabilities
-- posterior mean over model parameters including random effects.
WARNING: By default, the labels of components are based on artificial identifiability constraints based on ordering of the mixture means in the first margin. Very often, such identifiability constraint is not satisfactory!
frequency table for the MCMC sample of the number of mixture components in the distribution of the random effects.
posterior probabilities for the numbers of mixture components in the distribution of random effects.
a list with data.frames, one
data.frame per response profile. Each data.frame
with columns labeled id, observed,
fitted, stres,
eta.fixed and eta.random holding
identifier for clusters of grouped observations,
observed values and
posterior means for fitted values (response expectation given fixed and random effects),
standardized residuals (derived from fitted values),
fixed effect part of the linear predictor and the random effect part of
the linear predictor. In each column, there are first all values for
the first response, then all values for the second response etc.
a data.frame with columns labeled
b1, ..., bq, Logpb, Cond.Deviance, Deviance with
posterior means of random effects for each cluster, posterior
means of \(\log\bigl\{p(\boldsymbol{b})\bigr\}\),
conditional deviances, i.e., minus twice the conditional (given
random effects) log-likelihood for each cluster
and GLMM deviances, i.e., minus twice the marginal (random effects
integrated out) log-likelihoods for each cluster. The value of the
marginal (random effects integrated out) log-likelihood at each MCMC
iteration is obtained using the Laplacian approximation.
a numeric vector with posterior means of mixture weights after re-labeling. It is computed only if \(K_b\) is fixed and even then I am not convinced that these are useful posterior summary statistics (see label switching problem mentioned above). In any case, they should be used with care.
a matrix with posterior means of mixture means after re-labeling. It is computed only if \(K_b\) is fixed and even then I am not convinced that these are useful posterior summary statistics (see label switching problem mentioned above). In any case, they should be used with care.
a list with posterior means of mixture inverse variances after re-labeling. It is computed only if \(K_b\) is fixed and even then I am not convinced that these are useful posterior summary statistics (see label switching problem mentioned above). In any case, they should be used with care.
a list with posterior means of mixture variances after re-labeling. It is computed only if \(K_b\) is fixed and even then I am not convinced that these are useful posterior summary statistics (see label switching problem mentioned above). In any case, they should be used with care.
a list with posterior means of Cholesky decompositions of mixture inverse variances after re-labeling. It is computed only if \(K_b\) is fixed and even then I am not convinced that these are useful posterior summary statistics (see label switching problem mentioned above). In any case, they should be used with care.
numeric vector with a chain for the GLMM deviances, i.e., twice the marginal (random effects integrated out) log-likelihoods of the GLMM. The marginal log-likelihood is obtained using the Laplacian approximation at each iteration of MCMC.
numeric vector with a chain for the conditional deviances, i.e., twice the conditional (given random effects) log-likelihoods.
numeric vector with a chain for \(K_b\) (number of mixture components in the distribution of random effects).
numeric vector or matrix with a chain for \(w_b\) (mixture weights for the distribution of random effects). It is a matrix with \(K_b\) columns when \(K_b\) is fixed. Otherwise, it is a vector with weights put sequentially after each other.
numeric vector or matrix with a chain for \(\mu_b\) (mixture means for the distribution of random effects). It is a matrix with \(dimb\cdot K_b\) columns when \(K_b\) is fixed. Otherwise, it is a vector with means put sequentially after each other.
numeric vector or matrix with a chain for lower triangles of \(\boldsymbol{Q}_b\) (mixture inverse variances for the distribution of random effects). It is a matrix with \(\frac{dimb(dimb+1)}{2}\cdot K_b\) columns when \(K_b\) is fixed. Otherwise, it is a vector with lower triangles of \(\boldsymbol{Q}_b\) matrices put sequentially after each other.
numeric vector or matrix with a chain for lower triangles of \(\Sigma_b\) (mixture variances for the distribution of random effects). It is a matrix with \(\frac{dimb(dimb+1)}{2}\cdot K_b\) columns when \(K_b\) is fixed. Otherwise, it is a vector with lower triangles of \(\Sigma_b\) matrices put sequentially after each other.
numeric vector or matrix with a chain for lower triangles of Cholesky decompositions of \(\boldsymbol{Q}_b\) matrices. It is a matrix with \(\frac{dimb(dimb+1)}{2}\cdot K_b\) columns when \(K_b\) is fixed. Otherwise, it is a vector with lower triangles put sequentially after each other.
matrix with \(dimb\) columns with a chain for inverses of the hyperparameter \(\boldsymbol{\gamma}_b\).
numeric vector or matrix with order indeces of mixture components in the distribution of random effects related to artificial identifiability constraint defined by ordering of the first component of the mixture means.
It is a matrix with \(K_b\) columns when \(K_b\) is fixed. Otherwise it is a vector with orders put sequentially after each other.
numeric vector or matrix with rank indeces of mixture components in the distribution of random effects related to artificial identifiability constraint defined by ordering of the first component of the mixture means.
It is a matrix with \(K_b\) columns when \(K_b\) is fixed. Otherwise it is a vector with ranks put sequentially after each other.
data.frame with columns labeled
b.Mean.*, b.SD.*, b.Corr.*.*
containing the chains for the means, standard deviations and correlations of the
distribution of the random effects based on a normal mixture at each
iteration.
a matrix with the MCMC chains for random effects. It is
included only if store[b] is TRUE.
numeric vector or matrix with the MCMC chain(s) for fixed effects.
numeric vector or matrix with the MCMC chain(s) for standard deviations of the error terms in the (mixed) models for continuous responses.
matrix with \(dimb\) columns with MCMC chain(s) for inverses of the hyperparameter \(\boldsymbol{\gamma}_b\).
a list which specifies the algorithm used to re-label
the MCMC output to compute order_b, rank_b, poster.comp.prob_u,
poster.comp.prob_b, poster.mean.w_b,
poster.mean.mu_b, poster.mean.Q_b,
poster.mean.Sigma_b, poster.mean.Li_b.
a list with components useful to call underlying C++ functions (not interesting for ordinary users).
Object of class NMixMCMClist is the list having two components
of class NMixMCMC representing two parallel chains and
additionally the following components:
values of penalized expected deviance and related
quantities. It is a vector with five components: D.expect \(=\)
estimated expected deviance, where the estimate is based on two
parallel chains; popt \(=\) estimated penalty, where the
estimate is based on simple MCMC average based on two parallel
chains; PED \(=\) estimated penalized expected deviance
\(=\) D.expect \(+\) popt; wpopt \(=\)
estimated penalty, where the estimate is based on weighted MCMC average
(through importance sampling) based on two parallel chains;
wPED \(=\) estimated penalized expected deviance \(=\)
D.expect \(+\) wpopt.
posterior mean of the deviance for each subject.
contributions to the unweighted penalty from each subject.
contributions to the weighted penalty from each subject.
for each subject, number of iterations (in both chains), where the
deviance was in fact equal to infinity (when the corresponding
density was lower than dens.zero) and was not taken into account when
computing D.expect.
for each subject, number of iterations, where the
penalty was in fact equal to infinity and was not taken into account
when computing popt.
for each subject, number of iterations, where the
importance sampling weight was in fact equal to infinity and was not taken into account
when computing wpopt.
for each subject, sum of importance sampling weights.
sampled value of the observed data deviance from chain 1
sampled values of the obserbed data deviance from chain 2
sampled values of the deviance of data replicated according to the chain 1 evaluated under the parameters from chain 1
sampled values of the deviance of data replicated according to the chain 1 evaluated under the parameters from chain 2
sampled values of the deviance of data replicated according to the chain 2 evaluated under the parameters from chain 1
sampled values of the deviance of data replicated according to the chain 2 evaluated under the parameters from chain 2
Arnošt Komárek arnost.komarek@mff.cuni.cz
See accompanying papers (Komárek et al., 2010, Komárek and Komárková, 2013).
Komárek, A. and Komárková, L. (2013). Clustering for multivariate continuous and discrete longitudinal data. The Annals of Applied Statistics, 7(1), 177--200.
Komárek, A. and Komárková, L. (2014). Capabilities of R package mixAK for clustering based on multivariate continuous and discrete longitudinal data. Journal of Statistical Software, 59(12), 1--38. tools:::Rd_expr_doi("10.18637/jss.v059.i12").
Komárek, A., Hansen, B. E., Kuiper, E. M. M., van Buuren, H. R., and Lesaffre, E. (2010). Discriminant analysis using a multivariate linear mixed model with a normal mixture in the random effects distribution. Statistics in Medicine, 29(30), 3267--3283.
Plummer, M. (2008). Penalized loss functions for Bayesian model comparison. Biostatistics, 9(3), 523--539.
NMixMCMC.
## See also additional material available in
## YOUR_R_DIR/library/mixAK/doc/
## or YOUR_R_DIR/site-library/mixAK/doc/
## - files http://www.karlin.mff.cuni.cz/~komarek/software/mixAK/PBCseq.pdf,
## PBCseq.R
## ==============================================
Run the code above in your browser using DataLab