cumlogitRE {glmmAK} | R Documentation

## Logit and cumulative logit model with random effects

### Description

This function implements MCMC sampling for the logit model with binary response and the cumulative logit model for multinomial ordinal response. Details are given in Komárek and Lesaffre (2007). On as many places as possible, the same notation as in this paper is used also in this manual page.

In general, the following (cumulative) logit model for response Y is assumed:

log[P(Y>=1)/P(Y=0)] | = | eta[1]
log[P(Y>=2)/P(Y<=1)] | = | eta[2]
| ... |
log[P(Y=C)/P(Y<=C-1)] | = | eta[C],

where the form of the linear predictors eta[1],...,eta[C] depends on whether a hierarchical centering is used or not. In the following, beta denotes fixed effects and b random effects.

No hierarchical centering (DEFAULT)
The linear predictor for the cth logit has the following form

eta[c] = beta[c]'(v', v(b)') + beta[*]'(x', x(b)') + b'(v(b)', x(b)') (c=1,...,C),

where beta=(beta[1]',...,beta[C]',beta[*])' is the vector or the fixed-effects and b is a vector of random effects with zero location.

Hierarchical centering
The linear predictor for the cth logit has the following form

eta[c] = beta[c]'v + beta[*]'x + b[c]'v(b) + b[*]x(b) (c=1,...,C),

where beta=(beta[1]',...,beta[C]',beta[*])' is the vector or the fixed-effects and b=(b[1]',...,b[C]',b[*]')' is a vector of random effects with location alpha=(alpha[1]',...,alpha[C]',alpha[*])'.

Normal random effects (`drandom="normal"`)
A vector of random effects is assumed to follow a (multivariate) normal distribution.

That is, if there is no hierarchical centering we assume for b:

b ~ N(0, D[b]),

where D[b] is their variance-covariance matrix.

If the random effects are hierarchically centered then we assume for b=(b[1]',...,b[C]',b[*]')':

b ~ N(alpha, D[b]),

where alpha=(alpha[1]',...,alpha[C]',alpha[*])' is a vector of random effect locations (means) and D[b] is their variance-covariance matrix.

Further, in a Bayesian model, it is assumed that D[b]^(-1) has a Wishart W(nu[b], S[b]) prior with nu[b] degrees of freedom and scale matrix S[b]. That is, a priori

D[b]^(-1) ~ W(nu[b], S[b]),

E(D[b]^(-1)) = nu[b]*S[b].

Note that nu[b] must be higher than the number of random effects minus 1.

Alternatively, when there is only a univariate random effect with the variance d[b]^2, it is possible to specify a uniform prior for the standard deviation of the random effect. That is, a priori

d[b] ~ Unif(0, S).

G-spline distributed random effects (`drandom="gspline"`)
See `gspline1` and `gspline2` for description of the G-spline (penalized Gaussian mixture) distribution. Further, see Komárek and Lesaffre (2007) for description of the prior distribution on the G-spline. Brief description follows here as well.

Univariate G-spline
If there is only a univariate random effect in the model and its distribution is specified as a G-spline that it is assumed that

b ~ alpha + sum[j=-K]^K w[j] N(tau*mu[j], (tau*sigma)^2),

where alpha is a location parameter, tau is a scale parameter and w=(w[-K],...,w[K])' is a vector of G-spline weights. For hierarchically centered random effects, the location parameter alpha is fixed to zero.

Further, M=(mu[-K],...,mu[K])' is a vector of fixed equidistant knots (component means), where

mu[j] = j*delta (j=-K,...,K)

and sigma is a fixed basis standard deviation.

The constraints 0<w[j]<1 (j=-K,...,K) and sum[j=-K]^K w[j] = 1 are sufficient for G-spline to be a density. To avoid constraint estimation we will estimate transformed weights a=(a[-K],...,a[K])' instead which relates to the original weights by

w[j] = exp(a[j])/sum[k=-K]^K exp(a[k]) (j=-K,...,K).

a[j] = log(w[j]/w[0]) (j=-K,...,K).

In the estimation procedure a penalty on the a coefficients in the form of a Gaussian Markov random field prior is imposed.

Bivariate G-spline
If there is a bivariate random effect b=(b[1], b[2])' in the model and its distribution is specified as a G-spline that it is assumed that

b ~ (alpha[1], alpha[2])' + sum[j1=-K1]^K1 sum[j2=-K2]^K2 w[j1,j2] N((tau[1]*mu[1,j1], tau[2]*mu[2,j2])', diag((tau[1]*sigma[1])^2, (tau[2]*sigma[2])^2)),

