Hierarchal Generalized Linear (Mixed) Models

Basic HGLM

(Double) Hierarchical Generalized Linear (Mixed) Models are implemented in ASReml as a bivariate model where the first variate (Trait 1) is the response used for the Mean model (the trait being analysed) and the second variate (Trait 2) is the response used for the dispersion model and contains the deviance residuals from the mean model. The dispersion is typically modelled as a gamma distribution with log link fitted to the deviance contributions from the first variate.

In the simplest cases, the random effects are modelled with a distribution conjugate to that specified for the primary response variable. The model is specified as for a GLMM model but with !r replaced by !h in the model statement. HGLM is implied if some terms in the model are specified in a !h submodel.

For example

 Angle !GAMMA !LOG !DISP 1 ~ mu Recipe*Temp !h Rep Rep.Batch

models the Rep and Rep.Batch effects as inverse gamma rather than as Normal effects.

The distributions/links available are:

Distribution	Conjugate	Distribution qualifier	Link qualifier
Normal	Normal	!NORMAL	!IDENTITY
Binomial	Beta binomial	!BINOMIAL	!LOGIT
			!PROBIT
			!CLOG
Poisson	Gamma	!POISSON	!LOG
			!SQRT
Gamma	Inverse Gamma	!GAMMA	!INVERSE
			!LOG

Although by default the random effects have the distribution conjugate to that specified for the data, a different distribution/link may be specified using

   !HDIST  Distribution
   !HLINK  Link  qualifiers.

For example

  LM !BIN !PROBIT !TOT Np !DISP 1 !HDIST NORMAL !HLINK IDENTITY) ~ +
      mu Cultivar*Fungicide !h Block Block.Whole

General Cases

More complicated syntax is required to access the full potential of the ASReml implementation. The syntax presented above is for the basic HGLM model where we estimate a common variance within each random term. ASReml sets up the variates and structures and records needed to fit this standard HGLM model. This section shows how to set up the model for more general situations.

The explicit syntax is to specify the data line qualifier !HGLM HSection [terms] and to write the appropriate bivariate model statement using dev(Y) as the name of the second variable (trait 2) where Y is the variate being analysed (trait 1). The second variable is declared for syntax (model parsing) purposes but its values are calculated internally. HSection is a factor which may be defined by the user (a factor in an augmented data file) but if not, will be generated. The method used by Lee and Nelder to fit a HGLM is to augment the data with additional records which model the random effects of the basic model; a set of effects for each term. ASReml uses the same data structure with HSection being a factor which indexes the sections of data: 1 being the original data, 2 the first random factor, 3 the second random factor, etc.. These extra records are created by ASReml if the respective model terms are listed after !HGLM HSection and if HSection is not already defined by the user.

The easiest case is when the conjugate distribution is used for the augmented data. Then you need to specify only the distribution required for the primary data using the GLMM distribution/link qualifiers.

Two additional qualifiers are !HDIST [N|B|G|I] and !HLINK [IDE LOGT LOGN INV] to set the distribution explicitly for the HGLM random effects (i.e. random terms in the mean model which are fitted by augmenting the data). These are data line qualifiers and should be specified along with the !HGLM qualifier.

So, for normal data, the default distribution for the random effects is normal with identity link; for binomial data with logit link, the default distribution for the random effects is beta binomial with logit link; for poisson data with log link, the default distribution for the random effects is gamma with log link; for gamma data with inverse (i.e. reciprocal) link, the default distribution for the random effects is inverse gamma with inverse link; and for gamma data with log link, the default distribution for the random effects is inverse gamma with log link.

The default for the dispersion model ( dev(Y)) is a gamma distribution with log link. This is appropriate for all combinations of distributions and links for the primary variable; it does not need to be changed.

The !HDIST and !HLINK qualifiers allow for other combinations.

Specifying the appropriate model raises a few issues which are discussed by way of example. Consider the Cake example used by Lee and Nelder. The trial measured 'breaking angle' of cakes made using 3 recipes and baked at 6 temperatures. There were 15 replicates consisting of 3 batches, a total of 270 measures. The design is nested with the 3 batches being the 3 recipes within a Replicate and the 6 temperatures nested within Batch/Recipe. The fixed effects model is therefore Recipe*Temperature and the random effects model is Replicates + Rep.Batch + Residual.

