GLMM_MCMC: MCMC estimation of a (multivariate) generalized linear mixed model with a normal mixture in the distribution of random effects

Description

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.

Usage

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)
## S3 method for class 'GLMM_MCMC':
print(x, \dots)
## S3 method for class 'GLMM_MCMClist':
print(x, \dots)

Arguments

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.

dist

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

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

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.

random.intercept

logical (vector) which determines for which responses random intercept should be included.

prior.alpha

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.b can have the components listed below

init.alpha

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)

init2.alpha

a numeric vector with initial values of fixed effects for the second chain.

scale.b

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 u

prior.b

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

init.b

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.

init2.b

a list with initial values of the second chain for parameters related to the distribution of random effects and random effects themselves.

prior.eps

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

init.eps

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 d

init2.eps

a list with initial values of the second chain for parameters related to the distribution of error terms of continuous responses.

nMCMC

numeric vector of length 4 giving parameters of the MCMC simulation. Its components may be named (ordering is then unimportant) as: [object Object],[object Object],[object Object],[object Object] In total $(M_{burn} + M_{keep}) \cdot M_{thin}$

tuneMCMC

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: [object Object],[object Objec

store

logical vector indicating whether the chains of parameters should be stored. Its components may be named (ordering is then unimportant) as: [object Object]

PED

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

keep.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.

dens.zero

a small value used instead of zero when computing deviance related quantities.

parallel

a logical value which indicates whether parallel computation (based on packages snow and snowfall-a-package) should be used when running two chains for the purpos

...

additional arguments passed to the default print method.

Value

An object of class GLMM_MCMC. It can have the following components (some of them may be missing according to the context of the model):
iterindex of the last iteration performed.
nMCMCused value of the argument nMCMC.
dista character vector of length R corresponding to the dist argument.
Ra two component vector giving the number of continuous responses and the number of discrete responses.
pa numeric vector of length R giving the number of non-intercept alpha parameters for each response.
qa numeric vector of length R giving the number of non-intercept random effects for each response.
fixed.intercepta logical vector of length R which indicates inclusion of fixed intercept for each response.
random.intercepta logical vector of length R which indicates inclusion of random intercept for each response.
lalphalength of the vector of fixed effects.
dimbdimension of the distribution of random effects.
prior.alphaa list containing the used value of the argument prior.alpha.
prior.ba list containing the used value of the argument prior.b.
prior.epsa list containing the used value of the argument prior.eps.
init.alphaa numeric vector with the used value of the argument init.alpha.
init.ba list containing the used value of the argument init.b.
init.epsa list containing the used value of the argument init.eps.
state.first.alphaa numeric vector with the first stored (after burn-in) value of fixed effects $\alpha$.
state.last.alphaa numeric vector with the last sampled value of fixed effects $\alpha$. It can be used as argument init.alpha to restart MCMC.
state.first.ba 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.
state.last.ba 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.
state.first.epsa 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.
state.last.epsa 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.
prop.accept.alphaacceptance 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.
prop.accept.bacceptance 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.
scale.ba list containing the used value of the argument scale.b.
summ.Deviancea data.frame with posterior summary statistics for the deviance (approximated using the Laplacian approximation) and conditional (given random effects) devience.
summ.alphaa data.frame with posterior summary statistics for fixed effects.
summ.b.Meana matrix with posterior summary statistics for means of random effects.
summ.b.SDCorra matrix with posterior summary statistics for standard deviations of random effects and correlations of each pair of random effects.
summ.sigma_epsa matrix with posterior summary statistics for standard deviations of the error terms in the (mixed) models of continuous responses.
poster.comp.prob1a 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.prob1 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!
poster.comp.prob2a 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.prob2 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!
freqK_bfrequency table for the MCMC sample of the number of mixture components in the distribution of the random effects.
propK_bposterior probabilities for the numbers of mixture components in the distribution of random effects.
poster.mean.ya 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.
poster.mean.profilea 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.
poster.mean.w_ba 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.
poster.mean.mu_ba 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.
poster.mean.Q_ba 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.
poster.mean.Sigma_ba 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.
poster.mean.Li_ba 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.
Deviancenumeric 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.
Cond.Deviancenumeric vector with a chain for the conditional deviances, i.e., twice the conditional (given random effects) log-likelihoods.
K_bnumeric vector with a chain for $K_b$ (number of mixture components in the distribution of random effects).
w_bnumeric 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.
mu_bnumeric 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.
Q_bnumeric 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.
Sigma_bnumeric 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.
Li_bnumeric 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.
gammaInv_bmatrix with $dimb$ columns with a chain for inverses of the hyperparameter $\boldsymbol{\gamma}_b$.
order_bnumeric 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.
rank_bnumeric 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.
mixture_bdata.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.
ba matrix with the MCMC chains for random effects. It is included only if store[b] is TRUE.
alphanumeric vector or matrix with the MCMC chain(s) for fixed effects.
sigma_epsnumeric vector or matrix with the MCMC chain(s) for standard deviations of the error terms in the (mixed) models for continuous responses.
gammaInv_epsmatrix with $dimb$ columns with MCMC chain(s) for inverses of the hyperparameter $\boldsymbol{\gamma}_b$.
relabel_ba list which specifies the algorithm used to re-label the MCMC output to compute order_b, rank_b, poster.comp.prob1, poster.comp.prob2, poster.mean.w_b, poster.mean.mu_b, poster.mean.Q_b, poster.mean.Sigma_b, poster.mean.Li_b.
Cpara list with components useful to call underlying C++ functions (not interesting for ordinary users).

Details

See accompanying papers (Komárek et al., 2010, Komárek and Komárková, 2011).

References

Komárek, A. and Komárková, L. (2011). Clustering for Multivariate Continuous and Discrete Longitudinal Data. Submitted for publication.

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, 3267-3283.

Plummer, M. (2008). Penalized loss functions for Bayesian model comparison. Biostatistics, 9, 523-539.

Examples

Run this code

## See also additional material available in 
## YOUR_R_DIR/library/mixAK/doc/
## or YOUR_R_DIR/site-library/mixAK/doc/
## - files PBCseq.pdf,
##         PBCseq.R
## ==============================================

Run the code above in your browser using DataLab