where alpha[1], alpha[2] are location parameters, tau[1], tau[2] are scale parameters and W=(w[-K1,-K2],...,w[K1,K2])' is a matrix of G-spline weights. For hierarchically centered random effects, the location parameters alpha[1], alpha[2] are fixed to zero.

Further, M[1]=(mu[1,-K1],...,mu[1,K1])' is a vector of fixed equidistant knots (component means) in the first margin, where

mu[1,j1] = j1*delta[1] (j1=-K1,...,K1)

and sigma_1 is a fixed basis standard deviation in the first margin. Similarly for the second margin.

Similar reparametrization of the G-spline weights as in the univariate case is used to avoid constrained estimation.

### Usage

    cumlogitRE(y, v, x, vb, xb, cluster,
          intcpt.random=FALSE,
          hierar.center=FALSE,
          drandom=c("normal", "gspline"),
          C=1,
          logit.order=c("decreasing", "increasing"),
          prior.fixed,
          prior.random,
          prior.gspline,
          init.fixed,
          init.random,
          init.gspline,
          nsimul = list(niter=10, nthin=1, nburn=0, nwrite=10),
          store = list(prob=FALSE, b=FALSE, alloc=FALSE, acoef=FALSE),
          dir=getwd(),
          precision=8)

### Arguments

