4
$\begingroup$

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.
$\endgroup$
12
  • $\begingroup$ Related: stats.stackexchange.com/questions/593616 $\endgroup$ Commented Oct 26, 2022 at 14:21
  • $\begingroup$ "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? $\endgroup$
    – statmerkur
    Commented Oct 27, 2022 at 13:32
  • $\begingroup$ @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. $\endgroup$ Commented Oct 27, 2022 at 13:38
  • $\begingroup$ 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... $\endgroup$
    – statmerkur
    Commented Oct 30, 2022 at 17:19
  • $\begingroup$ ... (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. $\endgroup$
    – statmerkur
    Commented Oct 30, 2022 at 17:21

1 Answer 1

3
+50
$\begingroup$

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.

$\endgroup$
6
  • $\begingroup$ Thank you for you answer! It includes quite a few useful pointers and ideas. Regarding Q1, for starters it is encouraging to learn that Your model is basically well specified and identifiable. What I am looking for is sets of parameter values (including sample size and number of time periods) that would make convergence easier, if that is at all possible. Regarding Q2, I am following Tutz & Schmid (2016) chapter 9 on mixed effects survival models. So are you sure that frequentist methods can't estimate any object-level effects? $\endgroup$ Commented Nov 6, 2022 at 11:08
  • $\begingroup$ The logic of If you only look at each object until it "dies", you have one observation per object. This means that frequentist methods can't estimate any object-level effects. is very appealing. Yet I am reluctant to conclude that an entire chapter of an otherwise sensible textbook by seemingly expert authors is nonsense. I wonder what is missing here... $\endgroup$ Commented Nov 6, 2022 at 11:23
  • $\begingroup$ I meant I am reluctant to conclude it right away, as I suspect there might be other reasons behind the problem. I did not mean to say I would never do it. $\endgroup$ Commented Nov 6, 2022 at 12:07
  • $\begingroup$ Sorry i accidently hit enter and then got stuck editing, I'm still looking at frailty models $\endgroup$ Commented Nov 6, 2022 at 12:12
  • $\begingroup$ It is of course sensible to not discount a reputable source over my Statements. I have looked at a variety of sources about frailty in survival models but i can't find anything that is on the individual level. for Example here there is this data.princeton.edu/pop509/frailtyr but all the moms in their data had multiple kids. Sadly i don't have access to Trutz and Schmidt right now, but perhaps their REs had some structure to them as well. $\endgroup$ Commented Nov 6, 2022 at 12:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.