Inference of nonparametric tests as linear models

Question

I found some statements throughout the web that suggest that most common statistical tests can be performed as general(ized) linear models (cf. here). For a Wilcoxon test the author of the referenced source exemplarily suggests to first transform the original data into signed ranks and then to run a linear model. Imo this procedure won't give reliable p values because now the residuals are of course not normal but uniform. However, I wonder whether a generalized linear model could help here (e.g. after transforming the ranks to a range of [0,1] and then apply a GLM with beta link)? If so, which GLM could be applied here? And any idea how this analysis could be performed using R? (The Wilcoxon test was just an example - how could this look like for a Spearman correlation analysis?)

Answer 1

You might be interested in the proportional odds ordinal logistic regression (POOLR) model.

Wilcoxon-Mann-Whitney U is a score test of a POOLR model with one variable indicating group membership, similar to how an equal-variance t-test is a Wald test of a linear regression with one variable indicating group membership.

This is not the same as applying a regression model to the rank-transformed data, but that seems like a hack to use until one learns about the POOLR model.

Answer 2

The author writes the following:

# Built-in
a = wilcox.test(y)
Equivalent linear model
b = lm(signed_rank(y) ~ 1)  # See? Same model as above, just on signed ranks

But that is only the same model because the function wilcox.test is approximating the distribution of the mean of the signed ranks by a Normal distribution.

What happens is that the test estimates the error based on the residuals and computes a t-value and assumes that this is approximately t-distributed.

It uses the following formula:

$$t = \frac{\bar{X}}{\hat{\sigma} / \sqrt{n}}$$

in code you could use:

z = signed_rank(y)                         # compute signed ranks
t = mean(z) / (var(z)^0.5 / length(z)^0.5) # t-value
pt(t,n)*2                                  # p value

The Wilcox signed rank test does not follow a simple parametric distribution and it is difficult to compute. If we would compute the exact p value then we get a different result. The code below demonstrates this.

# data
set.seed(1)
signed_rank = function(x) sign(x) * rank(abs(x))
y = c(rnorm(15), exp(rnorm(15)), runif(20, min = -3, max = 0))  # Almost zero mean, not normal
Built-in
wilcox.test(y)                # p-value = 0.7282
wilcox.test(y, exact = TRUE)  # p-value = 0.7304

In the article that you refer to several models like anova and t-test can indeed be interpreted as linear models.

However, the non-parametric tests are only indirectly related to linear models. The relation happens because the distribution of the test statistic is approximately normal distributed under the null hypothesis (when the sample size is large).

---

Note: I tried to look at smaller sample sizes and found that the wilcox.test is actually estimating the standard deviation differently. It uses the function

        NTIES <- table(r)
        z <- STATISTIC - n * (n + 1)/4
        SIGMA <- sqrt(n * (n + 1) * (2 * n + 1) / 24
                      - sum(NTIES^3 - NTIES) / 48)

and it bases the estimate of the standard deviation purely on the sample size, and not on the data. But for large samples, they approach each other. That we get the same p-values might just be a coincidence and is due to the fact that all these approximations approach the same normal distribution.