`y` | response vector taking values 0,1,...,C.
`v` | vector, matrix or data.frame with covarites for fixed effects whose effect does not necessarily satisfy proportional odds assumption. If the argument `intcpt.random` is set to `FALSE` then the fixed intercept is included by default in the model. The intercept column should be included neither in `x`, nor in `v`.
`x` | vector, matrix or data.frame with covarites for fixed effects whose effect is assumed to satisfy proportional odds assumption.
`vb` | vector, matrix or data.frame with covarites for random effects whose effect does not necessarily satisfy proportional odds assumption. If you want to include random intercept, do it by setting the argument `intcpt.random` to `TRUE`. The intercept column should be included neither in `xb`, nor in `vb`.
`xb` | vector, matrix or data.frame with covarites for random effects whose effect is assumed to satisfy proportional odds assumption.
`cluster` | vector which determines clusters. Needed only when there are any random effects in the model.
`intcpt.random` | logical indicating whether a random intercept should be included in the model.
`hierar.center` | logical indicating whether a hierarchical centering of random effects should be used or not.
`drandom` | character indicating assumed distribution of random effects (if there are any).
`C` | number of response categories minus 1.
`logit.order` | either "decreasing" or "increasing" indicating in which direction the logits are formed. See the same argument in `cumlogit` for more details. Currently, only "decreasing" is implemented for `cumlogitRE`.
`prior.fixed` | list specifying the prior distribution for the fixed effects (regression coefficients). mean vector giving the prior mean for each regression coefficient. It can be a single number only in which case it is recycled. var vector giving the prior variance for each regression coefficient. It can be a single number only in which case it is recycled. In the case that prior mean and/or variances are different for different regression coefficients then they should be given in the following order: intercept for the first logit, beta(`v`) for the first logit, ..., intercept for the last logit, beta(`v`) for the last logit, beta(`x`),
`prior.random` | list specifying the prior distribution for the parameters of the distribution of random effects. Composition of this list depends on the chosen distribution of random effects (normal or Gspline). Mdistrib character specifying the prior distribution of the location of random effects. It is ignored if `hierar.center`=`FALSE`, in which case the location of the random effects is fixed to zero. It can be one of the following. “`fixed`” Locations of random effects are assumed to be fixed and are not updated. “`normal`” Locations of random effects are assumed to be apriori normally distributed. Parameters of the normal distribution are specified further by items `Mmean` and `Mvar`. This is also a default choice when `Mdistrib` is not specified. Mmean vector giving the prior means for the locations of random effects. It can be a single number only in which case it is recycled. It is ignored when `hierar.center` is `FALSE` in which case all random effects have zero location. Mvar vector giving the prior variances for the locations of random effects. It can be a single number only in which case it is recycled. In the case that prior means and/or variances are different for different locations of random effects then they should be given in the similar order as specified above for argument `prior.fixed`. It is ignored when `hierar.center` is `FALSE` in which case all random effects have zero location. Ddistrib character specifying the prior distribution of the covariance matrix of random effects. It can be one of the following. “`fixed`” covariance matrix of random effects is assumed to be fixed and is not updated. “`wishart`” inverse of the covariance matrix of normally distributed random effects is assumed to have a priori Wishart distribution. Parameters of the Wishart distribution are specified further by items `Ddf` and `DinvScale`. This is also a default choice when `Ddistrib` is not specified and random effects are normally distributed. This option is not allowed when the random effects distribution is modelled using G-splines. “`sduniform`” when random effects are normally distributed then this option may be used when there is only a univariate random effect in the model. Its standard deviation is then assumed to follow a uniform distribution on the interval (0,S). Value of S is specified further by an item `Dupper`. When random effects distribution is modelled using G-splines then `sduniform` can be used also for multivariate random effects. It is then assumed that the overall standard deviation tau[m] in the mth margin follows a uniform distribution on the interval (0,S[m]). “`gamma`” when random effects are normally distributed then this option may be used when there is only a univariate random effect in the model. Its inverse variance is then assumed to have a Gamma prior with the shape specified further by the item `Dshape` and the rate (inverse scale) specified by the item `DinvScale`. When random effects distribution is modelled using G-splines then `gamma` can be used also for multivariate random effects. It is then assumed that the overall inverse variance tau[m]^[-2] in the mth margin follows a Gamma prior with the shapes specified further by the item `Dshape` and the rates (inverse scales) specified by the item `DinvScale`. Ddf degrees of freedom nu[b] for the Wishart prior distribution of the inverse covariance matrix of normally distributed random effects. Dshape number or vector with shape parameters for the gamma priors of the variance components of the random effects. If it is a single number and random effects are multivariate it may be recycled. DinvScale number or matrix determining the inverse scale matrix S[b]^(-1) for the Wishart prior distribution of the inverse covariance matrix of normally distributed random effects. Or number or vector giving the rate (inverse scale) parameter(s) for the gamma priors of the variance components. If it is a single number and `Ddistrib` is wishart then it is assumed that S[b]^(-1) is diagonal with that single number on a diagonal. Dupper upper limit for the uniform distribution of the standard deviation of the random effect(s) when `Ddistrib` is equal to sduniform. If it is a single number and random effects are multivariate it may be recycled.
`prior.gspline` | list specifying the G-spline distribution of random effects and prior distribution of the G-spline parameters. This argument is required only when `drandom` is equal to gspline. In the following let q denote the number (dimension) of random effects. The list `prior.gspline` can have the following components. K vector of length q or a number (it is recycled) which specifies, for each marginal G-spline, the number of knots on each side of the zero knot. That is, the m-th marginal G-spline has 2K[m]+1 knots. It is set to 15 if not explicitely specified. delta vector of length q or a number (it is recycled) which specifies the distance between two consecutive knots for each marginal G-spline. That is, the m-th marginal G-spline has the following knots mu[m,j] = j*delta[m], j=-K[m],...,K[m]. It is set to 0.3 if not explicitely specified. sigma vector of length q or a number (it is recycled) which specifies the basis standard deviation of each marginal G-spline. sigma[m] is set to (2/3)*delta[m] if not explicitely specified. CARorder vector of length q or a number (it is recycled) giving the order of the intrinsic conditional autoregression in the Gaussian Markov random field prior for the transformed G-spline weights in each margin. It does not need to be specified when `neighbor.system` is different from `uniCAR`. It is set to 3 if not explicitely specified. neighbor.system character specifying the type of the Gaussian Markov random field in the prior for the transformed G-spline weights a of a bivariate G-spline. It does not have to be specified for univariate G-splines. It can be one of the following. “`uniCAR`” univariate (in each margin) conditional autoregression as described in Komárek and Lesaffre (2007). That is a priori p(a|lambda) propto exp[ -{ lambda[1]/2 sum[j1]...sum[jq](Delta[1]^d a[j1,...,jq])^2 + ... + lambda[q]/2 sum[j1]...sum[jq](Delta[q]^d a[j1,...,jq])^2 } ], where Delta[m]^d is a difference operator of order d in the mth margin, e.g., Delta[1]^3 a[j1,j2,...,jq] = a[j1,j2,...,jq] - 3a[j1-1,j2,...,jq] + 3a[j1-2,j2,...,jq] - a[j1-3,j2,...,jq], and lambda = (lambda[1], ..., lambda[q])' are smoothing hyperparameters. This is also a default choice when `neighbor.system` is not specified. “`eight.neighbors`” this prior is applicable for bivariate G-splines only and is based on eight nearest neighbors in a spatial meaning. That is, except on edges, each full conditional of a depends only on eight nearest neighbors and local quadratic smoothing. The prior is then defined as p(a|lambda) propto exp(-lambda/2 * sum[j1=-K1][K1-1]sum[j2=-K2][K2-1] (Delta a[j1,j2])^2), where Delta a[j1,j2] = a[j1,j2] - a[j1+1,j2] - a[j1,j2+1] + a[j1+1,j2+1]. Parameter lambda is a common smoothing hyperparameter. “`twelve.neighbors`” not (yet) implemented Ldistrib character specifying the prior distribution of the smoothing hyperparameters lambda[m], m=1,...,q (precision parameters of the Markov random fields in each margin) or of a common hyperparameter lambda It can be one of the following. “`fixed`” smoothing hyperparameters lambda are fixed to their initial values and are not updated. “`gamma`” each of the smoothing hyperparameters lambda[1],...,lambda[q] is assumed to follow a Gamma prior with the shapes specified further by the item `Lshape` and the rates (inverse scales) specified further by the item `LinvScale`. This is also a default choice when `Ldistrib` is not specified. “`sduniform`” square root of the inversion of each smoothing hyperparameter, i.e., sqrt(lambda[1]^{-1}),...,sqrt(lambda[q]^{-1}) is assumed to follow apriori a uniform distribution on the intervals (0, S[m,lambda]). Values of S[1,lambda],...,S[q,lambda] are specified further by the item `Lupper`. Lequal logical indicating whether all smoothing hyperparameters should be kept equal. It is set to `FALSE` if not explicitely specified and `neighbor.system` is `uniCAR`. It is always `TRUE` when `neighbor.system` is different from `uniCAR`. Lshape number or vector with shape parameters for the gamma priors of the smoothing hyperparameters lambda. If it is a single number and there is more than one smoothing hyperparameter lambda in the model it may be recycled. LinvScale number or vector with rate (inverse scale) parameters for the gamma priors of the smoothing hyperparameters lambda. If it is a single number and there is more than one smoothing hyperparameter lambda in the model it may be recycled. Lupper number or vector with upper limits for the uniform distribution on sqrt(lambda^{-1}) parameters when `Ldistrib` is sduniform. If it is a single number and there is more than one smoothing hyperparameter lambda in the model it may be recycled. Aident character specifying in which way the transformed G-spline weights (a coefficients) are identified. It can be one of the following. “`mean`” with this option, the a coefficients are forced to sum up to zero and have a zero mean. Note that this option usually causes problems during MCMC, especially with bivariate G-splines. The reason is that if there are many almost zero weights, they lead to many negative a coefficients and to satisfy the zero mean constrain, there must be some a's which are highly positive. When exponentiating them to get weights w, an overflow occur. “`reference`” with this option, one of the a coefficients in each margin is chosen as the reference one and is always equal to zero. Index of the reference coefficient is specified by the item `Areference` (see below). This is also a default choice when `Aident` is not specified. Areference vector or number (it is recycled) which specifies the index of the reference a coefficient in each margin in the case `Aident` is equal to reference. For the m-th margin, it must be an integer between -K[m],...,K[m]. To avoid numerical problems, the index of the reference a coefficient may change during the MCMC. AtypeUpdate character specifying in which way the transformed G-spline weights (a coefficients) are updated. It can be one of the following. “`slice`” slice sampler of Neal (2003). “`ars.quantile`” adaptive rejection sampling of Gilks and Wild (1992) with starting abscissae being quantiles of the envelop at the previous iteration. “`ars.mode`” adaptive rejection sampling of Gilks and Wild (1992) with starting abscissae being the mode plus/minus 3 times estimated standard deviation of the full conditional distribution. “`block`” all a coefficients are updated in 1 block using the Metropolis-Hastings algorithm. This is only available for univariate G-splines. Default is `slice`.
`init.fixed` | optional vector with initial values for fixed effects, supplied in the same order as described above in the argument `prior.fixed`. If not given, initials are determined from the maximum-likelihood fit to the model where random effects are replaced by corresponding fixed effects.
`init.random` | optional list with initial values for the random effects and parameters determining their distribution. It can have the following components. b vector or matrix with initial values of cluster specific random effects. Number of rows of the matrix or the length of the vector must be equal to `length(unique(cluster))`. Columns of the matrix correspond to the random effects in the same order as described above for argument `prior.random`. If not given, initials are taken to be equal over clusters and equal to the corresponding fixed effects from the maximum-likelihood fit to the model with fixed effects only. mean vector giving the initial values of the means of random effects. If not given, initials are taken to be equal to the corresponding fixed effects from the maximum-likelihood fit to the model with fixed effects only. var matrix giving the initial value for the covariance matrix of the random effects when `drandom` is normal. Vector giving the initial values of marginal overall variances of the random effects when `drandom` is gspline.
`init.gspline` | optional list with initial values related to the G-spline distribution of random effects (if `drandom` is equal to gspline). It can have the following components. lambda vector with initial values of the smoothing hyperparameters (precision parameters of the Markov random field). If not fixed and not given, the initials are sampled from the prior distribution. weights for univariate G-splines: vector with initial weights. For bivariate G-splines: matrix with initial weights. If the initial weights do not sum up to 1 they are re-scaled. alloc vector or matrix with initial component allocations for the individual random effects. For univariate random effects this should be a vector with numbers from {-K,...,K}. For bivariate random effects this should be a matrix with two columns where in the first column numbers from {-K1,...,K1} appear and in the second column numbers from {-K2,...,K2} appear.
`nsimul` | list indicating the length of the MCMC. It should have the following components. niter total number of the MCMC iterations after discarding the thinned values nthin thinning of the sample nburn length of the burn-in period nwrite frequency with which the iteration count changes. Further, during the burn-in, only every `nwrite`th sampled value is stored on the disk
`store` | list indicating which chains (out of these not stored by default) should be compulsory stored. The list has the logical components with the following names. prob if `TRUE` values of individual predictive probabilities are stored. b if `TRUE` values of cluster specific random effects are stored. alloc if `TRUE` values of allocation indicators are stored. acoef if `TRUE` and distribution of random effects is given as a bivariate G-spline values of log-G-spline weights (a coefficients) are stored for all components.
`dir` | character string specifying the directory in which the sampled values are stored.
`precision` | precision with which the sampled values are written in files.

