1

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?
Abo7aneen
  • 85

2 Answers2

2

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.

James Stanley
  • 2,736
  • Thank you James, both are numeric (checked by is.numeric); I rescaled these variable and the warning still apearing – Abo7aneen Mar 03 '23 at 00:00
  • It looks like the Warning message is actually from glmer() rather than mfp(): see https://stats.stackexchange.com/q/412185/16974 and https://stats.stackexchange.com/a/231721/16974 – James Stanley Mar 04 '23 at 00:31
  • The second link is answered by Ben Bolker (author of lme4 package) and for that example he suggests there might be complete separation occurring. – James Stanley Mar 04 '23 at 00:32
  • As to the error: if it's coming from glmer, you could try a few things like varying the degrees of freedom for the continuous covariates. It's also very possible that mfp() might be incompatible with glmer(), since mfp() essentially wraps around the underlying model estimation function...[continued] – James Stanley Mar 04 '23 at 00:35
0

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.

jd_io
  • 1
  • 1