[R-sig-ME] model convergence warnings with default optimizer but not with L-BFGS-B

[Steve Denham]

Hi Alessandro,
I would be concerned that you have a partial confounding issue due to fitting change from baseline AND using baseline as a continuous covariate.  See Frank Harrell's comments regarding analysis of repeated measures at https://biostat.app.vumc.org/wiki/Main/ManuscriptChecklist, especially regarding the use of change scores in parallel design analyses.
Steve Denham Senior Biostatistics Scientist, Charles River Laboratoies 

    On Tuesday, May 25, 2021, 10:16:14 AM EDT, Noci, Alessandro via R-sig-mixed-models <r-sig-mixed-models using r-project.org> wrote:  
 
 Hi all,

I am running some experiments where I need to fit a mixed-effects model for
repeated measures (MMRM) for two arms clinical trials. The outcome is the
mean change from baseline of a given response variable; the model assumes
unstructured mean (i.e. different mean across groups and across time) and
unstructured covariance which is assumed to be different for the two
groups. Moreover the model contains the baseline outcome as an additional
predictor. As discussed in another post (
https://stat.ethz.ch/pipermail/r-sig-mixed-models/2020q4/029135.html) I can
use the package glmmTMB for the model fit when assuming different
covariances for the two arms.

The issue I am encountering is related to the model convergence warnings (a
reproducible example is provided below). Fitting the model using the
default nlminb optimizer leads to the following convergence warning:

Warning message:In fitTMB(TMBStruc) :  Model convergence problem;
false convergence (8). See vignette('troubleshooting')


If I fit the same model on the same data using the L-BFGS-B optimizer I do
not receive any warning and the model convergence is detected as
successful. However the results from the model fit look very similar:
indeed regression coefficients, covariance matrices and log-likelihood
value are the same (at least up to the third decimal digits). Also, when
generating the data with different seeds or changing some parameters used
for the simulation (e.g. the sample size, population mean trajectories..),
I still obtain this issue; so I don't think this is a problem related to
the specific generated dataset.

I also tried to initialize the optimization problem close to the optimum,
i.e. with the results coming from the L-BFGS-B optimizer. I obtain that (1)
I get the same convergence warning and (2) I don't improve computational
time (I expect that if I start very close to the optimal value, I do have
convergence and I also reduce the computational time).

My questions are: why do I get the convergence warning only for the nlminb
optimizer but the optimal value seems almost identical to the results from
the L-BFGS-B optimizer? Why if I initialize the problem close to the
optimal value I do not improve computational times? Am I supplying the
starting values correctly?

Thank you very much.
Best regards,
Alessandro


library(tidyverse)
library(glmmTMB)
library(MASS)
library(lme4)

# set-up simulation
set.seed(123)
N = 300
n1 = n2 = N # equally sized arms

time = seq(0,12,by=2) # time variable

# create covariance matrix
sd_intercept = 10
sd_slope = 5
cor_slope_inter = 0.5

sd_error = 6

covRE = matrix(c(sd_intercept^2,cor_slope_inter*sd_intercept*sd_slope,
                cor_slope_inter*sd_intercept*sd_slope,sd_slope^2),ncol=2)
Sigma =
cbind(1,time/12)%*%covRE%*%rbind(1,time/12)+diag(sd_error^2,nrow=length(time))

# mean trajectory control
mu2 = 10+12/12*time

# mean trajectory intervention
mu1 = 10+6/12*time

m = length(time)

# simulate data
r2 = data.frame("patnum" = rep(1:n2,each=m),
                "group" = "Control",
                "visit" = rep(0:(m-1),n2),
                "time" = rep(time,n2),
                "x_bl" = NA,
                "x" = c(t(mvrnorm(n=n2,mu=mu2,Sigma=Sigma))))
# Intervention group
r1 = data.frame("patnum" = n2+rep(1:n1,each=m),
                "group" = "Intervention",
                "visit" = rep(0:(m-1),n1),
                "time" = rep(time,n1),
                "x_bl" = NA,
                "x" = c(t(mvrnorm(n=n1,mu=mu1,Sigma=Sigma))))
# Pool both groups and add baseline
r = rbind(r2,r1)
r$x_bl = rep(r$x[time==0],each=m)
r$group = factor(r$group,levels=c("Control","Intervention"))

data = as.data.frame(r)

# add change from baseline
data = data %>%
  group_by(patnum) %>%
  mutate(x_change_bl = x - x_bl)

data = subset(data, visit != 0)
data$visit = as.factor(data$visit)
data$time = as.factor(data$time)

# create dummy variable related to the two treatment arms
data$g_ref = lme4::dummy(data$group, "Control") # reference arm
data$g_treat = lme4::dummy(data$group, "Intervention") # treatment arm

# set model formula
formula = x_change_bl ~ x_bl + visit*group + us(0 + g_ref:visit | patnum) +
us(0 + g_treat:visit | patnum)

# fit model
init = Sys.time()
fit_mmrm = glmmTMB::glmmTMB(formula,
                            data = data,
                            dispformula = ~0,
                            REML = TRUE)
end = Sys.time()
end-init

# extract regression coefficients
betas = glmmTMB::fixef(fit_mmrm)$cond

# log likelihood
loglik = logLik(fit_mmrm)

# extract Cholesky decomposition of the covariance matrices
theta = getME(fit_mmrm, "theta")

# Change optimizer
control = glmmTMB::glmmTMBControl(optimizer = optim,
                                  optArgs = list(method="L-BFGS-B"))
init = Sys.time()
fit_mmrm = glmmTMB::glmmTMB(formula,
                            data = data,
                            dispformula = ~0,
                            REML = TRUE,
                            control = control)
end = Sys.time()
end-init

# extract regression coefficients, log-likelihood and Cholesky
decomposition of covariance matrices
betas_2 = glmmTMB::fixef(fit_mmrm)$cond
loglik_2 = logLik(fit_mmrm)
theta_2 = getME(fit_mmrm, "theta")

# check equality of results
all.equal(betas, betas_2, tolerance=1e-3)
all.equal(loglik, loglik_2, tolerance=1e-3)
all.equal(theta, theta_2, tolerance=1e-3)

# fit model
init = Sys.time()
fit_glmm = glmmTMB::glmmTMB(formula,
                            data = data,
                            dispformula = ~0,
                            REML = TRUE,
                            start = list("beta" = as.numeric(betas_2),
                                        "theta" = theta_2))
end = Sys.time()
end-init

    [[alternative HTML version deleted]]

_______________________________________________
R-sig-mixed-models using r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
  
	[[alternative HTML version deleted]]