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
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??