glmmTMB: start a new optimization at the optimum of previous model

Suppose I have the following model:

library(glmmTMB)
Owls <- transform(
  Owls,
  Nest = reorder(Nest, NegPerChick),
  NCalls = SiblingNegotiation,
  FT = FoodTreatment
)
fit_zipoisson <- glmmTMB(
  NCalls ~ (FT + ArrivalTime) * SexParent + offset(log(BroodSize)) + (1 | Nest),
  data = Owls,
  ziformula = ~1,
  family = poisson
)

Now suppose I want also a random slope for SexParent.

I want to start the new optimization at the optimum of the previous model (I found out this allows convergence for more complex models).

I tried:

glmmTMB_get_optimum <- function(mod) {
  # Fixed effects
  starting_point <- list(
    beta = glmmTMB::fixef(mod)$cond,
    betazi = glmmTMB::fixef(mod)$zi,
    betadisp = glmmTMB::fixef(mod)$disp
  )
  #Random effects
  ref <- glmmTMB::ranef(mod)
  if (length(ref$cond) > 0) {
    starting_point$theta = ref$cond
  }
  if (length(ref$zi) > 0) {
    starting_point$thetazi = ref$zi
  }
  if (length(ref$disp) > 0) {
    starting_point$thetadisp = ref$disp
  }
  return(starting_point)
}

fit_zipoisson_2 <- update(
  fit_zipoisson,
  formula = NCalls ~
    (FT + ArrivalTime) *
      SexParent +
      offset(log(BroodSize)) +
      (SexParent | Nest),
  start = glmmTMB_get_optimum(fit_zipoisson)
)

But I get:

Error: parameter vector length mismatch: in 'start', length(theta)==1, should be 3

What am I doing wrong?

Solution

The proximal thing you're doing wrong is failing to distinguish between the random effects values (BLUPs/conditional modes/whatever else you want to call them) and the random effects parameters which are "theta" by default. (Internally the RE parameters for the conditional, zero-inflation, and dispersion components of the model are called theta, thetazi, and thetadisp.) If you examine ref$cond you'll see that it's a one-element list containing a matrix of conditional modes.

More generally, though, you would be in trouble anyway because the theta vector for the model with (1|Nest) only has one parameter, while that for the model for the (SexParent|Nest) model has three parameters (two log-standard-deviations and a correlation parameter).

Without going into all the details,

retrieve all the parameters

You can do this slightly more compactly (although with more messing around in internal structures) as follows:

glmmTMB_get_optimum <- function(mod) {
  allpars <- mod$obj$env$parList(mod$fit$par, mod$fit$parfull)
  ## drop empty components
  allpars <- allpars[lengths(allpars) > 0]
  ## drop conditional modes
  allpars <- allpars[!names(allpars) %in% paste0("b", c("", "zi", "disp"))]
  allpars
}
pars <- glmmTMB_get_optimum(mod)

augment the theta vector to make it the right length. In this case c(theta_orig, 0, 0) will probably work OK. theta_orig is the estimated log-SD of the intercept; the second parameter is the log-SD of the slope; the third is a correlation parameter. See the covariance structure vignette for details of the parameterization.

pars$theta <- c(pars$theta, 0, 0)