Generalized Linear Mixed Models
Introduction
ASReml uses an approximate likelihood technique called
penalized quasi-likelihood (PQL)
to analyse
data sampled from one of the common members of the exponential family.
In this section we present a few examples to demonstrate the coding in ASReml.
Binomial analysis of Footrot score
Mohommad Alwan (pers comm) for his Master thesis at Massey University
scored the feet of 2513 lambs born in 1980 and 1981.
The lambs were from 5 mating groups: 7 Perendale rams over Perendale ewes in 1980, 6 Booroola by Romney rams over Perendale ewes in 1980,
3 Booroola rams over Romney ewes in 1980,
6 Perendale rams over Perendale ewes in 1981, and
12 Booroola by Romney rams (from froup 3) over Perendale ewes in 1981.
This data was analysed by Gilmour (1984) and Gilmour et al. (1987).
The data file LAMB.DAT contains grouped data for the 68 combinations
of Sex and Sire for two footshape classes: FS1, all four feet are normal,
FS2, one foot is deformed; and two indicator variables for the presence
of disease conditions Scald and Rot. No scald or rot was present in group 4 lambs and these responses have been set to missing.
The genetic relationships among sires are ignored in this analysis
although it would just require a sire relationship matrix to include them.
Our first analysis is of the incidence of foot rot on the Normal scale
as a weighted analysis to mimic analysis of the ungrouped data.
Using 56 of the 68 records (ignoring Group 4), there
are 1960 (=56 cross 35.00) observations and so we use the !DF 1904 (=1960-56) qualifier
to get the correct residual degrees of freedom for this analysis of the
proportion with footrot. The !YSS 62.54249 qualifier adds 62.54249=67-4.45751
to the Total Sum of Squares so that it includes the extra variation associated
with the extra degrees of freedom. There were 67 (=56*1.196) cases of foot
rot so the Total uncorrected Sum of Squares for a binary variable should be 67.
However the weighted sum of squares for the pRot values is only
4.45751 (for example the first record contributes 1/39=(1/39)2 cross 39 instead of 1.0. 4.45751
was discovered from the .asl file on the line
4.45751 SSPD before inserting the !YSS qualifier.
The transformations in the code which follows convert
Scald and Rot to 'missing' for group 4.
Lamb data from ARG thesis page 177-8
Year GRP 5 !V99=V2 !==4 !M1
SEX SIRE !I
Total
FS1 FS2 Scald !+V99 Rot !+V99
pRot !=Rot !/Total
# 1 1 1 101 39 33 6 6 1
LAMB.DAT !skip 1
!DF 1904 !YSS 62.54249
pRot !TOTAL=Total ~ mu SEX GRP !r SIRE
predict SEX 0 1 GRP 1 2 3 5
The pertinant results are
Univariate analysis of pRot
Summary of 56 records retained of 68 read
Model term Size #miss #zero MinNon0 Mean MaxNon0 StndDevn
1 Year 0 0 1.000 1.536 2.000 0.5032
2 GRP 5 0 0 1 3.1429 5
3 SEX 0 28 1.000 0.5000 1.000 0.5045
4 SIRE 34 0 0 1 17.0714 34
5 Total Weight 0 0 16.00 35.00 64.00 12.89
6 FS1 0 0 6.000 23.46 50.00 10.76
7 FS2 0 0 3.000 10.14 30.00 5.661
8 Scald 0 13 1.000 3.071 16.00 3.458
9 Rot 0 19 1.000 1.196 4.000 1.151
10 pRot Variate 0 19 0.1754E-01 0.3606E-01 0.1818 0.3833E-01
11 mu 1
12 SEX.GRP 5 3 SEX : 1 2 GRP : 5
Forming 46 equations: 12 dense.
Initial updates will be shrunk by factor 0.224
Notice: Algebraic Denominator DF calculation is not available
Numerical derivatives will be used.
Notice: 4 singularities detected in design matrix.
1 LogL= 2423.41 S2= 0.32397E-01 1952 df : 1 components restrained
2 LogL= 2431.71 S2= 0.32792E-01 1952 df 0.6325E-02 1.000
3 LogL= 2431.80 S2= 0.32737E-01 1952 df 0.9265E-02 1.000
4 LogL= 2431.80 S2= 0.32738E-01 1952 df 0.9200E-02 1.000
Final parameter values 0.92543E-02 1.0000
- - - Results from analysis of pRot - - -
Approximate stratum variance decomposition
Stratum Degrees-Freedom Variance Component Coefficients
SIRE 25.70 0.506971E-01 59.7 1.0
Residual Variance 15.83 0.327367E-01 0.0 1.0
Source Model terms Gamma Component Comp/SE % C
SIRE 34 34 0.918415E-02 0.300659E-03 0.98 -22 P
Variance 56 1952 1.00000 0.327367E-01 2.81 0 P
Wald F statistics
Source of Variation NumDF DenDF F-inc Prob
11 mu 1 19.9 42.79 <.001
3 SEX 1 16.2 0.02 0.882
2 GRP 3 21.9 2.04 0.139
12 SEX.GRP 3 16.1 0.39 0.763
Notice: The DenDF values are calculated ignoring fixed/boundary/singular
variance parameters using numerical derivatives.
4 SIRE 34 effects fitted ( 6 are zero)
Two things stand out in this analysis. From a genetic perspective, the
heritability estimate is 0.0364=4 *.0003007/(.0003007+.0327367)
This can be calculated
in ASReml with the .pin file commands
F GenVar 1*4
F TotVar 1 2
H heritability 3 4
Secondly, there is little evidence of significant difference between classes.
The predicted values are
| Sex | PxP 1980 | BRxP 1980 | BxR 1980 | BRxP 1981
|
| 0 | 0.0183 +/- 0.0130 |
0.0432 +/- 0.0126 |
0.0758 +/- 0.0268 |
0.0305 +/- 0.0111
|
| 1 | 0.0152 +/- 0.0132 |
0.0375 +/- 0.0124 |
0.0603 +/- 0.0244 |
0.0425 +/- 0.0108
|
An analysis of footrot as a binomial variable using the logistic link
is performed by the model line (and dropping the !DF qualifier).
Rot !bin !TOTAL=Total ~ mu SEX GRP SEX.GRP !r SIRE .16783
pertinant results are
Distribution and link: Binomial; Logit Mu=P=1/(1+exp(-XB))
V=Mu(1-Mu)/N
Warning: The LogL value is unsuitable for comparing GLM models
Notice: 4 singularities detected in design matrix.
1 LogL=-28.1544 S2= 1.0000 48 df Dev/DF= 0.9060
2 LogL=-28.7417 S2= 1.0000 48 df Dev/DF= 0.8897
3 LogL=-28.7186 S2= 1.0000 48 df Dev/DF= 0.8805
4 LogL=-28.6705 S2= 1.0000 48 df Dev/DF= 0.8551
5 LogL=-28.6494 S2= 1.0000 48 df Dev/DF= 0.8238
6 LogL=-28.6687 S2= 1.0000 48 df Dev/DF= 0.7959
7 LogL=-28.6774 S2= 1.0000 48 df Dev/DF= 0.7915
8 LogL=-28.6784 S2= 1.0000 48 df Dev/DF= 0.7909
9 LogL=-28.6785 S2= 1.0000 48 df Dev/DF= 0.7908
Final parameter values 0.26321 1.0000
Deviance from GLM fit 48 37.96
Variance heterogeneity factor [Deviance/DF] 0.79
- - - Results from analysis of Rot - - -
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.
Approximate stratum variance decomposition
Stratum Degrees-Freedom Variance Component Coefficients
SIRE 3.10 0.263207 1.0
Source Model terms Gamma Component Comp/SE % C
SIRE 34 34 0.263207 0.263207 1.25 0 P
Variance 56 48 1.00000 1.00000 0.00 0 F
Wald F statistics
Source of Variation NumDF DenDF F-inc Prob
11 mu 1 20.2 418.38 <.001
3 SEX 1 48.0 0.02 0.881
2 GRP 3 21.5 1.99 0.146
12 SEX.GRP 3 NA 0.36 NA
The effects in this analysis are on a logistic
scale with a variance of
3.28987=π2/3 and so the heritability on the underlying (logistic) scale
is 0.296=4 * 0.2632/( 3.28987+0.26321). This can be calculated
in ASReml with the .pin file commands
F GenVar 1*4
F TotVar 1 4*3.28987
H heritability 3 4
Repeating the analysis on the Probit scale by inserting !PROBIT after
!BIN in the model line produces a Sire component of 0.0514 on the
Probit scale which has an underlying variance of 1.0. The heritabily estimate is then 0.196.
Given the incidence (0.034), the heritability on the probit scale
is expected to be around 0.215=0.0364/(z2/pq) where z=0.0758 is the ordinate of a Normal(0,1) corresponding to p=1-q=0.034.
The preceding Wald F Statistics pertain to the working variable
created as part of the PQL analysis.
The SEX.GRP interaction is clearly not significant
even though ASReml was not able to calculate a plausible value
for the Denominator DF for this summarized data. The predicted
means shown below are not that different from those
obtained from analysis on the 0,1 scale but the standard errors are very different.
These predicted means have been backtransformed by \ASReml\ from the underlying (logistic)
scale to the probablity scale.
The initial analysis (on the 0,1 probability scale) ignores the variance differences
associated with binomial data.
| PxP 1980 | BRxP 1980 | BxR 1980 | BRxP 1981
| |
0 | 0.0180 +/- 0.0070 |
0.0430 +/- 0.0124 |
0.0748 +/- 0.0323 |
0.0281 +/- 0.0083
| |
1 | 0.0151 +/- 0.0063 |
0.0373 +/- 0.0110 |
0.0592 +/- 0.0257 |
0.0401 +/- 0.0103 | |
ASReml has an
'Analysis of Deviance' option which we now demonstrate.
In a mixed model, the variance components will change
depending on which fixed terms are in the model. This
will invalidate the Analysis of Deviance unless the variance components
are fixed at the full model solution.
So, fitting the model line
Rot !bin !TOT=Total !AODEV ~ mu SEX GRP SEX.GRP !r SIRE .2632 !GF}
produces the Analysis of Deviance
Analysis of Deviance Table for Rot
Source of Variation df Deviance Derived F
SEX 1 0.02 0.021
GRP 3 4.35 1.833
SEX.GRP 3 1.16 0.487
Deviance from GLM fit 48 37.96
Variance heterogeneity factor [Deviance/DF] 0.79
The Deviance is the deviance calculated from the
binomial part of the log-likelihood. This is distinct from the
log-likelihood obtained by the REML algorithm which pertains
to the working variable. Since the working variable changes
with the model fitted, the LogL values are not comparable between
models.
The heterogeneity factor is the Deviance / df
and gives some indication as to how well the discrete distribution
fits the data. A value greater than 1 suggests the data is over-dispersed,
that is the data values are more variable than expected
under the chosen distribution.
There is also a !DISPERSION [d] qualifier. If d is supplied, it
serves as a scaling factor for the weights in the analysis, changing the
reported variances and standard deviations. If d is not supplied,
it is estimated from the residual as the model is fitted to the working variable.
ASReml solves for the linear effects twice (see the !GLMM qualifier)
iteration of the variance components so that the variance component updates
are based on solutions obtained using the same variance parameters.
I.e. We start with a set of solutions and some parameters.
We use these to update the solutions. Then use the updated solutions
to update the variance parameter.
Bivariate analysis of binomial Foot score and Normal weight
The data file BINNOR.txt contains the expanded version (2513 records) of
the lamb data from the previous example augmented with an extra
simulated variable YVar.
It was
created from the summarized data without knowing which actual individuals
had which combinations of trait values.
The binary variable Score1 indicates
whether all four feet are sound. The following code produces a bivariate
analysis of Score1 on the underlying logistic scale and YVar
on the Normal scale.
Lamb data from ARG thesis page 177-8
Year GRP 5 !V99=V2 !==4 !M 1
SEX SIRE !I
Score1
Score2 Scald !+V99 Rot !+V99
YVar
binnor.txt !skip 1 !ASUV !MAXIT 40
Score1 YVar !bin ~ Trait.SEX Trait.GRP !r Trait.SIRE
1 2 1
2513
2 0 US !GFPP
1 .01 0.25
Trait.SIRE 2
Trait 0 US 0.015 0.01 1.05
SIRE
There are several issues addressed in this code.
!ASUV is required, and if there had been any missing values
in the data, the fixed model term mv would also be required.
ASReml constructs the R matrix by scaling the reported matrix by the
binomial variance calculated from the fitted value of the binomial variate.
Consequently, to avoid over/under dispersion being also fitted,
the residual 'variance' for the binomial trait is fixed at 1.0
by giving its initial value as 1.0 and using the qualifier !GFPP.
The response variables must be listed before the qualifiers.
If written as \as{Score !BIN YVar}, YVar would be parsed as an argument to !BIN
rather than as a response variable.
Only one categorical response is permitted, and it must be specified first.
output follows.
Distribution and link: Binomial; Logit Mu=P=1/(1+exp(-XB))
V=Mu(1-Mu)/N
Warning: The LogL value is unsuitable for comparing GLM models
1 LogL=-894.974 S2= 1.0000 5014 df Dev/DF= 0.6196
2 LogL=-894.554 S2= 1.0000 5014 df Dev/DF= 0.6194
3 LogL=-890.600 S2= 1.0000 5014 df Dev/DF= 0.6178
4 LogL=-884.431 S2= 1.0000 5014 df Dev/DF= 0.6144
5 LogL=-885.759 S2= 1.0000 5014 df Dev/DF= 0.6109
6 LogL=-892.413 S2= 1.0000 5014 df Dev/DF= 0.6085
7 LogL=-896.969 S2= 1.0000 5014 df Dev/DF= 0.6077
8 LogL=-897.941 S2= 1.0000 5014 df Dev/DF= 0.6076
9 LogL=-897.962 S2= 1.0000 5014 df : 1 components restrained
10 LogL=-897.962 S2= 1.0000 5014 df Dev/DF= 0.6076
11 LogL=-897.961 S2= 1.0000 5014 df Dev/DF= 0.6076
Deviance from GLM fit 5014 3046.50
Variance heterogeneity factor [Deviance/DF] 0.61
- - - Results from analysis of Score1 YVar - - -
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.
Source Model terms Gamma Component Comp/SE % C
Residual UnStructured 2 1 -0.162615E-03 -0.162615E-03 -0.03 0 P
Residual UnStructured 2 2 0.255609 0.255609 35.20 0 P
Trait.SIRE UnStructured 1 1 0.166092 0.166092 2.73 0 U
Trait.SIRE UnStructured 2 1 0.330313E-02 0.330313E-02 0.07 0 U
Trait.SIRE UnStructured 2 2 0.303900 0.303900 3.76 0 U
Covariance/Variance/Correlation Matrix UnStructured Residual
1.000 -0.3216E-03
-0.1626E-03 0.2556
Covariance/Variance/Correlation Matrix UnStructured Trait.SIRE
0.1661 0.1470E-01
0.3303E-02 0.3039
Wald F statistics
Source of Variation NumDF DenDF-con F-inc F-con M P-con
11 Trait.SEX 2 NA 393.15 76.10 A NA
12 Trait.GRP 10 40.9 1993.52 1993.52 A <.001
Notice: The DenDF values are calculated ignoring fixed/boundary/singular
variance parameters using numerical derivatives.
YVar data was artificially created and the SIRE variance is too large
represent purely genetic variance.
Ordinal analysis of Footrot score
By way of introduction to ordinal analysis in ASReml consider
the cheese data from page 175 of McCullagh and Nelder (1994).
Four cheeses were scored on a ninepoint scale by 52 tasters giving
Table 1 Response frequencies in a cheese tasting experiment
| I | II | III | IV | V | VI | VII | VIII | IX | Total
| | A | 0 | 0 | 1 | 7 | 8 | 8 | 19 | 8 | 1 | 52
| | B | 6 | 9 | 12 | 11 | 7 | 6 | 1 | 0 | 0 | 52
| | C | 1 | 1 | 6 | 8 | 23 | 7 | 5 | 1 | 0 | 52
| | D | 0 | 0 | 0 | 1 | 3 | 7 | 14 | 16 | 11 | 52
| |
There are several ways of supplying the data for multinomial analysis.
In this case, totals in the 9 classes are supplied in a single grouped
response.
It is analysed using a multiple (8) threshold model as in
McCullagh and Nelder (1994) with the \ASReml\ code
McCullagh and Nelder Cheese example p 175
Cheese !A
Rating !G 9 Total
Cheese.txt
Rating !MULT 9 !CUM ~ Trait Cheese
PREDICT Cheese
where Cheese.txt contains the data laid out as in Table 1.
The model term Trait fits the thresholds and interpreting the model
as a threshold model implies it should not be interacted with
other terms. Nevertheless, sometimes an interaction is fitted.
Note that \ASReml\ does not have a procedure for multinomial data
which is not ordered (except as fitted with a log linear model),
nd fitting a bivariate analysis involving a multinomial trait
is not possible.
The output is
Univariate analysis of Rating
Summary of 4 records retained of 4 read
Model term Size #miss #zero MinNon0 Mean MaxNon0 StndDevn
1 Cheese 4 0 0 1 2.5000 4
2 Rating Variate 0 2 1.000 1.750 6.000 2.872
2 Rating Variate 0 2 1.000 2.500 9.000 4.359
2 Rating Variate 0 1 1.000 4.750 12.00 5.500
2 Rating Variate 0 0 1.000 6.750 11.00 4.193
2 Rating Variate 0 0 3.000 10.25 23.00 8.770
2 Rating Variate 0 0 6.000 7.000 8.000 0.8165
2 Rating Variate 0 0 1.000 9.750 19.00 8.221
2 Rating Variate 0 1 1.000 6.250 16.00 7.411
2 Rating Variate 0 2 1.000 3.000 11.00 5.354
3 Total 0 0 52.00 52.00 52.00 0.000
4 Trait 8
Forming 12 equations: 12 dense.
Initial updates will be shrunk by factor 0.010
Distribution and link: Cum. Multinomial; Logit P=1/(1+exp(-XB))
Warning: The LogL value is unsuitable for comparing GLM models
Notice: 1 singularities detected in design matrix.
1 LogL=-26.4243 S2= 1.0000 21 df Dev/DF= 0.3356
2 LogL=-26.4503 S2= 1.0000 21 df Dev/DF= 0.3376
3 LogL=-26.4506 S2= 1.0000 21 df Dev/DF= 0.3376
4 LogL=-26.4506 S2= 1.0000 21 df Dev/DF= 0.3376
5 LogL=-26.4506 S2= 1.0000 21 df Dev/DF= 0.3376
Deviance from GLM fit 21 20.31
Variance heterogeneity factor [Deviance/DF] 0.97
- - - Results from analysis of Rating - - -
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.
Source Model terms Gamma Component Comp/SE % C
Wald F statistics
Source of Variation NumDF F-inc
4 Trait 8 17.45
1 Cheese 3 38.38
Warning: These Wald F statistics are based on the working variable and are
not equivalent to an Analysis of Deviance. Standard errors are scaled
by the variance of the working variable, not the residual deviance.
Finished: 17 Jun 2008 13:19:51.484 LogL Converged
Reverting to the collapsed lamb data, the two response variables FS1
and FS2 contain counts of the lambs with all feet sound, and with
one foot deformed, respectively. The count for those with two or more deformed
is given by difference from Total. A threshold model analysis
of this data is given by the model line
FS1 FS2 !mult 3 !TOTAL=Total ~ Trait SEX GRP !r SIRE
output
Notice: 1 singularities detected in design matrix.
1 LogL=-105.631 S2= 1.0000 129 df Dev/DF= 1.082
2 LogL=-105.632 S2= 1.0000 129 df Dev/DF= 1.082
3 LogL=-105.631 S2= 1.0000 129 df Dev/DF= 1.081
4 LogL=-105.628 S2= 1.0000 129 df Dev/DF= 1.080
5 LogL=-105.627 S2= 1.0000 129 df Dev/DF= 1.079
6 LogL=-105.627 S2= 1.0000 129 df Dev/DF= 1.078
Deviance from GLM fit 129 139.09
Variance heterogeneity factor [Deviance/DF] 1.08
- - - Results from analysis of FS1 FS2 - - -
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.
Source Model terms Gamma Component Comp/SE % C
SIRE 34 34 0.174697 0.174697 2.80 0 P
Wald F statistics
Source of Variation NumDF DenDF F-inc Prob
11 Trait 2 77.8 405.40 <.001
3 SEX 1 129.0 5.61 0.020
2 GRP 4 30.0 8.03 <.001
Notice: The DenDF values are calculated ignoring fixed/boundary/singular
variance parameters using numerical derivatives.
Warning: These Wald F statistics are based on the working variable and are
not equivalent to an Analysis of Deviance. Standard errors are scaled
by the variance of the working variable, not the residual deviance.
Solution Standard Error T-value T-prev
2 GRP
2 -0.727155 0.273336 -2.66
3 -1.76491 0.356573 -4.95 -2.93
4 -1.19399 0.273168 -4.37 1.61
5 -0.915605 0.242677 -3.77 1.16
3 SEX
1 -0.197719 0.856093E-01 -2.31
11 Trait
1 1.54993 0.200125 7.74
2 3.82051 0.216314 17.66 27.12
4 SIRE 34 effects fitted
Finished: 18 Jun 2008 12:35:09.062 LogL Converged
Return to start