Fitting problem only when using random effects

Dear Maaslin2 developers,
Thanks for your great tool!
I’m getting a weird error- WARNING::Fitting problem for feature X returning NA for all features when using random effects.

fit_data <- Maaslin2(
  input_data, input_metadata, outdir, transform = 'Log', max_significance = 0.2,min_abundance = 0, min_prevalence = 0 , normalization = 'None',
  fixed_effects = c('Status'), random_effects = c('BMI', 'Age')  ,standardize = FALSE,
  cores=1)

Which does not happen when not using random effects-

 fit_data <- Maaslin2(
  input_data, input_metadata, outdir, transform = 'Log', max_significance = 0.2,min_abundance = 0, min_prevalence = 0 , normalization = 'None', fixed_effects = c('Status','BMI', 'Age'), random_effects = c()  ,standardize = FALSE,
  cores=1)

Any idea why?

This is probably caused by the fact that BMI and age are (I assume) continuous. I think you’re probably not correctly specifying the model you want to fit. Unless you have some unusual experimental design where you have multiple measurements for each combination of a discrete set of BMI & age categories, I think you probably mean to include those as fixed effects.

Thanks Andrew! I agree that these should be included as fixed effects, but wanted to test the results while using random effects.
The problem is that even when using categorical data for either BMI or age (or both), the fitting returns NA anyway, and I try to understand why.

when using categorical data for either BMI or age (or both)

Did you manually bin BMI/age into a small number of categories or just use something like input_metadata$BMI = as.character(input_metadata$BMI)? Because the latter won’t work.

I think Maaslin2 is using lmerTest::lmer() under the hood in your case. That function checks if the number of group levels is < the number of observations and throws an error if not, which bubbles up to the Maaslin2 error you saw.

You could try using cut() if you’d like to discretize BMI/age. Example:

dat = data.frame(x = rnorm(100))

dat$y = rnorm(100, mean = dat$x)

m <- lmerTest::lmer(y ~ (1|x), data = dat) # Doesn't work

dat$z = cut(dat$x, breaks = 4)

m <- lmerTest::lmer(y ~ (1|z), data = dat) # Works