# Simulating and estimating generalized mixed effects models: how to avoid convergence problems?

Basic generalized mixed effects
I am simulating data from a generalized mixed effects model $$P(y_{it}=1\mid x_{it})=h(\delta_0+\alpha_i+\beta_t+\gamma x_{it})$$ where the dependent variable is binary, $$y_{it}\in\{0,1\}$$, and where the inverse link function is logistic, $$h(z)=\frac{1}{1+e^{-z}}$$. I set $$\alpha_1=\beta_1=0$$ for identifiability.

I then fit the same model on the simulated data. (Well, not exactly. I do not force $$\hat\alpha_1=\hat\beta_1=0$$. I think $$\hat\beta_1=0$$ is forced automatically, though. When R discovers that the model matrix suffers from multicollinearity, it kicks out some terms to get rid of it, and the term corresponding to $$\beta_1$$ is one of them.) For some configurations of parameters $$\alpha$$ and $$\beta$$, sample size (number of individuals) $$n$$ and the number of time intervals $$q$$, I tend to get decent estimation results. For others, I get convergence problems and/or poor estimates.
Question 1: What should I tweak in the simulation to avoid convergence problems and to obtain better estimation results?

Generalized mixed effects for modelling survival
I also try a variation of the same setup. After generating the data, I examine each individual $$i=1,\dots,n$$ separately. I discard all data points after the first $$y_{it}=1$$ is recorded. E.g. if $$q=5$$ and $$y_{i,\cdot}=(0,0,1,0,1)$$, I discard the last two points and only retain $$\tilde y_{i,\cdot}=(0,0,1)$$. The subject-matter interpretation is that I am generating survival data where $$0$$ denotes survival and $$1$$ denotes death in a given period. Thus the individual above survived the first two time periods and died in the third one. When I try estimating the model based on $$\tilde y_{it}$$, I always get either convergence problems or a perfect fit (or both). I have tried various parameter configurations, but I never got decent estimates.
Question 2: What should I tweak in the simulation to avoid convergence problems and to obtain better estimation results?

R code:


#============================== Simulation

# Set the number of individuals and the number of time periods
n=100; q=40

# Generate parameter values:
set.seed(1); alphas=runif(n=n  ,min=-1,max=1) #-0.025*q
set.seed(2); betas =runif(n=q  ,min=-1,max=1)
gamma =1
delta0=-0.1 # the intercept

# Generate covariate x
set.seed(3); x     =runif(n=n*q,min=-1,max=1); X=matrix(x,nrow=n,byrow=TRUE)

# Specify the desired hazard rates
lambdas=matrix(NA,ncol=q+1,nrow=n)
for(i in 1:n){
for(t in 1:q){
lambdas[i,t]=plogis(delta0+alphas[i]+betas[t]+gamma*X[i,t]) # plogis yields a logistic transformation
}
}; lambdas[,q+1]=1 # by definition

# Obtain the binary representation of survival
Y=matrix(NA,ncol=q,nrow=n)
for(i in 1:n){
for(t in 1:q){
set.seed(n*i+1e6*t); Y[i,t]=rbinom(n=1,size=1,prob=lambdas[i,t])
#if(Y[i,t]==1) break; # !!! If this line is commented out, we get the basic GLM. Otherwise, we get the survival model !!!
}
}
y=c(t(Y))

# Create identifiers of individual (obj) and time (time)
obj=c(1,rep(0,q-1)); obj=rep(obj,n); obj=cumsum(obj)
time=c(1:q); time=rep(time,n)

# For each individual, delete observations after death
subset=which(!is.na(y)); y=y[subset]; x=x[subset]; obj=obj[subset]; time=time[subset]

# Print out some data

#============================== Estimation

# We will use timeout to prevent the calculation from taking forever and R hanging up
tout=20

