I would like to know what would be the key differentiator assessing change over time between a wilcoxon paired test and a linear mixed effect model with random intercept.

I did my analysis in R:

The wilcoxon paired test was a simple wilcox.test(df$T1.v1, df$T1.v2,data = df, paired = TRUE), where v1 stands for visit 1 and v2 stands for visit 2. While the linear mixed effect model was nlme(T1 ~ Time + sex + age, random=ID, data=df). Where time is a variable formed by 2 factors: V1 and V2.

In some cases I find significant differences in the wilcoxon paired test that do not stay significant in the linear mixed effect model. I know this can be right. Because the wilcoxon test just informs whether there is significant change between two metrics. However the LME models do a prediction of change.

However, I would like to understand if there is a better way of explaining this, because in the interpretation of LME, I used the words "significant change over time" when maybe I should have been using "predicted change of B=xxx, SE=xxx and p value=xxxx"

Is that right? could you give me advise on this please?

Thank you very much in advance Best regards,

EDIT Comment: I am wondering if I could say that the Wilcoxon test gives me information about the difference between T1 at timepoint 1 and T1 at timepoint 2, whereas the LME, in a multivariate setting, gives me information about the linear relationship of Time with T1 and whether T1 si significantly different depending on the sex or the age. I think this is right but I am not 100% sure


1 Answer 1


A linear change across two timepoints is actually exactly the same as a difference in means between the two timepoints, so there is no need to change that language. You still have found a "significant difference" and it's perfectly fine to describe it as that. However, all that significance means is that your dataset would be fairly unusual if there were truly no difference between time V1 and V2. Which actually isn't that informative! So, for sure, it's a good idea to report the B and its SE too. But, for the reason that they tell us much more than a $p$-value; not for any reason that is peculiar to the mixed model.

(For example: your Wilcoxon provides you with a pseudomedian difference and its confidence interval too, which are prodding at very similar things to the B and SE of the mixed model. Again, they are much more informative than a $p$-value).

Your Controls

Regarding the controls, yes I think you describe it well! You are seeing how each of the three things (Time, sex, age) affect T1. (But, in case I've read your sentence wrong: you are not asking if time-differences are impacted by sex and age; that would require you to add an interaction term to your model).

In most cases, adding those extra controls in the model also would help it to tidy up the T1 variable, so that once sex and age effects have been removed, the model gets a clearer view of the differences across time. Because you have paired samples, I think I'm right in saying, age and sex don't add anything in this sense and if you removed them, the $p$-value on Time should stay about the same. Because you already have a random effect that controls for the average of each individual, and the effect of their age and sex is already captured in that. (This is assuming that ID captures participants and the ages/sexes are the ages/sexes of these same participants).

So Why Do the $p$-Values Differ?

If the controls are indeed not controlling for anything, then the time-component of your mixed model is now roughly equivalent to a paired $t$-test (I suspect a paired $t$-test should give the same p-value as the mixed model). And the differences in $p$-value (compared with Wilcoxon) is really down to the parametric vs non-parametric thing.

Both tests try to answer fairly similar questions in (in my opinion) quite similar but, nonetheless, different ways. That one is significant and the other not is no contradiction, because a non-significant result is not the opposite of a significant one -- it is just one that didn't quite reach a threshold. It is a reflection of the fact that a slightly different angle can colour things slightly differently. And even more, it might be a reflection of the fact that a hard line at $p=.05$ isn't really meaningful. If you compared the two models in terms of B and its CI for the mixed model (you can get that by running intervals on your mixed model) and the pseudomedian difference and its CI (which you can get by setting confint=TRUE in wilcox.test) it's quite possible that the results won't differ that much at all...

(The best approach would, of course, be to use just one approach in the first place, so that we are not tempted to choose the $p$-values that suit us -- a.k.a. $p$-hacking).


By the way, I'm not super familiar with the nlme package, but I think you want the lme function in the nlme library. (You are using the nlme function at present, which is for non-linear models, I think.) Also, I think your random factor needs to be something more like random=~1|ID, but again, not super-familiar, so that might be wrong.

  • $\begingroup$ I don't understand the "expected to get the same p-value" statement. $\endgroup$ May 3, 2021 at 11:55
  • $\begingroup$ @FrankHarrell Clearly my brain is in reverse today. I missed the Wilcoxon bit. Will delete my answer and start again!! $\endgroup$ May 3, 2021 at 12:37
  • 1
    $\begingroup$ (And I should say, thanks for alerting me to that!) $\endgroup$ May 3, 2021 at 20:44
  • $\begingroup$ Thank you so much Phil!! very clear :) - with regards to the model in R, nlme works with gaussian and non gaussian model distributions. The random effects as you had expressed them are correct for lme, but nlme has a different syntaxis. - thanks again!! $\endgroup$
    – Lili
    May 4, 2021 at 8:27
  • $\begingroup$ @Lili you're welcome! But, I'm not sure about your interpretation of nlme: from what I can see it is mainly motivated by modelling nonlinearities, I don't see anything in the help about non-Gaussian error components, and I don't see anything in your call to the function asking for non-Gaussian errors. I would definitely spend a minute checking that it's doing what you want it to do. $\endgroup$ May 4, 2021 at 10:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.