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

Question

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
#head(round(lambdas,2))

# 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
#data=cbind(obj,time,y,x,c(t(lambdas))); head(round(data,2),2*q+2)

#============================== 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[1]
# 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[1]
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.

Answer 1

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)))
  cat("grad:\n"); print(grad <- grad(devfun, unlist(pars)))
  cat("scaled gradient:\n")
  print(scgrad <- solve(chol(hess), grad))
}
## 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)

Further reading: https://en.wikipedia.org/wiki/Conditional_logistic_regression

Question 2

If you only look at each object until it "dies", you have one observation per object. It's survival time. This means that frequentist methods can't estimate any object-level effects. You sometimes have observation level random effects in over dispersed counting data, but those are basically just data-transformations and don't apply here.

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.