# Estimate a random effects (RE) model using mgcv package
m1=R.utils::withTimeout( mgcv::gam(y~s(factor(obj),bs="re")+factor(time )+x, family="binomial"), timeout=tout )
summary(m1)
delta0_hat1=m1$coef # How to extact estimates of individual random effects from this model object? # Extract estimates of time fixed effects betas_hat1=coef(m1)[2:q] cor(betas[-1], betas_hat1 )^2 # Plot them against the true parameter values plot(x=betas[-1],y=unlist(betas_hat1),xlab="true",ylab="fitted",main="Time fixed effects"); abline(a=0,b=1) # Because of the multicollinearity induced by individual and time effects, the estimates may be biased by a constant. # Estimate a random effects (RE) model using lme4 package m2=R.utils::withTimeout( lme4::glmer(y~(1|obj)+factor(time)+x, family="binomial"), timeout=tout ) # This does not converge for me if q=40 but does converge if q=20. You can set q at the top of this script. summary(m2) # Extract estimates of individual random effects and time fixed effects #delta0_hat2=as.numeric(coef(m2)$$obj[1,-c( q+1)]) ??? alphas_hat2=unlist(lme4::ranef(m2)) betas_hat2 =as.numeric(coef(m2)$$obj[1,-c(1,q+1)]) cor(alphas ,alphas_hat2)^2 cor(betas[-1],betas_hat2 )^2 # Plot them against the true parameter values dev.new(); mar1=c(4,4,3,0.5); par(mfrow=c(2,1),mar=mar1) plot(x=alphas ,y=alphas_hat2,xlab="true",ylab="fitted",main="Individual random effects"); abline(a=0,b=1) plot(x=betas[-1] ,y=betas_hat2 ,xlab="true",ylab="fitted",main="Time fixed effects"); abline(a=0,b=1) par(mfrow=c(1,1)) # Because of the multicollinearity induced by individual and time effects, the estimates may be biased by a constant. # Estimate a fixed effects (FE) model m3=R.utils::withTimeout( glm(y~factor(obj)+factor(time)+x, family="binomial"), timeout=tout ) summary(m3) # Extract estimates of individual fixed effects and time fixed effects delta0_hat3=m3$$coef alphas_hat3=m3$$coef[2:n] betas_hat3 =m3$coef[(n+1):(n+q-1)]
cor(alphas[-1],alphas_hat3)^2
cor(betas[-1] ,betas_hat3 )^2
# Plot them against the true parameter values
dev.new(); mar1=c(4,4,3,0.5); par(mfrow=c(2,1),mar=mar1)
plot(x=alphas[-1],y=alphas_hat3,xlab="true",ylab="fitted",main="Individual fixed effects"); abline(a=0,b=1)
plot(x=betas[-1] ,y=betas_hat3 ,xlab="true",ylab="fitted",main="Time fixed effects"); abline(a=0,b=1)
par(mfrow=c(1,1))
# Because of the multicollinearity induced by individual and time effects, the estimates may be biased by a constant.

• Oct 26, 2022 at 14:21
• "For some configurations of parameters $\alpha$ and $\beta$, sample size (number of individuals) $n$ and the number of time intervals $q$, I tend to get decent estimation results. For others, I get convergence problems and/or poor estimates." Can you be a bit more specific here? When do get convergence problems and/or poor estimates? Oct 27, 2022 at 13:32
• @statmerkur, actually, I have mostly kept $\alpha$ and $\beta$ as they are except for shifting their locations (adding or subtracting up to 2.5). I have also varied $n$ (100, 200, 1000) and $q$ (10, 20, 40, 100). For quite some combinations I got poor results. If you check out my code with the settings there, models m1 and m3 work OK but m2 fails. For some other combinations, models m1 and m3 also fail. And they all fail almost universally if I do survival rather than basic GLM. Oct 27, 2022 at 13:38
• I don't have time for a full answer at the moment. But: (i) I think if you set method = "ML" in mgcv::gam the objective functions of m1 and m2 are identical, but different approximations are used to calculate the MLEs; (ii) In the GLM m3 we have the same issues regarding the model matrix and parameters as in the Q you linked to (but now $\alpha_1=0=\beta_1$). In the GLMMs only $\beta_1=0$ is forced; (iii) The way in which the $\alpha_i$ are penalized in mgcv::gam and lme4::glmer is equivalent to assuming that the $\alpha_i$ are (iid) normally distributed with mean zero... Oct 30, 2022 at 17:19
• ... (iv) the problem when you "discard all data points after the first $y_{it}=1$ is recorded" is probably complete separation, which you could double-check with the detectseparation package. Oct 30, 2022 at 17:21

