Issue with GxE_interaction_test and lme4 - same AIC/BIC for competing models

[RobbieRECAP]

Hi there,

Sorry, this is my first time posting an issue for an R package so I am unsure what is standard protocol. If you need me to submit a full markdown file then I could probably do that.

I am trying to use the GxE_interaction_test() command. I have combined data from 4 different cohorts. I have one measure of the environment and one measure of their temperament (the "genes" component for my analysis).

Originally I used GxE_interaction_test() without accounting for cohort effects (first image- this all looks fine to me). This is partially because I had standardised all the cohort's data before pooling so it had a mean of 0 and an SD of 1 but also because I couldn't get the lme4=TRUE option to work correctly. Now I have had feedback from peer reviewers where they are asking for a linear mixed model.

The problem I have is that when using LME4=TRUE I get the exact same AIC/BIC values for multiple competing models (second image), which I do not think should be possible. In particular, it looks like the main effects of the genes and environment are being included in all models as the intercept-only model/no interaction models have the same AIC/BIC values as the G+E model, and it seems to be saying that when I look at the coeficients for these models (third image).

Am I doing something wrong in the formulation of the formula? I have tried modifying the specification of the random effects- I would ideally only test for random slopes by cohort- and I have tried not including the covariates but I keep getting the same AIC for all models not including interactions when using LME4.

Any help would be much appreciated!

Rob

[AlexiaJM]

Hi Rob,

Yes this ok, there is no protocol, you can ask your questions here.

The issue comes is because you include "G*E" into your formula F. As the name state, you need to provide a "formula_noGxE", and here the "G*E" term causes R to include both G, G:E, and E terms. Hence why the models that add E or G to the formula don't change anything and all BICs are the same.

Here G*E means G + E + G*E so this what you use for the lme4 random effects, and you can keep it like this. Just make sure to remove the "G*E" from the fixed effects so
instead of "y tilda G*E + (0+G*E | blabla)"
do "y tilda 1+ (0+G*E | blabla)"
or "y tilda (0+G*E | blabla)".

[RobbieRECAP]

Hi Alex,

Thanks very much for your help, that did the trick! I do have another question because the results I am getting using the new methodology look a bit strange to me. In the original analyses, I could use LEGIT's GxE_interaction_test and then run a normal linear regression/produce simple slopes and Johnson Neyman plots and everything really looked like it matched up well. For example here is an original analysis that LEGIT thought demonstrated diathesis stress. This matches up well with the linear regression and the resulting Johnson-Neyman plot.

Original analysis using multiple regression:

New analysis using LME4 (random slopes for G, E, and G:E)
I have specified the random slopes as desired and the GxE_interaction_test is now suggesting evidence for differential susceptibility, with the crossover at 0.4. I find this very surprising because when I run a linear mixed model in LME4 with random slopes and then plot the simple slopes/Johnson-Neyamn plots then it looks like to me I have no crossover and so evidence for either E only or diathesis stress. Any ideas on why the conclusions seem to now differ so drastically? Anything I should look to change with the specification of random effects or how GxE_interaction_test decides which optimal model is supported?

Comparison of regression models (left original regression, right LME4 with random slopes)

One final thing I noticed about the new analyses: in the two differential susceptibility models (strong and weak), both are saying that about 60% of observations are below the crossover, I am not sure whether that is logical based on the fact that one is at 0.4 and the other is at -0.16. When I look at my environmental factor, 0.4 would match up well with 60% but not -0.16

[AlexiaJM]

Hi Rob,

Its a bit weird to have G, E, and G:E in the random effects when you have a different parametrization in the fixed effects, but I am not an expert in mixed models, so I can't say for sure if its okay. You could try with just G:E and E instead of G*E but i don't know.

It's also possible that the lme4 interfere with things somehow. What you could do in this case is remove lme4 models which have a non-significant GxE. There are ways to get p-values from these models given the lme 4 model fit. Or I would recommend using BIC, because AIC is way too lenient at allowing extra variables even when they contribute very little.

[RobbieRECAP]

Hi Alex,

I think your point about the potential problem of including a random slope when you do not have a corresponding fixed effect is important. I do not know whether this would influence the AIC/BIC calculations for these models (maybe it just subsequently influences the fixed effects and thus the estimated crossover point) but it is said here that it is not recommended. https://stats.stackexchange.com/questions/199377/is-it-reasonable-to-include-a-random-slope-term-in-an-lmer-model-without-the-cor

Would it be possible that in a future update of LEGIT the random effects also change depending on the model of interest? For example, the "E-only" model would only have a random component of (0+E|cohort_num) while differential susceptibility would have (0+G*E|cohort_num).

I think at the moment, if the LME4 component of LEGIT is still in an experimental phase, I will wait before using it in the manuscript. However, thank you very much for the quick responses which have helped a lot.