The ASReml code is

 Lee and Nelder Cake example;
 #Replicate,Batch,Recipe,Temperature,Angle
 #1,1,1,175,42

  Replicate  *
  Batch  *
  Recipe  *
  Temperature  !I
  Angle  !REP 0 1
 Cake.csv !SKIP 1 !MAXIT 100

 !PART 1  # Simple syntax
 Angle !GAMMA !INV !DISP 1 ~ mu Recipe*Temp !h Rep Rep.Batch

 !PART 2  # Equivalent Long form syntax
 !HGLM HSection Repl Repl.Batch

 Angle !GAMMA !INV dev(Angle)  !DISP 1 ~ +
 at(Trait,1).(at(HSection,1)  Recipe Temp Recipe.Temp),
         at(Trait,1).(Repl.hs(HS,2)  Repl.Batch.hs(HS,3)),
         at(Trait,2).HSection
 residual units.id(2)

Discussion of Code

There are three random components in the model to be fitted: Residual, Repl and Repl.Batch. The algorithm requires augmenting the data, one extra record for each of the 60 extra random effects (15 for Repl and 45 for Repl.Batch). Thus the full data has 3 sections and these are coded in a variable called here HSection. The first section is the 270 records of the original data. The second section is the 15 records for Repl. The third section is the 45 records for Repl.Batch.

In the augmented data, all variables should be zero except the response variable which should have a value of 1, and the variables need to code for the particular random effects. In this example, ASReml augments the data but with values of zero for Angle. Therefore, a transformation Angle !REP 0 1 is specified to convert angles of zero to one.

If the data is manually augmented, it must also include the factor HSection and it is also necesary to adjust the coding of labelled ( !I or !A) factors. In this example, the actual temperatures are 175, 185, 195, 205, 215, 225. A zero in the augmented data would be treated as a seventh temperature which is not what we want. So, instead of writing Temperature !I, write Temperature !I 6 !REP 7 0 !PRUNE to force the seventh class to zero.

!HGLM HSection Repl Repl.Batch
declares we want to fit a HGLM with random terms Repl and Repl.Batch. If HSection has not been declared as a factor and the data have not been augmented, ASReml will create HSection and augment the data for random terms Repl and Repl.Batch

Angle !GAMMA !INV
declares we are to analyse the variate Angle as a Gamma variable with inverse link function

dev(Angle) is a label (variate name) for the dispersion.

!DISP 1
sets the dispersion scaling parameter to 1 as the dispersion is actually included in the model through the first level of at(Trait,2).Hsection

at(Trait,1).(at(HSection,1) Recipe Temp Recipe.Temp),
gives the model for the fixed terms of the mean model. at(Trait,1) indicates they are fitted just for the mean model. The augmented data is set up so that these factors are missing in the augmented data. Note that at(HSection,1) is equivalent to mu.hs(HS,1) and fits the intercept for section 1 only since the intercept applies just to the data records.

at(Trait,1).(Repl.hs(HS,2) Repl.Batch.hs(HS,3))
These are the two random terms fitted in the mean model. The augmented data is set up so that Repl is defined for Section 2, Repl and Batch are defined in section 3. The covariate hs(HS,2) prevents the design matrix for Repl being formed for Section 3.

at(Trait,2).HSection
declares the dispersion model. In this case it is a constant for each section, the constant being the log of the section variance.

residual units.id(2)
formally gives the residual variance structure but this can be omitted.

The hs(HS,2) model term ensures that the Repl term is fitted only for sections 1 and 2.

When creating the design matrix, the random term Repl must be created only for sections 1 and 2, not in section 3 which pertains to the Repl.Batch term. To form the Repl.Batch term, the augmented data for section 3 has codes for replicate and batch and so the design matrix for repl will also be formed for section 3. We do not want that, so we use hs(HS,2) which is based on the factor HSection and acts like a mask to keep the repl design only for sections 1 and 2. It creates a variable which is 1 for sections 1 (always) and 2 (nominated) which when interacted with repl, wipes out the design for repl with respect to section 3.

Interpreting the output

