I started with this initial model:
m1 <- lmer(response ~ treatment + (1|subjectID), data = data)
However, the residuals of the model are heavy-tailed (presumably enough to violate the normality of residuals assumption?).
So, I tried a log transformation:
m2 <- lmer(log(response) ~ treatment + (1|subjectID), data = data)
Unfortunately, the residuals are virtually the same:
I also tried sqrt, cube root etc., with not much luck.
Finally, I tried a linear quantile mixed model, as it makes no assumptions about normality of residuals:
m3<- lqmm(response ~ treatment, random = ~1, group = subjectID, data = data)
Since m1 and m2 seem to violate the assumptions, I would have thought that lqmm (m3) is the most appropriate model. However, when comparing AICs, m2 (linear mixed model with log transformation) is HUGELY lower than m1 and m3:
AIC(m1, m2, m3)
I suppose my question is - would you choose the lqmm (m3) regardless of AIC, since m1 and m2 are arguably both invalid due to the non-normality of residuals? Or, would the considerably lower AIC of m2 give you 'pause' to suspect that it is the more appropriate model, despite deviation of residuals from normality?
For what it's worth, p values for one treatment group (compared to the reference group) are the same regardless of which model I use, but differ for two of the treatment groups compared to the reference group; m3 gives a significant p value for both of these comparisons, while m1 and m2 do not. In other words, the group 1 - group 4 comparison is always significant regardless of the model, but group 1 - group 2 and group 1 - group 3 comparisons are only significant in m3.
Apologies for any misunderstanding I may have and thank you for your help! I am not sure if I can provide the data due to confidentiality. Even if I provide only dummy variable names to ensure complete anonymity, I have to repeat this process many times over for different response variable and am therefore more interested in the conceptual reasoning for choosing one model over another, rather than an actual answer in this case.
Thank you again. :)
As requested, here are the residual vs fitted plots for m1 and m2:
Additionally, here are the QQ plots for the response variable split by group:
The '2 hour' group has more observations than the other three groups as all participants had a 2-hour observation, while each participant was randomised to one of the other groups (essentially within-subjects comparison when comparing 2 hours to all other groups and between-subjects when comparing the other three groups to each other, hence the random effect). There are 104 total observations across 55 participants.
From these plots, it looks like only the two-hour group is an issue - does that change the interpretation?