modular fitting

library(lme4)
library(glmmTMB)
library(ggplot2); theme_set(theme_bw())
library(parameters)
library(mgcv)

Goals

• allow more flexible constructions, work around package limitations
• e.g. work-in-progress phyloglmm, multimembership models, examples from lmer vignette appendix

General steps

1. construct X, Z, Sigma in appropriate form
2. construct objective function (negative log-likelihood/deviance)
3. optimize objective function
4. package results

Modifications: (a) alter components between steps 1 and 2; (b) write a wrapper function for the objective function and optimize it instead

• glmer uses a two-stage optimization (first optimizing over \(\theta\) only, then \(\{\theta, \beta\}\))

lme4 bare-bones modular example

lmod <- lFormula(Reaction ~ Days + (Days|Subject), sleepstudy)
names(lmod)

## [1] "fr"      "X"       "reTrms"  "REML"    "formula" "wmsgs"

devfun <- do.call(mkLmerDevfun, lmod)
opt <- optimizeLmer(devfun)
fit <- mkMerMod(environment(devfun), opt, lmod$reTrms, fr = lmod$fr)

example: AR1 correlation structure

ar1cor <- function(rho, p) {
    R <- matrix(0, p, p)
    R[1,] = rho^(0:(p-1))        ## formula for 1st row
    c <- sqrt(1 - rho^2)         ## scaling factor: c^2 + rho^2 = 1
    R2 <- c * R[1,]              ## formula for 2nd row
    for (j in  2:p) {            ## shift elements in 2nd row for remaining rows
        R[j, j:p] = R2[1:(p-j+1)] 
    }
    return(R)
}
ar1cor(0.8, 5)

##      [,1] [,2] [,3]  [,4]   [,5]
## [1,]    1  0.8 0.64 0.512 0.4096
## [2,]    0  0.6 0.48 0.384 0.3072
## [3,]    0  0.0 0.60 0.480 0.3840
## [4,]    0  0.0 0.00 0.600 0.4800
## [5,]    0  0.0 0.00 0.000 0.6000

devfun2 <- function(sigma, rho) {
    cholfac <- sigma*t(ar1cor(rho, p))
    theta <- cholfac[lower.tri(cholfac, diag = TRUE)]
    devfun(theta)
}
d <- (read.csv("data/Cod_daily_depth_data.csv")
    |> na.omit()
    |> transform(ctemp = Temperature_1m-mean(Temperature_1m),
                 fdate = factor(as.Date(date)))
)

m1 <- glmmTMB(log(depth_mean_day) ~ 1 + ctemp +  (1 + ctemp | fish) +
            ar1(0 + fdate | fish),
        data = d)

## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; false convergence (8). See vignette('troubleshooting'),
## help('diagnose')

m2 <- update(m1, dispformula = ~0)
AIC(m1, m2)

##    df      AIC
## m1  8 1548.772
## m2  7 1546.772

ggplot(d, aes(ctemp, -log(depth_mean_day))) + 
  geom_point(alpha = 0.3) +
  facet_wrap(~fish)

m1 <- lmer(log(depth_mean_day) ~ 1 + ctemp +  (1 + ctemp | fish), data = d) # exercise
## library(mgcv)
m2 <- gam(log(depth_mean_day) ~ s(ctemp, factor(fish), bs = "fs"), data = d)
   
## performance::check_model(m1)
## performance::check_model(m2)
## nloptwrap

glmmTMB bare-bones modular example

library(glmmTMB)
## construct components
m1 <- glmmTMB(count ~ mined + (1|site),
              family=poisson, data=Salamanders, doFit = FALSE)
## make TMB object (don't optimize)
m2 <- fitTMB(m1, doOptim = FALSE)
## optimize
m3 <- with(m2, nlminb(par, objective = fn, gr = gr))
## construct fitted model object
m4 <- finalizeTMB(m1, m2, m3)

glmmTMB, rebuilding TMB object

If we modify the components of m1$env$data at this point … rebuild TMB structure (may be necessary)

m2 <- with(m2$env,
           TMB::MakeADFun(data,
                          parameters,
                          map = map,
                          random = random,
                          silent = silent,
                          DLL = "glmmTMB"))

use of the map argument in glmmTMB

• set parameters equal to starting values or fix values equal to each other
• e.g. for an AR1 model, map = list(theta = c(1, NA)), start = list(theta = 0, 1) would estimate the first parameter (log-SD) and fix the correlation to its starting value
• could pass a known correlation matrix (e.g. pedigree) by transforming into the glmmTMB parameterization and using map