The primary output from the simple syntax is

   64 LogL= 899.034     S2=  1.0000        579 df  Dev/DF=   1.000

           - - - Results from analysis of Angle - - -
  Notice: While convergence of the LogL value indicates that the model
          has stabilized, its value CANNOT be used to formally test
          differences between Generalized Linear (Mixed) Models.

  Model_Term                        Sigma         Sigma   Sigma/SE   % C
  Residual           SCA_V  660   1.00000       1.00000       0.00   0 F

                                 Wald F statistics
      Source of Variation        NumDF     DenDF    F-inc         P-inc
   20 mu_Mean                        1     579.0 13016.00         <.001
   21 Recipe_Mean                    2     579.0     1.57         0.209
   22 Temperature_Mean               5     579.0    23.03         <.001
   23 Recipe.Temp_Mean              10     579.0     1.03         0.414
   24 Replicate_HS-2_Mean           15     579.0 41217.50         <.001
   25 Rep.Batch_HS-3_Mean           45     579.0 0.22E+06         <.001
   14 HSection_Dispn                 3     579.0  1235.09         <.001

Our primary interest is in the effects for HSection_Dispn which are

   Model_Term           Level           Effect    seEffect
   HSection_Dispn        1.002          -3.920      0.9483E-01
   HSection_Dispn        2.002          -10.63      0.3982
   HSection_Dispn        3.002          -12.29      0.3429

The 3 dispersion values for Residual, Repl and Repl.Batch are the antilog of these three effects (0.01984, 0.00002418 and 0.000004597).

Variance modelling

More sophisticated modelling of dispersion is achieved by writing an appropriate model which may require use of the MBF machinery. In this example, we model Blood Pressure data where the variance increases with Age. This leads to a DHGLM model. The ASReml code is

  !WORKSPACE 100 !RENAME !ARGS 3 !out
 Respiratory disorder data (Section 7.3.3, pages 221-224) !DOPART $1
 #patient,treatment,sex,age,center,baseline,past,y
 #1,placebo,male,46,1,0,0,0

  patient  *       # 1
  treatment  !A    # placebo
  sex  !A          # male
  age              # 46
  center  *        # 1
  baseline  *      # 0
  past             # 0
  y                # 0
  TOT !=1          # Formally defines NBIN else discards extra records

 Respiratory.csv  !SKIP 1

 !PART 1
 y !BIN !LOGIT ~ mu treatment center sex age baseline past,
       !r   patient.us(pol(age))
 residual units

 !PART 2   3
 TAB age ~ patient      #  Used to create MBF file
 !HGLM HSection patient

 !MBF  mbf(patient,1) PatAge.txt  !KEY 1 !REN PatAge

 !PART 2  # with dispersion
 y !BIN !LOGIT !TOT TOT dev(y) ~ +
       at(Tr,1).(at(HS,1) treatment center sex age baseline past patient) ,
       at(Tr,2).(at(HS,1)  at(HS,2)  at(HS,2).PatAge )
 residual units.id(2)

 !PART 3  # without dispersion
 y !BIN !LOGIT !TOT TOT dev(y) ~ +
       at(Tr,1).(at(HS,1) treatment center sex age baseline past patient) ,
       at(Tr,2).(-at(HS,1)  at(HS,2)  at(HS,2).PatAge )
 residual units.id(2)

Note that the distributions of the random effects are not formally stated but here default to Normal. They need not be in general.

Double HGLM

A common form of Double HGLM assumes the random effects are normal, in which case we actually fit a double GLMM, and do not require augmented data. But we do model the dispersion and can in fact correlate dispersion effects with mean effects. A popular example is

 Double HGLM algorithm by Majbritt Felleki and Lars Ronnegard
   Yval Xterm * Zterm *
 Hsection * !=1

 MF13.dat #dhglm.dat

 !HGLM Hsection    !MAXIT 100  !SIGMAPAR  !ASUV

 Yval  dev(Yval) ~ Trait Trait.X !R us(Trait !init 1 -.5 .8 !GH).Z  !H
 !PART 1// residual units.id(Trait)
 !PART 2// residual units.diag(Trait 1 1 !GFP)

The two response variables analysed in this analysis are the variable Yval and the deviance (residual squared) from the first trait. The first trait is weighted using weights predicted from the second (Trait_2) model. The second trait is weighted as a gamma variable using weights based on the same predicted values.

Note that we have initialized the Hsection variable at 1. and referred to it in the !HGLM Hsection qualifier. The !ASUV qualifier is required because the residual structure is uncorrelated (not units.us(Trait)). The !SIGMAPAR qualifier is required as the variance for the response is incorparated into the weights.

The 'gamma constraint' parameter H causes the variance components to be held at their initial values for the first 3 iterations (while the residual variance structure is estimated). The number of iterations the parameters are held at the initial values can be set with the dataline !FREEGH i qualifier.

Return to index