Model selection for Mixed Effects Regression - Correcting heteroscedasticity

Question

My experiment consist on measuring the root mean squared error (RMSE) between a path and what the participant produces. To produce the trace each participant experiences 3 feedbacks (1st factor - categorical), and each feedback has a bandwidth of either 2 or 5 Hz (2nd factor - numerical). Finally, each trace can have a mean value of 10 or 15 (3rd factor - numerical). Each participant experiences all the possible combinations (3 feedbacks * 2 frequencies * 2 mean levels = 12 combinations) twice. Hence, there are 24 RMSE measurements per participants. I have 7 participants.

The model I am fitting is

model_error <-lme(error ~ feedback * frequency * force_level * motor_units,
               random = ~ 1 | subject, data=df)

I am using the nlme library in R to implement it and emmeans to test for differences between the marginal means.

When checking for the soundness of the model, as expected from a dependent variable with a lower bound (i.e. non negative), the residual plot presents heteroscedasticity (see image below).

After using the log transform for the error and fitting the model again, the situation normalises.

What bothers me is that when looking for significant differences, such as

emmeans(model_error, pairwise ~ feedback)

some results change drastically, going from p-values of 0.0001 to 0.3.

Can someone please help me on understanding whether the log transform is necessary and how to proceed in any case?

UPDATE

I also tried the Box-Cox transform with lambda = -0.6 (the value was given as the optimal lambda by the boxcoxmix package). The QQ plot looks even nicer with respect to the log transform, but still: which results from emmeans should I believe? Before or after the transformation? Why?

UPDATE 2 These are the residuals after Box-Cox transformation

Answer 1

My suggestion would be to not focus so much on asterisks and instead focus on getting the best possible understanding of the results you have. It appears that the model with the Box=-Cox transformed response is a much better fit, so that is probably what you should go with.

A model with a transformed response is a different model than the one with the same right-hand side but no transformation. So it shouldn't be surprising that those two models yield different predictions and that those predictions compare differently.

What you do for post-hoc analysis depends on which interactions are still in play. Going straight to doing marginal comparisons of the levels of feedback is likely a mistake -- that is appropriate only if there are no strong interactions with feedback in the final model. For starters, you should something like this:

require(car)
Anova(final_model)

require(emmeans)
emmip(final_model, feedback ~ frequency * force_level | motor_units)

(where final_model is the one with the Box-Cox-transformed response). The Anova results provide an overview based on hierarchical tests of model effects. The emmip plot provides a more subjective view. You may follow up with emmip() runs with some of the factors excluded, if they aren't involved in interactions.

Let's suppose, say, that feedback interacts with frequency but interacts very little with the other two factors. Then it would be appropriate to do post-hoc tests of feedback separately for each frequency -- don't do marginal comparisons of feedback. This can be done via

emm = emmeans(final_model, ~ feedback | frequency)
emm          # view the EMMs
pairs(emm)   # compare the emms

Finally, on the transformation. It is often simpler to explain things on the original response scale. It is possible to back-transform the results to the original response scale:

emm = update(emm, tran = make.tran("boxcox", -0.6))
confint(emm, type = "response")

For pairwise comparisons, I suggest sticking with the pairs() output above, because the main purpose there is testing the comparisons, not quantifying them on the response scale.

There are several other considerations having to do with interactions and transformations. Look at a good experimental design text, and for the emmeans implementations, read the vignettes in that package.

Answer 2

From the plots in the OP it is clear that the original model has a heteroscedastic residuals.

This means that p-values and confidence intervals will likely be biased.

The usual way to proceed is to log-transform the response variable, which improves matters substantially, although there still appears to be larger variance with increased fitted values.

It is unsurprising that inference on the log-transformed model is different from the original model, since it is a different model.

The Box-Cox transform has successfully removed the heteroscedasticity.

The downside to proceeding with the Box-Cox transformed model is that it makes interpretation of the model coefficients much more difficult. With the untransformed, or log-transformed model, there is a very obvious interpretation, but this is not usually so after applying the Box-Cox transform. There is presumably some underlying theory of how the data are generated, so a good question to ask is whether these transform are consistent with that theory? Missing/unmeasured variables could also lead to heteroscedasticity and making arbitrary transformations in the pursuit of p-values may be very misguided.

An alternative approach is not to be too concerned with p-values, and instead just go with the log-transformed model and focus on effect sizes, not p-values. Another alternative is to use a heteroscedastic-robust estimator such as Huber-White. In this case it may not be possible to fit random intercepts, so another way to control for non-independence is to fit fixed effect for subjects instead.