# Model selection for Mixed Effects Regression - Correcting heteroscedasticity

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 • @FransRodenburg Yes (it is written in the question). I also tried the Box-Cox transform with lambda = -0.6. In any case, the question still remains: which statistical results should I believe? With or without the transformation? Why? Dec 10, 2019 at 4:02
• @FransRodenburg Please have a look at update 2 with the new graph. I would say it definitely holds with the Box-Cox transform. My question is: what is the interpretation of why the p-values change? Am I measuring something different after transforming the outcome variable? Dec 10, 2019 at 4:29
• I agree that the Box-Cox transformed response looks much better in the diagnostic plot. The reason the $p$-value changes is that the model assumes the residual variance to be constant. If it is not, then neither is the standard error of the estimate. Dec 10, 2019 at 4:33
• You could also use a function for the variance, such as a power function: model_error <-lme(error ~ feedback * frequency * force_level * motor_units, random = ~ 1 | subject, weights = varPower(), data=df). This would avoid transforming the variables, which potentially decreases the interpretability of the coefficients. Dec 10, 2019 at 20:47
• Relevant for the question about why significance changes: When is performing back-transformation of inferences on transformed variables Ok, and when is it not Ok?. Notably for frequentist inference, because $f(\sigma^{2}_{x})\ne\sigma^{2}_{f(x)}$, except for $f(x)=x$ generally, and because both CIs and hypothesis tests rely on the estimated variance back-transformed inference on the transformed variable is not inference on the (untransformed) variable. Dec 11, 2019 at 16:53

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.

• Just to emphasize: That last part with the emmeans() stuff is just an illustration of the sorts of things you might consider. What is actually appropriate depends very much on what interactions you need to account for. Dec 11, 2019 at 16:28

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.

• I disagree that it is useful to understand the regression coefficients in a four-factor model. Even without a response transformation, that is a pretty convoluted enterprise. And I disagree on settling for the log model when the Box-Cox one is clearly the better model. Dec 11, 2019 at 16:25
• @rvl I am not saying that settling with the log model is the right thing to do. I am simply discussing the alternatives as I see them. A more parsimonious model can have it’s advantages. Without having good knowledge of the study, data and domain it is folly to say what is best, let alone suggest there is one right model. In any case, all models are wrong. It’s up to the analyst to decide what is right for them. I have seen several real examples where an analyst used a transformation in order to normalise residuals, when the actual problem was a missing variable, or missing nonlinear term... Dec 11, 2019 at 21:17
• .... (and it’s easy to construct such an example), resulting in completely misleading inference. Also, it might be the case the there are few or even no interesting interactions, in which case interpreting the coefficients would not be convoluted. Many times in statistical modelling there are trade-offs with interpretation, model fit, parsimony and other aspects. It's always worth considering the alternatives. Dec 11, 2019 at 21:19
• Good points, for sure. I can’t cover all of statistics in one answer and I did recommend further reading. Dec 11, 2019 at 22:03