Effect size derived from LME longitudinal model: the statistical findings projected back down onto a group of people

Question

I have been studying the change in a metric X with a linear mixed effect model. I have built this model in a multivariate setting, so I can see how each of my covariates (Time, sex, age) affect X. Subsequently, I assessed how the time-differences were impacted by sex and age using simple interaction models.

I paste here my code and results here:

mod1.lme <- lme(X ~ Time + bmi+ sex+ smoker + drinker + diabetes + hypertension + race , random = ~ 1|ID, data = df1, na.action = na.omit, method = "ML")

summary(mod1.lme)

mod3.lme <- lme(X ~ Time + sex + Time * sex , random = ~ 1|ID, data = df1, na.action = na.omit, method = "ML")

summary(mod3.lme)

Then I have substracted the mean to the metric BMI, because I was told this would eliminate the regression to the mean for this metric

df1$bmi_acc = df1$bmi-df1$mean_bmi

and assessed its longitudinal impact like this (maybe I should have done it with an interaction model?)

THE IMPORTANT QUESTION IN THIS POST:

From these results I could say that there is change over time in the metric X (based on mod1), and that there are several variables affecting this the variable X crossectionally, but after testing them all in interaction models, I see that only sex (women(==0) progress more than men(==1)) and people with high BMI have a higher progression over time.

EFFECT SIZE: For example, given these results, I would like to know if I can say (based on my cofficients) that for example women with high BMI would progress at xy% per year, while a an average men with BMI xy, and a starting X of xy, would progress at xy% per year.

The way I normally interpret these coefficients is: (based on mod1) For every year, you will increase 1.72 in X (so 581.48 + 1.72). The other covariates in that model are crossectional and I am not sure how to interpret them in this context.

I am confused though as to how to include the interaction models in this context. Shall I forget about mod1 and focus on the specific interaction models?

My interpretation:

SEX: Based on model 3:

For every increase in a year, you will suffer an increase of 2.76 in X, that if you are a woman, will be 1.78 in addition to men. So the final model will be:

• Women: 687.23 + 2.76 +1.78 = 691.77 (being a woman you would expect an increase in X of 0.66% per year)
• Men: 687.23 + 2.76 = 689.99 (being a man you would expect an increase in X of 0.40% per year)

For BMI (based on model 2)

where every increase in a BMI unit would increase X by 3.79...

Going from a BMI of 25 to 31 Kg/m2 in one year would mean:

6*3.79= 22.7 697.89 +22.7 = 720.59 (an increase of 3.25%)

if we merge this with the previous context. We would say:

• Women: 0.66% + 3.25% = 3.91%
• Men: 0.40% + 3.25% = 3.65%

Is this all right?? I am trying to make my results easy to read for clinicians.

Thanks!!

Answer 1

With this scale of a data set, you should be combining all your predictors into a single model. Doing multiple mini-models runs a risk of omitted-variable bias, not allowing you to take the omitted predictors into account in each mini-model.

The best approach is explained for example in Harrell's course notes, in particular Chapter 4 on multivariable modeling strategies. Consider how much data you have, decide on how many predictors that might allow you to evaluate, and build a rich model that captures as much of the data as is reasonable without overfitting.

You have on the order of 5000 observations total, with a bit under 2 observations per individual on the average. With your continuous outcome you could probably accommodate 100 or more predictors, counting each interaction term as a predictor. It seems that you could easily incorporate all potential 2-way fixed-effect interactions into your model, with a random intercept accounting for (estimated) baseline differences among individuals. Why not do so, at least to start? If you want a simpler model that is "easier to read for clinicians" you can always pull back from that, trading off the loss of predictive ability against the parsimony you achieve.

In terms of centering continuous variables like BMI, that will make interpretation of the intercept and the coefficients of predictors with which it interacts easier to interpret, but it won't affect any predictions from your model.

In terms of interpreting your model, you are predicting individual values of X as a function of all of your predictors, including Time. So you need to evaluate interactions to see whether there are combined effects of individual predictors that are different from what you would expect from each predictor individually. That's not only true for things like Time and sex as in one of your mini-models to see if the effect of Time on outcome depends sex, but also for BMI and sex if you think that the association of BMI with outcome depends on sex, and so on.

With the default treatment coding used in R, the interpretation of coefficients and interaction terms is straightforward if tedious. The intercept is the estimated value when all continuous predictors have values of 0 and all categorical predictors are at their baseline levels. The individual coefficient for a predictor involved in interactions is its extra association with outcome when all its interacting predictors are at 0 or reference values. Two-way interaction coefficients represent the differences in outcome between what you would predict from the individual coefficients and what you get with both together.

In response to comment:

The principles in the previous paragraph can be applied to your model mod3.lme. With female as the reference condition (sex = 0), then we can say the following:

The intercept of $687.2$ is the estimated outcome for a female at Time =0.

For a female, the outcome increases by $2.77$ per unit of Time, the coefficient for Time. (I think you might have misinterpreted the Time coefficient; it's for the reference "female" condition of the interacting sex predictor, while you seem to have taken it to be the coefficient for "male.")

For a male at Time = 0, the estimated outcome is $687.2 + 17.3 = 704.5$ (intercept plus coefficient for sex = 1). (Men start at that level at Time = 0, not at $687.2$ as you seem to imply in part of the question.)

For a male, the estimated outcome increases by $2.77-1.78=0.99$ per unit time, the Time coefficient plus the Time:sex interaction. (You were correct that the change per unit time is $1.78$ greater for females than males according to the interaction term, but see above for the misinterpretation of the Time coefficient itself.)

You apply the same interpretations to all predictors involved in interactions. I find it's best to think things through one coefficient at a time, asking just what that coefficient represents and the corresponding value(s) of the predictor(s) by which to multiply the coefficient, then adding all the terms up. That's tedious, but it's the best way I know to avoid errors. Or you can use software designed to make such estimates easier. Many find the emmeans package in R to be a help.