### Value

The function returns a complete list of parameters of the prior distribution and initial values.
The main task of this function is to sample from the posterior distribution using MCMC. Sampled values are stored in various files which are described below.

### Files created

iteration.sim
  one column labeled `iteration` with indeces of MCMC iterations to which the stored sampled values correspond.
betaF.sim
  sampled values of the fixed effects beta=(beta[1]',...,beta[C]',beta[*])'.

  Note that in models with G-spline distributed random effects which are not hierarchically centered, the average effect of the covariates involved in the random effects (needed for inference) is obtained as a sum of the corresponding beta coefficient and a scaled mean of the G-spline. beta coefficients adjusted in this way are stored in the file ‘betaRadj.sim’ (see below).

betaR.sim
  sampled values of the location parameters alpha=(alpha[1]',...,alpha[C]',alpha[*])' of the random effects when the hierarchical centering was used.

  Note that in models with G-spline distributed random effects which are hierarchically centered, the average effect of the covariates involved in the random effects (needed for inference) is obtained as a sum of the corresponding alpha coefficient and a mean of the G-spline. alpha coefficients adjusted in this way are stored in the file ‘betaRadj.sim’ (see below).

varR.sim
  variance components of the random effects. Format of the file depends on the assumed distribution of the random effects.
  Normal random effects
    Let q be the dimension of the random effects. The first 0.5q(q+1) columns of ‘varR.sim’ contain a lower triangle (in column major order) of the matrix D[b], the second 0.5q(q+1) columns of ‘varR.sim’ contain a lower triangle of the matrix D[b]^(-1).
  Univariate G-spline random effects
    The first column of ‘varR.sim’ contains the G-spline variance parameter tau^2, the second column its inverse.
  Bivariate G-spline random effects
    The first two columns of ‘varR.sim’ contain the G-spline variance parameters tau[1]^2, tau[2]^2, the second two columns their inverse.