# Question 1

There are a few things to say here. First of random effect coefficients are supposed to be (multivariate)normal. So(keeping the variance the same) you should write: set.seed(1); alphas=rnorm(n=n, sd = sqrt(1/3)) Source(among many) the original lme4-Paper, page 2: https://arxiv.org/pdf/1406.5823.pdf

But fixing this does not remove your convergence problem so i looked at the GLMM-FAQ(https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#convergence-warnings) which points to the R-command ?convergence(only works if you have lme4 loaded)

## What actually is your Problem

Keeping your code as is the convergence warning on m2 is:

optimizer (Nelder_Mead) convergence code: 4 (failure to converge in 10000 evaluations) Model failed to converge with max|grad| = 0.00857345 (tol = 0.002, component 1) failure to converge in 10000 evaluations

0.0085 is pretty close to the tolerance of 0.002 and checking the Hessian as instructed by ?convergence reveals no major deviations:

if (isLMM(m2)) {
pars <- getME(m2,"theta")
} else {
## GLMM: requires both random and fixed parameters
pars <- getME(m2, c("theta","fixef"))
}
if (require("numDeriv")) {
cat("hess:\n"); print(hess <- hessian(devfun, unlist(pars)))
}
## compare with internal calculations:
m2@optinfo$derivs summary(m2@optinfo$derivs$Hessian - hess)  So just keep the fitting process running for more than 10,000 steps and you get there: m2_restart <- update(m2, start=getME(m2, c("theta","fixef"))) #continues fitting from the m2 parameters m2_restart@optinfo$feval # 1692 steps
summary(m2_restart)


You might need multiple restarts, make sure to use the latest model sor the new start-parameters. Weirdly enough running it in one step did not work (I took the code from here: https://stackoverflow.com/a/19479820/7840119). But it actaully self-stopped at just over 12,000 iterations with the warnig-code of "ok".

m2_long <- glmer(y~(1|obj)+factor(time)+x, family="binomial", control = glmerControl(optCtrl=list(maxfun = 30000)))


### Why?

Your model is basically well specified and identifiable, but the huge number of coefficients and the low informational value of a binary response create a very large region of fairly flat likelihood that posses a considerable challenge to numerical optimization. This can also pretty unpredictable.

### How to fix this

I can think of multiple "solutions":

1. Ignore it and just accept the slight deviation from the true ML-estimator. The coefficients between m2 and m2_restart only differ by 0.02.Of course fitting 10,000 iterations almost everytime is slow
2. Embrace the slowness and fit to completion every time.
3. Simplify and here i actually have two ideas:

a) You can use a simpler approximation of the likelihood by setting the nAGQ parameter to 0:

m2_simplified <-  glmer(y~(1|obj)+factor(time)+x, family="binomial", verbose = 1, nAGQ = 0) # fitted in just 19 steps
summary(m2_simplified)


Check help("glmer") for what that actually means

b) If you don't like having your $$\alpha$$ being normal than you can use a conditional logit model. The downside is, that it doesn't "care" about estimating the $$\alpha$$ correctly, but it does respect the stratification:

library(survival)
m4 <- clogit(y ~ strata(obj)+factor(time)+x)
summary(m4)


You could maybe accomplish something by utilizing strong priors in stan_glmer but i will warn you: The rstanarm-Package takes a long time to install and fitting models with stan is very computationally expensive.