Multifractional polynomial with mixed-effects logistic model

Question

Dataset

I am working on a dataset collected from more than 20 hospitals. I used the glmer command to estimate a mixed-effects logistic regression model with the following as patient-level predictors:

• age (patient's age; a continuous variable),
• cci (Charlson Comorbidity Index; a continuous variable),
• sex (patient's sex; man vs. woman),
• race (patient's race; Asian vs. Black vs. South Asian vs. White), and
• htn (has hypertension vs. no hypertension).

I used the hospital variable (the hospital that recruited patients) to estimate the random intercept.

Problem

My code for mixed-model (worked very smoothly):

mixed_died30in <- glmer(died30in ~ age + cci + sex + race + htn + infx_source + (1 | hospital), data = passem, family = binomial, control = glmerControl(optimizer = "Nelder_Mead"), nAGQ = 10)

Then I used the mfp command:

mfp_died30in <- mfp(glmer(died30in ~ fp(age, df = 4, select = 0.05) + fp(cci, df = 4, select = 0.05) + sex + race + htn + infx_source + (1 | hospital), family = binomial, data = passem, control = glmerControl(optimizer = "Nelder_Mead"), nAGQ = 10))

The following error appeared:

Error in fixed.only && random.only : invalid 'x' type in 'x && y'
In addition: Warning message:
In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model is nearly unidentifiable: large eigenvalue ratio
 - Rescale variables?

Questions

• Does the mfp command from the mfp package in R support mixed-effect models?
• Is (are) there any problem(s) in my code?

Answer 1

You have both an error:

Error in fixed.only && random.only : invalid 'x' type in 'x && y'

For the former, you could check that both age and cci are numeric rather than factor objects? I can't see any other clues as to this error message, which seems to crop up in lots of different settings.

And a warning:

In addition: Warning message:

In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :

Model is nearly unidentifiable: large eigenvalue ratio

- Rescale variables?

This may be due to age in particular needing centering/scaling before analysis? Though apparently rescale=TRUE is the default option.

Answer 2

If your primary interest is in modeling the continuous variables as smooth, potentially nonlinear effects, you might try fitting the model as a GAMM instead. mgcv can implement a similar model to the one you're trying to fit with mfp() by using its random effect capabilities. Here's an example of a similar model:

library(mgcv)
library(gamair) #for wesdr dataset

data(wesdr) #diabetic retinopathy dataset, binary outcome
set.seed(42)
#Make up fake random effect ids for 20 hospitals
wesdr$hospital <- factor(sprintf("H%i", floor(20*runif(dim(wesdr)[1])))) 
wesdr$is_male <- round(runif(dim(wesdr)[1])) #Fake variable for sex

#number of knots (k) chosen using the approximate rule of sqrt(# unique observations of data)
# bs="cr" indicates cubic regression splines
mod <- gam(ret ~ s(dur, bs="cr", k=16) #Smooth effect for duration of disease
           + s(bmi, bs="cr", k=12) #Smooth effect for bmi
           + s(gly, bs="cr", k=12) #Smooth effect for glycosylated hemoglobin
           + s(hospital, bs="re") #Random effect for hospital
           + is_male, #fixed effect for sex
           family = binomial(), data=wesdr, method="REML")
summary(mod)
plot(mod, pages=1, shade=TRUE)

As long as you use "enough" knots the exact number is not very important, thanks to the REML-based smoothness penalty.