“Downdated VtV” example

This example represents a particular failure mode

This data set is for lme4 testing purposes only …

load("propData.RData")

The example works with lme4.0:

library("lme4.0")
g0 <- glmer(prophi ~ co2 * sex + (1 | plate), data = propData, weights = n, 
    family = binomial)
(f0 <- fixef(g0))
## (Intercept)        co2d        co2p        sexm   co2d:sexm   co2p:sexm 
##    -2.59536     5.35975     4.38614     0.07091     0.11542    16.41142
res <- c(getME(g0, "theta"), getME(g0, "beta"))
detach("package:lme4.0")

It works with lme4 if nAGQ=0: the results are not too far off (but we would like to do better)

library("lme4")
packageVersion("lme4")
## [1] '0.99999911.2'
g1 <- glmer(prophi ~ co2 * sex + (1 | plate), data = propData, weights = n, 
    family = binomial, nAGQ = 0)
rbind(f0, f1 <- fixef(g1))
##    (Intercept)  co2d  co2p    sexm co2d:sexm co2p:sexm
## f0      -2.595 5.360 4.386 0.07091   0.11542     16.41
##         -2.531 5.231 4.259 0.07171   0.09158     17.79

A naive nAGQ=1 fit fails:

glmer(prophi ~ co2 * sex + (1 | plate), data = propData, weights = n, family = binomial, 
    verbose = 10)
## (NM) 1: f = inf at  0.795257  -2.53102   5.23126   4.25852 0.0717101  0.091576    17.791
## (NM) init_pos <= d_n
## (NM) 2: f = 49.9185 at  0.795257  -2.53102   5.23126   4.25852 0.0717101  0.091576    17.791
## (NM) init_pos <= d_n
## (NM) 3: f = 49.9185 at  0.795257  -2.53102   5.23126   4.25852 0.0717101  0.091576    17.791
## (NM) init_pos <= d_n
## (NM) 4: f = 49.9185 at  0.795257  -2.53102   5.23126   4.25852 0.0717101  0.091576    17.791
## (NM) init_pos <= d_n
## (NM) 5: f = 49.8816 at  0.795257  -2.53102   5.39571   4.25852 0.0717101  0.091576    17.791
## (NM) init_pos <= d_n
## (NM) 6: f = 49.8816 at  0.795257  -2.53102   5.39571   4.25852 0.0717101  0.091576    17.791
## (NM) init_pos <= d_n
## (NM) 7: f = 49.8816 at  0.795257  -2.53102   5.39571   4.25852 0.0717101  0.091576    17.791
## (NM) init_pos <= d_n
## Error: Downdated VtV is not positive definite
vv <- unlist(getME(g1, c("theta", "beta")))
ff2 <- function(p, reset = TRUE, debug = FALSE) {
    if (reset) 
        assign("ff", glmer(prophi ~ co2 * sex + (1 | plate), data = propData, 
            weights = n, family = binomial, devFunOnly = TRUE), .GlobalEnv)

    r <- ff(p)
    if (debug) 
        cat(p, r, "\n")
    r
}
opt1 <- optim(fn = ff2, par = vv)
cbind(res, opt1$par)
##                        res         
## plate.(Intercept)  0.80183  0.80047
##                   -2.59536 -2.59564
##                    5.35975  5.35760
##                    4.38614  4.38453
##                    0.07091  0.07116
##                    0.11542  0.11699
##                   16.41142 18.65621