loglik.sim
  sampled values of the log-likelihood (conditioned by the values of random effects).
probability.sim
  sampled values of category probabilities for each observations.

  Created only if `store$prob` is `TRUE`.

b.sim
  sampled values of individual random effects.

  Stores complete chains only if `store$b` is `TRUE`.

### Files created for models with G-spline distributed random effects

gspline.sim
  information concerning the fixed parameters of the G-spline which includes: dimension q of the G-spline, numbers of knots on each side of the reference knot for each margin (K[1],...,K[q]), basis standard deviations sigma[1],...,sigma[q] for each margin and knots mu[1,-K[1]],...,mu[1,K[1]], ..., mu[q,-K[q]],...,mu[q,K[q]] for each margin.
weight.sim
  this file is created only for bivariate G-splines and stores the weights w of the G-spline which are higher than a certain threshold value. That is, the weights that are numerically equal to zero are not recorded here. The link between the weights and G-spline components is provided by the file ‘knotInd.sim’.
knotInd.sim
  this file is created only for bivariate G-splines. In its first column, it stores the number of G-spline components for which the weights are recorded on a corresponding row of the file ‘weight.sim’. Subsequently, it stores indeces of the G-spline components for which the weights are given in the file ‘weight.sim’. The indeces are stored as single indeces on the scale 0,...,(2K[1]+1)*(2K[2]+1)-1 such that index 0 corresponds to the component (K[1],K[2]), index 1 corresponds to the component (K[1]+1,K[2]), ..., index K[1]-1 corresponds to the component (K[1],-K[2]), etc.
