example of profiling a gamma model in lme4

lme4_gamma_profile

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