Fits a generalized linear mixed model (S.J. Welham).
Options
PRINT = string token |
What output to display (model, monitoring, components, vcovariance, means, backmeans, effects, waldtests); default mode, moni, comp, vcov, mean, back, effe |
|---|---|
DISTRIBUTION = string token |
Error distribution (binomial, poisson, normal, gamma, negativebinomial); default bino |
LINK = string token |
Link function (identity, logarithm, logit, reciprocal, probit, complementaryloglog, logratio); default * gives the canonical link |
DISPERSION = scalar |
Value at which to fix the residual variance, if missing the variance is estimated; default 1 |
RANDOM = formula |
Random model excluding bottom stratum; this must be set |
FIXED = formula |
Fixed model; default * |
ABSORB = factor |
Absorbing factor to be used at the REML step of the iterations |
CONSTANT = string token |
Whether to estimate or omit constant term in fixed model (omit, estimate); default esti |
FACTORIAL = scalar |
Limit on number of factors/covariates in a model term; default 3 |
PTERMS = formula |
Formula specifying fixed terms for which means or back-transformed means are to be printed; default * prints all the fixed model terms |
PSE = string token |
Standard errors to print with tables of means (differences, estimates, alldifferences, allestimates, vcovariance); default diff, vcov |
MVINCLUDE = string tokens |
Whether to include units with missing values in the explanatory factors and variates and/or the y-variates (explanatory, yvariate); default * i.e. omit units with missing values in either explanatory factors or variates or y-variates |
MAXCYCLE = scalar |
Maximum number of iterations of the GLMM algorithm; default 20 |
TOLERANCE = scalar |
Convergence criterion for iterative procedure; default 0.0001 |
FMETHOD = string token |
Specifies fitting method (all, fixed): all indicates the method of Schall (1991); fixed indicates the marginal method of Breslow & Clayton (1993) ; default all |
OFFSET = variate |
Variate holding values to be used as an offset on the linear predictor scale; default * |
CADJUST = string token |
What adjustment to make to covariates for the REML analysis (mean, none); default mean |
AGGREGATION = scalar |
Fixed parameter for negative binomial distribution (parameter k as in variance function var = mean + mean2/k); default 1 |
KLOGRATIO = scalar |
Parameter k for logratio link, in form log(mean / (mean + k)); default as set in AGGREGATION option |
OWNDIST = text |
For non-standard distributions only: text specifying the variance function to be used with dummy variable DUM, e.g. OWNDIST='DUM' |
OWNLINK = text |
For non-standard link functions only: text specifying 3 functions using dummy variable DUM – the link function, its inverse and its derivative, e.g. OWNLINK = !T('log(DUM)','exp(DUM)','1/DUM') |
CDEFINITIONS = text |
Statements to execute to define correlation models; default * i.e. none |
CVECTORS = pointer |
Data structures involved in the correlation models |
WORKSPACE = scalar |
Number of blocks of internal memory to be set up for use by the REML algorithm; default 1 |
Parameters
Y = variates |
Dependent variates |
|---|---|
NBINOMIAL = scalars or variates |
Number of binomial trials for each unit (must be set if DISTRIBUTION=binomial) |
FITTEDVALUES = variates |
Variates to save fitted values |
COMPONENTS = variates |
Variate to save estimated variance components |
VCOVARIANCE = symmetric matrices |
Variance-covariance matrix for the variance components |
MEANS = pointers |
Pointer to save tables of means for each Y variate |
VARMEANS = pointers |
Pointer to save covariance matrices of tables of means for each Y variate |
BACKMEANS = pointers |
Pointer to save tables of back-transformed means for each Y variate |
ITERATIVEWEIGHTS = variates |
Saves the iterative weights from the generalized linear model fitting |
INITIALFITTEDVALUES = variates |
Defines initial values for the fitted values; if unset, these are formed automatically |
SAVE = REML save structures |
Saves details of the REML analysis used to fit the model |
Description
Procedure GLMM estimates the parameters of a generalized linear mixed model using either the method of Schall (1991) or the marginal method of Breslow & Clayton (1993), as described in the Methods Section.
The procedure assumes a generalized linear mixed model, that is a generalized linear model with both fixed and Normally-distributed random effects on the scale of the linear predictor. The procedure estimates the fixed effects together with the variance components associated with the random effects.
The DISTRIBUTION option sets the error distribution; the default is to assume a binomial distribution but the poisson, gamma and negative-binomial distributions are also available. Other distributions can be used via the OWNDIST option; this should be set to a text containing the formula for calculating the variance function for the required distribution, in terms of dummy variable DUM. The link can be set using the LINK option; the default takes the canonical link. Identity, logarithm, logit, reciprocal, probit, complementaryloglog or logratio link functions are also provided, and alternative link functions can be used via the OWNLINK option. In this case, OWNLINK must be set to a text with three values containing formulae (in terms of dummy variable DUM) for calculating the link function, its inverse and its first derivative. For example, instead of specifying a Poisson distribution with log link, the OWNDIST and OWNLINK options could be set as
OWNDIST='DUM'; OWNLINK=!T(LOG(DUM),EXP(DUM),'1/DUM')Where necessary, these expressions should be constructed so that invalid results (eg. divide by zero or log(zero)) are avoided.
The AGGREGATION option supplies the aggregation parameter for the negative-binomial distribution; default 1. The KLOGRATIO option supplies the parameter k to be used in the logratio link, and takes its default from AGGREGATION.
The dispersion parameter is assumed to be 1 unless otherwise specified by the DISPERSION option. Setting DISPERSION=* requests that the dispersion parameter be estimated.
The fixed and random models are specified by the FIXED and RANDOM options. The number of factors in the terms of the fixed model can be limited using the FACTORIAL option. The ABSORB option can specify an absorbing factor for use in the REML steps of the GLMM algorithm. However, if the absorbing factor appears in any of the terms of the FIXED model, no estimates of error will be available for these terms (see the Guide to the Genstat Command Language, Part 2, Sections 5.3.3 and 5.3.7). By default, a constant term is included in the model; this can be suppressed by setting option CONSTANT=omit. An offset can be included in the linear predictor by setting option OFFSET. By default any covariates are centred for the REML fitting by subtracting their means, weighted according to the iterative weights of the generalized linear model. You can save the iterative weights using the ITERATIVEWEIGHTS parameter, or you can set option CADJUST=none to request that the uncentred covariates are used instead.
It is also possible to define correlation models on the random terms, although the results should be used with caution as their properties are not yet well understood. To do this, you should set the CDEFINITIONS option to a text containing the Genstat statements required to define the models (e.g. using VSTRUCTURE). You also need to set the CVECTORS option to a pointer containing the data structures involved in the statements. Then, in the statements themselves, you should refer to each of these as CVECTORS[n], where n is the position of the relevant data structure in the pointer. For example:
TEXT cdef; VALUE=\
'VSTRUCTURE [CVECTORS[1].CVECTORS[2]] ar,ar; FACTOR=CVECTORS[1,2]; ORDER=1'
GLMM [DISTRIBUTION=gamma; LINK=log; FIXED=variety;\
RANDOM=fieldrow*fieldcolumn; CDEFINITION=cdef;\
CVECTORS=!p(fieldrow,fieldcolumn)] yield
The MVINCLUDE option allows the inclusion of units with missing values, as in the REML directive. By default, units where there is a missing value in the y-variate or in any of the factors or variates in the model terms are excluded. The setting explanatory allows units with missing values in factors or variates in the model to be included. For missing covariate values, this is equivalent to substituting the mean value. The setting yvariate includes units with missing values in the y-variate. This can be useful to retain the balanced structure of the data for use with direct product covariance matrices (see VSTRUCTURE), or to produce predictions of data values for given values of explanatory factors and/or variates.
The FMETHOD option specifies the method used to form the fitted values and therefore determines the fitting method to be used. The default setting all specifies that both fixed and random terms should be used to form fitted values which gives the method of Schall (1991); setting fixed indicates that only fixed terms are used to form fitted values which gives the marginal method of Breslow & Clayton (1993).
Output is controlled by options PRINT, PTERMS and PSE. PRINT allows printing of the current model, monitoring information, estimates of the variance components, their variance-covariance matrix, Wald tests, tables of means on the scale of the linear predictor (with standard errors), tables of back-transformed means (i.e. on the original scale) and tables of effects. If there is an offset, the predicted means are for an offset value of zero. Option PTERMS can select which tables of fixed effect means are to be printed; by default, tables of means are produced for all the terms in the fixed model. Option PSE controls the standard errors that are printed with tables of means: differences produces a summary of standard errors of differences between means; estimates produces a summary of standard errors of the means; allestimates produces a standard error for every mean; vcovariance produces the variance-covariance matrix for the table; alldifferences produces the full matrix of standard errors of differences between means. Setting PSE=* alone suppresses printing of error estimates. More than one setting can be used and, by default, a summary of seds and the variance covariance matrix are printed for each table.
Some control over the iterative GLMM algorithm is provided by option MAXCYCLE which sets the maximum number of iterations (default 20), and by option TOLERANCE which specifies the criterion for determining convergence of the algorithm (default 0.0001). Convergence is judged to have been attained once the maximum change in the ratio (variance component)/(residual variance) and the change in the residual variance are less than the specified TOLERANCE.
The dependent variate is specified using the Y parameter. The NBINOMIAL parameter must be set when DISTRIBUTION=binomial to specify the total number of trials on each unit, as a variate if the number varies from unit to unit or as a scalar if it is constant over all the units.
The other parameters are used to save results. The variance components and residual variance can be saved in a variate using parameter VCOMPONENTS, with their variance-covariance matrix stored in a symmetric matrix specified by parameter VCOVARIANCE. The tables of means to be saved are determined by the setting of PTERMS. The tables are stored in a pointer specified by parameter MEANS, in the order in which they appear in the FIXED model. Their variance matrices and tables of back-transformed means are stored similarly in pointers specified by parameters VARMEANS and BACKMEANS.
VDISPLAY and VKEEP can be used after procedure GLMM to redisplay or store other results from the internal REML estimation. You can use the SAVE parameter to save the associated REML save structure, so that the information will still be available if REML is used for another analysis in the interim.
Options: PRINT, DISTRIBUTION, LINK, DISPERSION, RANDOM, FIXED, ABSORB, CONSTANT, FACTORIAL, PTERMS, PSE, MVINCLUDE, MAXCYCLE, TOLERANCE, FMETHOD, OFFSET, CADJUST, AGGREGATION, KLOGRATIO, OWNDIST, OWNLINK, CDEFINITIONS, CVECTORS, WORKSPACE.
Parameters: Y, NBINOMIAL, FITTEDVALUES, COMPONENTS, VCOVARIANCE, MEANS, VARMEANS, BACKMEANS, ITERATIVEWEIGHTS, INITIALFITTEDVALUES, SAVE.
Method
GLMM estimates the parameters of the Generalized Linear Mixed Model using either the method of Schall (1991) or the marginal method of Breslow & Clayton (1993). The method used is determined by the setting of option FMETHOD.
The data y arises from some specified distribution with variance function sV and expected value μ. The link function g (with inverse h) is such that
g(μ) = η = X a + Z b
where X is the design matrix for the vector a of fixed effects and Z is the design matrix for the vector b of random effects. The random effects b can be attributed to c random factors which are assumed to have zero mean and to be uncorrelated with each other and with e:
Cov(b) = D = Diag{ σ12 I1 … σc2 Ic }
The method used by Schall (1991) develops the algorithm by analogy with the algorithm for estimating conventional generalized linear models. The link function applied to the data is linearized about μ to give the adjusted dependent variate z,
z = X a + Z b + e g′(μ)
where e=y-μ and g′ = dg/dμ.
Then
E(z) = X a; Cov(b) = D;
Cov(e g′(μ)) = sV(μ) × (dη/dμ) × (dη/dμ) = s × W(μ)-1
where s is the dispersion parameter. Hence
Cov(z) = s × W(μ)-1 + Z D Z′.
This has the same form as the general linear mixed model, and the fixed effects and variance components can be estimated by REML with (iterative) weights W.
This leads to the following algorithm:
Step 1) Using initial estimates of the variance components and of μ, calculate the adjusted variate z and weights W.
Step 2) Get new estimates of the variance components and of μ by REML on adjusted variate z with weights W.
Step 3) Convergence in estimates ⇒ exit algorithm.
Step 4) Use new estimates to update adjusted variate z and weights vector W.
Step 5) Go to Step 2.
The marginal model used by Breslow and Clayton is derived from a first order approximation (linearisation about Xa) to give
y ∼ h(Xa) + h′(Xa)Zb + e
where ∼ indicates approximation, h is the inverse of the link fuction g and e is y-μ. They then work in terms of the marginal mean, M=h(Xa). Quasi-likelihood estimation leads to an algorithm similar to the one above, but the working variate becomes
z = Xa + (y–M)g′(M) = Xa + Eg′(M)
where E=y–M. The working variate z then has variance
Cov(z) = s × W(M)-1 + Z D Z′.
The same algorithm is used to fit the model, replacing μ by M and e by E.
The only difference between the two algorithms is then in the method used to form the mean μ or M and the “error” variate e or E. The option RMETHOD of REML controls the method of forming fitted values after REML estimation (i.e. including just fixed terms, or all terms except the residual) and this option is used inside the procedure to determine which of the models is fitted.
Initial values for the variance components are calculated by REML estimation using the fixed and random models on the data transformed by the link function. Initial values for the fixed effects are calculated by fitting the fixed model only to a generalized linear model with the specified link and error distribution. The WORKSPACE option specifies the number of blocks of internal memory to be set up for use by the REML algorithm; see the REML directive for more details.
Action with RESTRICT
If the Y-variate is restricted, only the units not excluded by the restriction will be analysed.
References
Breslow, N.E. & Clayton, D.G. (1993). Approximate inference in generalized linear mixed models. Journal of the American Statistical Association, 88, 421, 9-25.
McCullagh, P. & Nelder, J.A. (1989). Generalized Linear Models (second edition). Chapman & Hall, London.
Schall, R. (1991) Estimation in generalized linear models with random effects. Biometrika, 78, 719-727.
See also
Commands for: Regression analysis.
Example
CAPTION 'GLMM example',\
!t('Data from McCullagh & Nelder (1989, Table 14.4),',\
'also see Schall (1991).'); STYLE=meta,plain
FACTOR [NVALUES=120; LEVELS=20] Female, Male
& [LEVELS=4; LABELS=!t(RR,RW,WR,WW)] Cross
VARIATE [NVALUES=120] Mate1
READ Cross,Male,Female; FREPRESENTATION=labels,2(levels)
RR 1 1 RW 14 1 RR 5 1 RW 11 1 RR 4 1 RW 15 1 RR 5 2 RW 15 2
RR 3 2 RW 13 2 RR 1 2 RW 12 2 RR 2 3 RW 11 3 RR 1 3 RW 14 3
RR 3 3 RW 13 3 RR 4 4 RW 12 4 RR 2 4 RW 15 4 RR 5 4 RW 14 4
RR 3 5 RW 13 5 RR 4 5 RW 12 5 RR 2 5 RW 11 5 RW 19 6 RR 9 6
RW 20 6 RR 7 6 RW 16 6 RR 8 6 RW 18 7 RR 8 7 RW 19 7 RR 9 7
RW 17 7 RR 6 7 RW 16 8 RR 6 8 RW 17 8 RR 10 8 RW 20 8 RR 9 8
RW 20 9 RR 7 9 RW 18 9 RR 6 9 RW 19 9 RR 10 9 RW 17 10 RR 10 10
RW 16 10 RR 8 10 RW 18 10 RR 7 10 WR 9 11 WW 19 11 WR 7 11 WW 20 11
WR 10 11 WW 18 11 WR 7 12 WW 16 12 WR 9 12 WW 17 12 WR 6 12 WW 20 12
WR 8 13 WW 17 13 WR 6 13 WW 19 13 WR 7 13 WW 16 13 WR 10 14 WW 20 14
WR 8 14 WW 18 14 WR 9 14 WW 19 14 WR 6 15 WW 18 15 WR 10 15 WW 16 15
WR 8 15 WW 17 15 WW 15 16 WR 2 16 WW 13 16 WR 4 16 WW 12 16 WR 1 16
WW 14 17 WR 1 17 WW 15 17 WR 2 17 WW 11 17 WR 5 17 WW 11 18 WR 4 18
WW 12 18 WR 5 18 WW 15 18 WR 3 18 WW 13 19 WR 3 19 WW 11 19 WR 1 19
WW 14 19 WR 4 19 WW 12 20 WR 5 20 WW 14 20 WR 3 20 WW 13 20 WR 2 20:
READ Mate1
1 1 1 0 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1
1 1 1 0 1 1 0 0 0 1 0 0 0 1 0 0 1 1 0 0 1 1 1 1 0 0 1 0 1 0
0 1 1 1 0 1 0 0 0 1 0 0 0 0 0 0 0 1 0 1 1 1 0 1 0 1 0 0 0 0
0 0 1 0 0 0 1 1 1 0 1 1 1 0 1 0 1 0 1 1 1 1 1 0 1 0 0 1 1 0 :
GLMM [DISTRIBUTION=binomial; LINK=logit; FIXED=Cross; RANDOM=Female+Male]\
Mate1; NBINOMIAL=1