logweight.sim
  sampled values of the transformed G-spline weights a.

  For univariate G-spline, this file stores always complete chains, irrespective of the value of `store$acoef`. For bivariate G-splines,this file stores complete chains only if `store$acoef` is `TRUE`.

gmoment.sim
  first two moments of the unshifted and unscaled G-spline at each iteration.

  For univariate G-spline the column labeled gmean is equal to

  gmean=sum[j=-K]^K w[j]*mu[j]

  and the column labeled gvar is equal to

  gmean=sum[j=-K]^K w[j]*(mu[j]-gmean)^2 + sigma^2.

  See Komárek, A. and Lesaffre, E. (2007) for formulas that apply in the bivariate case.

  Values stored here are the values of beta[1]^*,...,beta[q]^* and d[1,1]^*, d[2,1]^*, ..., d[q,q]^* as defined in Komárek and Lesaffre (2007).

betaRadj.sim
  sampled values of the average (overall) effects of the random effects. See notes under the files ‘betaF.sim’ and ‘betaR.sim’ above.

  Values stored here are the values of gamma[1],...,gamma[q] as defined in Komárek and Lesaffre (2007).

varRadj.sim
  sampled components of the variance-covariance matrix of the random-effects.

  Values stored here are the values of d[1,1], d[2,1], ..., d[q,q] as defined in Komárek and Lesaffre (2007).

alloc.sim
  sampled values of the component allocations (in Komárek and Lesaffre (2007) denoted by r[i]) for individual random effects.

  For univariate G-spline, the allocations are stored on the scale -K,...,K.

  For bivariate G-spline, the allocation are stored as single indeces on the scale 0,...,(2K[1]+1)*(2K[2]+1)-1 where the link between the single and double indeces is the same as in the file ‘knotInd.sim’.

  Stores complete chains only if `store$alloc` is `TRUE`.

lambda.sim
  sampled values of smoothing hyperparameter(s) lambda.

### Author(s)

Arnošt Komárek arnost.komarek[AT]mff.cuni.cz

### References

Agresti, A. (2002). Categorical Data Analysis. Second edition. Hoboken: John Wiley & Sons.

Gelfand, A. E., Sahu, S. K., and Carlin, B. P. (1995). Efficient parametrisations for normal linear mixed models. Biometrika, 82, 479–488.

Gilks, W. R. and Wild, P. (1992). Adaptive rejection sampling for Gibbs sampling. Applied Statistics, 41, 337–348.

Neal, R. M. (2003). Slice sampling (with Discussion). The Annals of Statistics, 31, 705–767.

Komárek, A. and Lesaffre, E. (2008). Generalized linear mixed model with a penalized Gaussian mixture as a random-effects distribution. Computational Statistics and Data Analysis, 52, 3441–3458.

Molenberghs, G. and Verbeke, G. (2005). Models for Discrete Longitudinal Data. New York: Springer Science+Business Media.

### See Also

`cumlogit`, `logpoissonRE`, `glm`, `polr`.

### Examples

    ### See ex-Toenail.pdf and ex-Toenail.R
    ### available in the documentation
    ### to the package

[Package glmmAK version 1.3-1 Index]