# Can we perform a MANCOVA for paired data?

I have data for a cross-over randomized controlled trial. I wonder if I can use MANCOVA to compare my outcomes (they are continuous) while having paired data. If I can't, what would be an appropriate test here? I have some covariates that I want to account for. I already did some Wilcoxon tests, but it didn't account for covariates.

The design of the study is as follows:

There are two groups: Group is given treatment (a drug) and then control (Sequence AP). Another group is given control then treatment (the same drug) (PA sequence). Between the treatment and control, there is a washout period.

I have four DVs that are continuous. Each variable is measured two times. For example, the AP sequence had the variables measured one time after the intervention (t0) and another time after the place (t1). Same thing for the PA sequence. They had the placebo, then had the variables measured (t0), then a washout period, then the intervention, and then the DVs were measured again (t1). I have covariates that are measured at baseline (before the study starts). So think of my data as a dataframe with the following columns:

1- Sequence (AP vs PA) 2- t0 (DV1) 3- t1 (DV1) 4- Covariate 1 5- Covariate 2 and so on

I would like to:

1- Compare between t0 and t1 for each sequence while controlling for covariates. I started with Wilcoxon test to compare t0 and t1 considering that this is paired data and some of the DVs are not normally distributed. However, I know that Wilcoxon doesn't account for my covariates. I thought about using MANCOVA (but am not sure if it is the right one because some of my DVs are not normally distributed and my data is paired data). Ignoring the normality thing, can we use MANCOVA for this scenario? In case of Mancova, I will have to convert my wide data to long data (i.e., instead of having each DV presented as column t0 once and t1 once, I will put them on top of each other). My DVs would be my four continuous DVs, then we have an independent categorical variable (drug vs placebo) and continuous covariates. I can do that for each sequence separately, then combine them.

Btw, Wilcoxon didn't find any significant difference.

2- If MANCOVA is not valid, what would be a simple but robust solution?

The correlation method below would work? I assume if the difference between t0 and t1 is not correlated with my covariates, then it is unlikely that they might be confounding variables.

Update: I am thinking of using correlation between the difference in my DVs (t1-t0) with my covariates to create a correlation matrix. If they are not significant, then it is safe to say that these covariates are not confounding my results. Does that make sense?

• It's easy to get hung up on terminology. There certainly is some linear model that can handle your data properly. Whether that model is technically called MANCOVA doesn't really matter. Please edit the question to say more about your study, its design, and your data; you should then be able to get some useful advice.
– EdM
Nov 20, 2023 at 14:57
• Thanks for the update. Please edit the question again to say more about the 4 DVs and the design. Are the DVs continuous? Do you have exactly one measurement of each DV per individual for the treatment condition and another one for the placebo condition? Is there more than one treatment being evaluated, with emphasis on treatment-placebo comparisons among the treatments? Was the order of treatment/placebo randomized among individuals?
– EdM
Nov 20, 2023 at 18:27
• Thank you so much for following up with me. I have updated my question and provided a comprehensive overview of the study. @EdM Nov 20, 2023 at 18:51

In this specific scenario, you can extend a paired t-test approach to a multivariate (in the sense of multiple outcomes) analysis of covariance (MANCOVA). That's just a name for a particular type of multivariate multiple-regression model.

The trick is to use, for each individual, the treatment - control difference as the outcome variable for each of the 4 dependent variables (DV) that you measured after each treatment or control period. That's what paired t-tests would do, anyway, comparing the differences for each DV against a value of 0.

Include as a covariate a binary indicator of the presentation order, for example taking the value of 0 for control first or 1 for treatment first. That will evaluate whether the treatment - control difference is affected by the order of presentation.

Then include whatever additional covariates that you would like to control for, and perform MANCOVA. The regression coefficient estmates will be the same as you would get for doing separate regressions for each of the DV differences, but the covariances of the estimates will take into account correlations among outcomes.

You will have to be a bit careful in evaluating the results. If you use standard treatment coding for covariates and handle the treatment order as above, the intercepts will be the treatment - control differences for the control-first group, at reference or 0 values for all of the covariates.

