Skip to content

Instantly share code, notes, and snippets.

@bbolker
Created August 21, 2018 16:53
Show Gist options
  • Save bbolker/43fca2ac7fcd52f1f4971224ae9dc096 to your computer and use it in GitHub Desktop.
Save bbolker/43fca2ac7fcd52f1f4971224ae9dc096 to your computer and use it in GitHub Desktop.
example of profiling a gamma model in lme4
library(lme4)
set.seed(101)
dd <- data.frame(x=rnorm(200),f=factor(rep(1:10,each=20)))
dd$y <- simulate(~x+(1|f),
newdata=dd,
newparams=list(beta=c(10,1),
theta=1,
sigma=5),
family=Gamma(link="identity"))[[1]]
## NB not sure what sigma really is here??? must be shape parameter
## (CV = (1/sqrt(shape)))?
plot(y~x,data=dd)
## fit
fit1 <- glmer(y~x+(1|f), data=dd, family=Gamma(link="identity"))
try(pp <- profile(fit1)) ## oops, not implemented
## Let's implement a crude version.
## Not sure why this wasn't done, except that profiling the dispersion
## (scale) parameter is non-trivial [requires re-parameterizing the
## function and figuring how to fix the scale parameter ...]
## if we have a simple model (only scalar random-effects terms)
## OR if we are only interested in profiling the fixed effects
## we don't have to worry about the reparameterization
## (because the 'theta' parameters are just standard deviations anyway
## in this case)
## extract deviance function
dfun <- update(fit1,devFunOnly=TRUE)
## best-fit params
p0 <- unlist(getME(fit1,c("theta","beta")))
d0 <- c(-2*logLik(fit1))
all.equal(dfun(p0),d0)
proffun <- function(p,par=1,target=d0+qchisq(0.95,1)) {
pvec <- p0
pvec[par] <- p
dfun(pvec)-target
}
proffun(p0[1])
## check
curve(Vectorize(proffun)(x),from=0,to=p0[1]+10)
##
## compute std devs to get an estimate of scale
sdvec <- sqrt(diag(solve(fit1@optinfo$derivs$Hessian)))
res <- matrix(NA,ncol=2,nrow=length(p0),
dimnames=list(c(names(getME(fit1,"theta")),
names(fixef(fit1))),
c("lwr","upr")))
for (i in seq(nrow(res))) {
lwr <- p0[i]-4*sdvec[i]
if (i<=length(fit1@lower)) {
lwr <- max(fit1@lower[i],lwr)
}
upr <- p0[i]+4*sdvec[i]
u_lwr <- try(uniroot(proffun,interval=c(lwr,p0[i]),par=i))
if (!inherits(u_lwr,"try-error")) {
res[i,"lwr"] <- u_lwr$root
}
u_upr <- try(uniroot(proffun,interval=c(p0[i],upr),par=i))
if (!inherits(u_upr,"try-error")) {
res[i,"upr"] <- u_upr$root
}
}
res
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment