# 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)

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)

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


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.