Results are identical except for the last parameter, which is essentially an approximation to infinity (the difference between the parameters corresponds to a difference between 1-7.5435 × 10-8 and 1-7.563 × 10-9.

Trying to understand what's going on with changes in the environment of the deviance function. Now I can't get it to break …

ff <- update(g1, nAGQ = 1, devFunOnly = TRUE)
fun <- function() environment(ff)$pp$delu
res2 <- c(0.795257, -2.53102, 5.39571, 4.25852, 0.0717101, 0.091576, 17.791)
fun()
## [1] -1.0991  0.4341  1.3864 -0.2498 -0.4715
(deluCopy <- environment(ff)$pp$copy(shallow = FALSE)$delu)
## [1] 0 0 0 0 0

Copying zeros out the value of delu… ??? Or copy the naive way …

(deluCopy2 <- environment(ff)$pp$delu)
## [1] -1.0991  0.4341  1.3864 -0.2498 -0.4715
ff(res2)
## [1] 49.88
fun()
## [1] -1.1583  0.4034  1.3835 -0.2864 -0.4974
deluCopy  ## still zero
## [1] 0 0 0 0 0
deluCopy2  ## has changed in concert
## [1] -1.1583  0.4034  1.3835 -0.2864 -0.4974

Further calls don't change the result any more:

ff(res2)
## [1] 49.88
fun()
## [1] -1.1583  0.4034  1.3835 -0.2864 -0.4974
ff(res2)
## [1] 49.88
fun()
## [1] -1.1583  0.4034  1.3835 -0.2864 -0.4974

Alternatively, rather than completely resetting the deviance function at each step, we can just zero out delu at each step – this is much quicker!

ff <- update(g1, nAGQ = 1, devFunOnly = TRUE)
ff3 <- function(p, reset = TRUE, debug = FALSE) {
    if (reset) 
        environment(ff)$pp$setDelu(rep(0, 5))
    r <- ff(p)
    if (debug) 
        cat(p, r, "\n")
    r
}
opt2 <- optim(fn = ff3, par = vv)

Compare all three approaches:

cbind(lme4.0 = res, opt1 = opt1$par, opt2 = opt2$par)
##                     lme4.0     opt1     opt2
## plate.(Intercept)  0.80183  0.80047  0.80167
##                   -2.59536 -2.59564 -2.59536
##                    5.35975  5.35760  5.35999
##                    4.38614  4.38453  4.38657
##                    0.07091  0.07116  0.07075
##                    0.11542  0.11699  0.11543
##                   16.41142 18.65621 17.98401

Part II

Example from issues list: can we fix it by resetting delu?

gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd), data = cbpp, 
    family = binomial)
options(digits = 7)
deviance(gm1)
## [1] 184.0527

Extract deviance function and parameter values:

ff <- update(gm1, devFunOnly = TRUE)
tt <- c(getME(gm1, "theta"), getME(gm1, "beta"))

A brief look at the internals:

environment(ff)$pp$delu
##  [1]  0.86771368 -0.51335392  0.57488085  0.01855007 -0.34804988
##  [6] -0.67178539  1.33192812  0.88905042 -0.40504049 -0.89134801
## [11] -0.18867379 -0.13868771 -1.12183911  1.46727803 -0.87062285
environment(ff)$pp$delb
## [1] -1.3604757 -0.9761732 -1.1110728 -1.5596757

Running the deviance function gives a slightly different answer from the deviance recorded in the model object, and modifies the internal copy of delb (but not delu, in this case):

ff(tt)  ## 184.0531
## [1] 184.0531
environment(ff)$pp$delu
##  [1]  0.91870073 -0.46537131  0.63257996  0.06117162 -0.29583923
##  [6] -0.62320984  1.38484199  0.93326691 -0.37001415 -0.84223209
## [11] -0.13176164 -0.10090476 -1.07421014  1.51145455 -0.82594416
environment(ff)$pp$delb  ## has been zeroed out
## [1] 0 0 0 0

Calling the function again gives the recorded value, and doesn't further change delu and delb … further calls continue to give the same value.

ff(tt)  ## 184.0527
## [1] 184.0527
environment(ff)$pp$delu
##  [1]  0.91870072 -0.46537131  0.63257996  0.06117162 -0.29583924
##  [6] -0.62320985  1.38484199  0.93326691 -0.37001415 -0.84223209
## [11] -0.13176164 -0.10090477 -1.07421015  1.51145455 -0.82594416
environment(ff)$pp$delb  ## has been zeroed out
## [1] 0 0 0 0
ff(tt)
## [1] 184.0527

But: re-computing the deviance function and then forcing delb to a vector of zeros doesn't give us back the original value:

ff <- update(gm1, devFunOnly = TRUE)
environment(ff)$pp$setDelb(rep(0, 4))
## NULL
ff(tt)
## [1] 184.0531

What is the most convenient way to make a full, deep copy of an environment (including reference class objects contained therein) and compare it with (the contents of) another environment??