glmer
Someone asked whether it was possible to fit a mixed model in lme4
with box constraints on the fixed-effect parameters. It is, although (1) it requires a little bit of extra hacking (see below) and (2) it works most easily for generalized (non-Gaussian) linear MMs rather than LMMs (but see the second section).
library("lme4")
library("ggplot2"); theme_set(theme_bw())
Fit the basic (unconstrained) model:
gm0 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
data = cbpp, family = binomial)
(beta0 <- fixef(gm0))
## (Intercept) period2 period3 period4
## -1.398343 -0.991925 -1.128216 -1.579745
Start following the sequence laid out in ?modular
:
glmod <- glFormula(cbind(incidence, size - incidence) ~ period + (1 | herd),
data = cbpp, family = binomial)
devfun <- do.call(mkGlmerDevfun, glmod)
glmer
uses a two-stage optimization process. During the first stage (nAGQ=0
) we can’t constrain the fixed effect parameters (they’re profiled out of the optimization), so we just run optimizeGlmer
in the usual way from ?modular
:
opt1 <- optimizeGlmer(devfun)
Now set up the next optimization:
theta.lwr <- environment(devfun)$lower ## this changes after update!
devfun <- updateGlmerDevfun(devfun, glmod$reTrms)
rho <- environment(devfun)
nbeta <- ncol(rho$pp$X)
theta <- rho$pp$theta
Now do the second-stage optimization, by hand rather than using optimizeGlmer(...,stage=2)
as in ?modular
. The code below is actually a little bit simpler than the guts of optimizeGlmer
, because it’s not handling as wide a range of cases (e.g. user-specified starting values, additional control parameters, …). Instead of the usual constraints (theta parameters bounded below by rho$lower
, which corresponds to 0 for diagonal and -Inf
for off-diagonal elements, fixed-effect parameters unbounded), we constrain the fixed-effects parameters other than the intercept to [-1,1].
opt2 <- nloptwrap(par=c(rho$pp$theta,rep(0,nbeta)),
fn=devfun,
lower=c(theta.lwr,-Inf,rep(-1,nbeta-1)),
upper=c(rep(Inf,length(theta)),Inf,rep(1,nbeta-1)))
opt2$par
## [1] 0.6682799 -1.5037866 -0.8798190 -1.0000000 -1.0000000
opt2$fval
## [1] 186.114
## compare with 2*log(L) of original model
-2*c(logLik(gm0))
## [1] 184.0531
Now we have the estimated parameters and log-likelihood, but it’s a little bit prettier if we put the results into a merMod
object. (The only place we take a shortcut below is in the mc
argument, for which we’ll just substitute the matched call from the original model. This will only be a problem if we try to use update()
on the fitted model …)
(gm1 <- mkMerMod(rho, opt2, glmod$reTrms, glmod$fr, mc=gm0@call))
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(incidence, size - incidence) ~ period + (1 | herd)
## Data: cbpp
## AIC BIC logLik deviance df.resid
## 196.1140 206.2408 -93.0570 186.1140 51
## Random effects:
## Groups Name Std.Dev.
## herd (Intercept) 0.6683
## Number of obs: 56, groups: herd, 15
## Fixed Effects:
## (Intercept) period2 period3 period4
## -1.5038 -0.8798 -1.0000 -1.0000
Alternatively, it is (at least in principle) possible to skip the nAGQ=0
stage entirely and just fit the final model; this might be a good idea if the reason for constraining the fixed-effect parameters was to keep them away from a bad region.
The stuff below here is not really working properly. You probably shouldn’t even bother reading on unless you’re interested in the guts of lme4.
This is harder for LMMs because by default lmer
profiles out the fixed effect parameters, so we have to work a little bit harder to run a GLMM with a Gaussian response. All of the latter part of the machinery stays the same, so we’ll package it into a function:
glmmConstr <- function(devfun,mod,mc,beta.lwr,beta.upr,debug=FALSE) {
opt0 <- minqa::bobyqa(fn=devfun,par=1,lower=0,upper=Inf)
if (debug) cat("initial theta: ",opt0$par,"\n")
opt1 <- optimizeGlmer(devfun)
if (debug) cat("next theta: ",opt1$par,"\n")
rho <- environment(devfun)
nbeta <- ncol(rho$pp$X)
theta <- rho$pp$theta
theta.lwr <- rho$lower ## this changes after update!
devfun <- updateGlmerDevfun(devfun, mod$reTrms)
opt2 <- nloptwrap(par=c(theta,rep(0,nbeta)),
fn=devfun,
lower=c(theta.lwr,beta.lwr),
upper=c(rep(Inf,length(theta.lwr)),beta.upr))
if (debug) cat("final theta: ",opt2$par,"\n")
mkMerMod(rho, opt2, mod$reTrms, mod$fr, mc=mc)
}
Note that we also have to be careful because we’re working with reference class objects, so it’s quite possible to mess up an object by modifying it or a copy of it, even within a function …
Fit a basic lmer
example and set up a deviance function (slightly trickier)
fm0 <- lmer(Reaction~Days+(1|Subject),sleepstudy,REML=FALSE)
fmod <- glFormula(Reaction~Days+(1|Subject),
data=sleepstudy,family=gaussian,REML=FALSE)
## note we need family=gaussian() here -- an actual family object,
## not a string ("gaussian") or a family function (gaussian)
ldevfun <- do.call(mkGlmerDevfun, c(fmod, list(family=gaussian())))
Fit the unconstrained model:
fm1 <- glmmConstr(ldevfun,fmod,getCall(fm0),rep(-Inf,2),rep(Inf,2),
debug=TRUE)
## initial theta: 37.30982
## next theta: 37.30983
## final theta: 37.2637 250.2542 10.50123
fm0
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: Reaction ~ Days + (1 | Subject)
## Data: sleepstudy
## AIC BIC logLik deviance df.resid
## 1802.0786 1814.8505 -897.0393 1794.0786 176
## Random effects:
## Groups Name Std.Dev.
## Subject (Intercept) 36.01
## Residual 30.90
## Number of obs: 180, groups: Subject, 18
## Fixed Effects:
## (Intercept) Days
## 251.41 10.47
fm1
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: gaussian ( identity )
## Formula: Reaction ~ Days + (1 | Subject)
## Data: sleepstudy
## AIC BIC logLik deviance df.resid
## 1926.6293 1939.4011 -959.3147 1918.6293 176
## Random effects:
## Groups Name Std.Dev.
## Subject (Intercept) 1092.27
## Residual 29.31
## Number of obs: 180, groups: Subject, 18
## Fixed Effects:
## (Intercept) Days
## 250.3 10.5
Estimated fixed effects, std. dev. from lmer
:
c(fixef(fm0),attr(VarCorr(fm0)[[1]],"stddev"))
## (Intercept) Days (Intercept)
## 251.40510 10.46729 36.01208
from hacked glmer
:
c(fixef(fm1),getME(fm1,"theta"))
## (Intercept) Days Subject.(Intercept)
## 250.25415 10.50123 37.26370
ldevfun <- do.call(mkGlmerDevfun, c(fmod, list(family=gaussian())))
tvec <- seq(0,2,length.out=51)
lvec <- sapply(tvec,ldevfun)
qplot(tvec,lvec,geom=c("line","point"))
Hmmm. This seems upside-down, and has its peak in the wrong place (max \(\approx 0.3\) rather than getME(fm1,"theta")
?)