The coefficients for the binary indicator of presentation order will be the differences of those differences between the treatment-first and control-first groups. One hopes that those coefficients for presentation order will be close to 0, supporting a lack of treatment-order effect. In that case you should be able to get a combined estimate of each of the 4 treatment - control differences by pooling the results (both coefficient estimates and their covariances) from the two presentation orders.

The coefficients for the other covariates will indicate how much each affects each of the treatment - control differences.

That will be a better approach than looking at pairwise correlations of differences against each of the covariates. Unlike the pairwise correlations, the multiple regression underlying MANCOVA takes all of the covariates into account simultaneously.

You don't need to worry about normal distributions of each of the DV at all in the above approach, or even about normal distributions of the treatment - control differences. What might be of concern is a reasonably normal distribution of residuals between observations and the predictions from the MANCOVA model, which take the covariates into account. Even that isn't absolutely required for a valid model, but a reasonably normal distribution of residuals is very reassuring.

If you ever need to include multiple covariates or continuous predictors in a Wilcoxon-type analysis, consider ordinal regression. Frank Harrell provides resources for ordinal regression models that only evaluate the rank-ordering of outcomes. Wilcoxon tests can be considered simple, special cases of more general ordinal regression models.

A mixed model is one way to analyze your data while taking the paired observations among individuals and covariates into account, but lmer() works only for one DV at a time. It's not a true "multivariate" model, which would assess correlations among the multiple 4 DV. There are ways to do true multivariate mixed models, as described on this page, but they require extra care. That said, a set fo 4 separate mixed models (one for each DV) might accomplish your goals.

The way that you attempted to fit a multivariate model

man_cov_model <- manova(cbind(DV1 + DV2) ~ independent_variable + covariate , data = data)


did NOT take the pairing of observations within individuals into account. I suspect that's why you didn't find "significant" results with it. With paired treatment and control values for each DV and each individual in your case, you could get around that by pre-calculating the treatment-control differences for each DV and individual. Use those differences at the outcomes in the multivariate model. Then you no longer need independent_variable as a predictor in your model; the intercept is the estimated treatment-control difference after accounting for covariates. Re-read the answer starting with the second paragraph.

If you take that true multivariate approach, I'd recommend avoiding the manova() function and instead following the approach suggested by Fox and Weisberg in the appendix to their book on regression models. Use lm() to get the model, then analyze results with the Anova() function (capital "A") in the R car package. If you have a perfectly balanced design manova() might work OK, but otherwise the results from its Type I tests can be misleading.

Finally, your mixed models and your attempt at a multivariate model don't seem to account for the order of presentation (treatment versus control first). If you have complete balance in the order of presentation that's probably OK, but otherwise you should include the order of presentation as a covariate.

• THANK YOU SO MUCH. Final question: Would you consider linear mixed models as an appropriate approach? model <- lmer(DV1~ independent_variable + covariate + (1 | ID), data = data) In this case, I analyzed each sequence separately. Instead of having t0 and t1 as separate columns, they were stacked on top of each other in one variable (DV1). The "independent_variable" indicates treatment or control, and then I have my covariates. I am getting same results as my wilcoxon test, which is interesting. Nov 21, 2023 at 16:57
• Also, when I use mancova (but using the same data as in the lmm in my comment above, man_cov_model <- manova(cbind(DV1 + DV2) ~ independent_variable + covariate , data = data) I get insignificant results, although DV1 and DV2 were significant in the wilcoxon and lmm - any explaination? Nov 21, 2023 at 17:17
• @yusefsoliman I think that manova() uses Type I sums of squares, which can be misleading. I'd recommend following the advice of Fox and Weisberg in this appendix. Do lm() to get the multivariate model and then use the Anova() function (capital "A") and others in the R car package for post-model evaluation.
– EdM
Nov 21, 2023 at 22:54
• Thank you. Final question and sorry for my many questions. Would you consider the linear mixed model approach correct? Because I told you, it gave me similar results as in the Wilcoxon test, but this time accounting for covariates. Also, I read that it accounts for paired findings and it is commonly used in the design of cross-over studies. model <- lmer(DV1~ independent_variable + covariate + (1 | ID), data = data) I explained these in details in the comment above. Nov 21, 2023 at 23:44