# Why does PQL vs MCMLEs optimisation give wildly different variance estimates?

In order to better understand modelling GLMM's using R I decided to re-do the example given in this Introduction to GLMM Package using the salamander dataset (provided in the glmm package).

The package has two optimisation methods: PQL (default) and MCMLE. The two methods give nearly identical fixed effect estimates, but the random effect (variance) estimates are wildly different. Here is the output of the two (all equal except for "doPQL=TRUE/FALSE").

doPQL = TRUE

> summary(mod)

Call:
glmm(fixed = Mate ~ 0 + Cross, random = list(~0 + Female, ~0 +
Male), varcomps.names = c("F", "M"), data = dat2, family.glmm = bernoulli.glmm,
m = 10000, debug = TRUE)

Link is: "logit (log odds)"

Fixed Effects:
Estimate Std. Error z value Pr(>|z|)
CrossR/R   1.0481     0.4454   2.353   0.0186 *
CrossR/W   0.3687     0.4645   0.794   0.4273
CrossW/R  -1.8569     0.4593  -4.043 5.29e-05 ***
CrossW/W   1.0000     0.4723   2.118   0.0342 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Variance Components for Random Effects (P-values are one-tailed):
Estimate Std. Error z value Pr(>|z|)/2
F   1.1177     0.3397   3.291     0.0005 ***
M   1.2190     0.5537   2.202     0.0138 *


doPQL = FALSE

> summary(mod)

Call:
glmm(fixed = Mate ~ 0 + Cross, random = list(~0 + Female, ~0 +
Male), varcomps.names = c("F", "M"), data = dat2, family.glmm = bernoulli.glmm,
m = 10000, doPQL = FALSE, debug = TRUE)

Link is: "logit (log odds)"

Fixed Effects:
Estimate Std. Error z value Pr(>|z|)
CrossR/R   0.9593     0.2396   4.004 6.24e-05 ***
CrossR/W   0.3180     0.2291   1.388 0.165139
CrossW/R  -1.3653     0.2713  -5.033 4.84e-07 ***
CrossW/W   0.8287     0.2404   3.447 0.000567 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Variance Components for Random Effects (P-values are one-tailed):
Estimate Std. Error z value Pr(>|z|)/2
F  0.36805    0.06722   5.475   2.19e-08 ***
M  0.44576    0.08144   5.473   2.21e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


What's going on here? Why do the two techniques give such wildly different estimates, and which one should be used?

The dataset in this example has binary response variable and is of this form:

> salamander[1:10,]
Mate Cross Female Male
1     1   R/R     10   10
2     1   R/R     11   14
3     1   R/R     12   11
4     1   R/R     13   13
5     1   R/R     14   12
6     1   R/W     15   28
7     0   R/W     16   27
8     0   R/W     17   25
9     0   R/W     18   29
10    0   R/W     19   26


## 1 Answer

When you want to fit a Generalized Linear Mixed Model under maximum likelihood, you need to numerically calculate the integral over the random effects in the definition of the observed data log-likelihood.

The glmm package does this integration using Monte Carlo and the importance sampling algorithm. The doPQL controls the importance sampling distribution for the random effects. I.e., if it is TRUE, then this distribution is centered around the PQL estimates of the random effects (and with the appropriate hessian matrix of the penalized-quasi likelihood), whereas if it is FALSE it is centered around 0 and with the identity matrix as the covariance matrix of the random effects. Hence, the optimization algorithm remains the same, the only thing that changes is the importance sampling distribution for the Monte Carlo procedure.

As you have seen, this can have an impact on the resulting estimates. The reason why this occurs is that for some sample units the conditional distribution of the random effects given the observed data is centered far away from 0 and with a variance different than 1. In general, the adaptive Gaussian quadrature is considered to be the most accurate approach to approximate these integrals over the random effects.