# Wilcoxon signed-rank test with an additional covariate

## Problem formulation

I have 5 strains on which I administered treatment and repeated the experiment two times independently. So my data looks something like this:

strain <- factor(1:5)
ctrl <- c(1, 1, 1, 1, 1)
trt1 <- c(60.82812657, 6.175973947, 7.395597282, 203.6573398, 24.42014734)
trt2 <- c(758.3219468, 26.84658479, 43.21119828, 103.2501452, 9.557721081)


with the output

strain ctrl       trt1       trt2
1    1  60.828127 758.321947
2    1   6.175974  26.846585
3    1   7.395597  43.211198
4    1 203.657340 103.250145
5    1  24.420147   9.557721


The values are normalised to the control (that's why the controls are always 1). I want put statistical significance to the question: "Does the treatment increase the measurement values?"

Problem: What is the most appropriate test to answer this question, given this specific structure of the experiment?

## What I have tried so far

### Ignoring the strain structure

The values don't seem to be normally distributed, so I thought a non-parametric test would be more appropriate. If I perform a Wilcoxon t-ttest individually for each treatment group

wilcox.test(ctrl, trt1)
wilcox.test(ctrl, trt2)


I get a p-value of 0.007495 each (together with an warning message: “cannot compute exact p-value with ties”).
Pooling the two treatments together increases the significance level as expected, i.e.

wilcox.test(ctrl, c(trt1, trt2))


yields a p-value of 0.002245 (together with a warning message: “cannot compute exact p-value with ties”).

### Including the strain structure

What came as a surprise to me is that the paired version of this test is not statistically significant anymore, e.g. the following code

wilcox.test(ctrl, trt1, paired = TRUE)
wilcox.test(ctrl, trt2, paired = TRUE)


yields a p-value of 0.0625 each (without a warning message). Now I would like to pool the treatment groups together (similar to the case above) but still somehow retain the strain-structure of the experiment. But this is where I am stuck right now.

Anyone knows if I am on the right track, or how to include the categorical variable strain as some sort of covariate into a non-parametric statistical model?

• If you had normality, what would you do. Depending on your response, there could be a fairly straightforward answer, as proportional odds models generalize Wilcoxon and Kruskal-Wallis analogous to how least squares linear regression generalized the usual t-test and ANOVA F-test.
– Dave
Jan 20, 2023 at 15:43
• A paired test is inappropriate in this case, as you have already normalized everything to the control. Furthermore, since all the control values are the same, you will have a problem with ties, greatly weakening the power of the test. Jan 20, 2023 at 16:18
• @Dave I am not sure myself what I would have done with normality. Probably something with a linear model à la lm(value ~ treatment + strain, data = data_df) and then looked at the p-value for the treatment coefficient. Jan 20, 2023 at 16:22
• Proportional odds modeling is analogous, and you might want to look into the topic.
– Dave
Jan 20, 2023 at 16:46
• Because of the normalization of controls, this is not a 2-sample problem, but rather a 1-sample problem. You want to test whether the mean of the treated units in each strain differs from 1. You don't need the control units to do this. But testing whether the mean of 5 numbers differs from a given number is not straightforward because the sample is so small, way smaller than most statistical analyses require. You have no idea whether the values are normally distributed when you only have 5 values. Why do inference at all? How could you expect 5 samples to generalize to anything meaningful?
– Noah
Jan 20, 2023 at 18:01

As mentioned in the comments, since your data is normalized to the control, you essentially have a one-sample problem, comparing the treatment data to the control value of 1.

You could conduct analyses for each treatment separately, with, for example, a one-sample t-test (if appropriate), one-sample Wilcoxon test, or one-sample sign test. As mentioned in the comments with five observations per treatment, these analyses are not likely to be very helpful.

But as you mention, you'd like to consider the two treatments together with the acknowledgement that they are paired by Strain.

For this particular case, here's what I might do. This isn't necessarily an authoritative answer, but it makes sense to me given the limitations of the study.

1. Log transform the observations. This may make sense depending on how you normalized the values. But it also makes the observations relatively symmetric, makes the distribution of values easier to examine.

2. Construct a model whereby the structure of the observations being paired by Strain.

3. From there, construct confidence intervals for the estimated marginal means of the treatments. If these confidence intervals don't include 0 (log(1)=0), when you have evidence that the e.m. for the treatment is different from the control.

3b) As it turns out, for these data, it won't matter much if you ignore Strain. You could simply do the analysis by just calculating the mean and confidence intervals for the log of the observations.

3c) I decided to use a mixed-effects model with Strain as a random effect. Since there are only five levels of this variable, it may not be the most justified approach to treat this as a random effect.

These are the results of this approach. All values are log-transformed. The upshot: The confidence interval for the e.m. means do not include 0, so this suggests the values for the treatments differ from the control. But the confidence intervals for the treatments overlap quite a bit, suggesting that the two treatments aren't different.

 Treatment  emmean    SE   df lower.CL upper.CL
Treatment1   3.29 0.698 6.35     1.60     4.97
Treatment2   4.12 0.698 6.35     2.43     5.80


As mentioned, just ignoring Strain, and calculating the confidence intervals for the means of the logged values is similar. Here, these are done by bootstrap.

   Treatment n Mean Conf.level Percentile.lower Percentile.upper
1 Treatment1 5 3.29       0.95             2.17             4.47
2 Treatment2 5 4.12       0.95             2.97             5.49


The R code I used:

Data = read.table(header=TRUE, stringsAsFactors=TRUE, text="
strain ctrl       trt1       trt2
1    1  60.828127 758.321947
2    1   6.175974  26.846585
3    1   7.395597  43.211198
4    1 203.657340 103.250145
5    1  24.420147   9.557721
")

Y = c(Data$$trt1, Data$$trt2)
Treatment = factor(c(rep("Treatment1", 5), rep("Treatment2", 5)))
Strain   = factor(c(1:5, 1:5))
Data2 = data.frame(Strain, Treatment, Y, LogY = log(Y))

plot(log(Y)~Treatment, data=Data2)

library(lme4)

library(lmerTest)

model = lmer(LogY ~ Treatment + (1|Strain), data=Data2)

hist(residuals(model))

plot(residuals(model) ~ predict(model))

anova(model)

library(emmeans)

emmeans(model, ~ Treatment)

library(rcompanion)

groupwiseMean(LogY ~ Treatment, data=Data2, percentile=TRUE)