Can I use multiple regression on a ranked response variable as a significance test for multiple covariates?

This blog post illustrates the relationship between inference tests on groups (t-test, ANOVA, etc.) and equivalent linear models. It also claims that for reasonable sample size, regression of a ranked variable approaches the nonparametric versions of these tests. The author links to some simulations.

For example, the author claims that for a non-normal response variable, and N > 11,

lm(rank(y) ~ X1 + X2 + X3 + ...)

would be roughly equivalent to the Kruskall-Wallis test. I was under the impression that KW could only handle two groups.

This would be fantastic because I am working with highly skewed genomic data but have multiple confounding demographic variables. For example, we are testing the significance of a response variable to disease state. However, other variables such as Age and Gender not only correlate with disease state, they also independently correlate with the response variable.

Would a journal accept this approach? Are there some references I could back it up with?

Kruskal-Wallis takes N > 2 independent samples. With two groups, it reduces to the Mann-Whitney U. So the blog is accurate enough.

To my knowledge, there is no analytical non-parametric solution when you begin adding covariates. In this case, I would consider:

• Ordinal logistic regression. It's an extension of non-parametric models to more predictor variables- and structures. However, it also requires some expertise to understand the model. See Frequentist version here though I'd always go Bayesian and use something like brm::brms.
• Quantile regression. Simpler and easier to interpret.
• Perhaps there's a meaningful transform of your dependent variable that makes it accessible to parametric (though making conclusions about raw data from analyses on transformed data often come with its own problems, unfortunately).
• For personal use, I would certainly not be afraid to do what you propose, but stakes are higher when it's research (as they should).

I'm afraid I don't know of academic literature backing up the parametric-on-ranks approach other than at the most general level (see this blog for some literature on this in the context of analysis of Likert-type response formats). I often find answers to the question "would a journal accept this approach" to be more about politics and habits than substance.

Your question has many different aspects:

Kruskall Wallis test is not for multiple independent variables

Note that the Kruskall Wallis H-test is a one-way test (it relates to a single independent variable $$X_1$$, albeit with multiple levels). This is already a crucial differentiation with the situation that you are looking at (multiple regressors $$X_1, X_2, X_3, ...$$).

The expression used in the blog $$rank(y) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3x_3 + ..$$ is not a linear model that corresponds to one-way anova (By this I mean 'not in general'. It only turns out that way in the specific case that the $$x_i$$ are dummy coded variables for a single categorical variable).

The numerical example that is used in the blog is a one-way ANOVA. Do not confuse the formula expression with a linear model that has multiple variables.

Comparing Kruskall Wallis test with linear model and ANOVA

I believe that the relationship/equivalence between KW and linear regression is, I guess, more made with the purpose of education/intuition in mind.

However, indeed, the test statistic in the KW H-test is closely related to the ratio of variances, just like the statistic in the F-test, as used in the ANOVA. Namely: you partition the sum of squared residuals into the total sum of squares (TSS), the explained sum of squares (ESS), and the residual sum of squares (RSS) and take the ratio for two of those as the statistic. With $$N$$ the total number of observations and $$G$$ the number of groups you have:

$$\begin{array}{} H &=& \frac{ESS}{TSS/(N-1)} & \quad \text{which is the statistic used in KW test}\\ F &=& \frac{ESS/(G-1)}{RSS/(N-G)} & \quad \text{which is the statistic used in ANOVA} \end{array}$$

and with $$TSS=ESS+RSS$$ you can relate these two as: $$(N-1) H^{-1} = 1 + \frac{G-1}{N-G} F^{-1}$$

But, while the two are related in the concept (a ratio of the squared residuals), they differ in two points:

• The Kruskall Wallis test assumes that the residuals are distributed according to a uniform distribution (ranks are not normal distributed). The F-test/ANOVA assumes that the residuals are normal distributed. The uniform distribution has shorter tails than the normal distribution and the occurrence of large values for ESS will be less often than one expects based on the variance.
• The Kruskall Wallis test and the ANOVA F-test use different types of ratios. Namely ESS/TSS versus ESS/RSS. The distribution of the sum of squares, ESS and RSS will approach a Chi-squared distribution, but the TSS will be constant (if each rank occurs exactly once and there are no ties then TSS is constant). So the ratio ESS/TSS (a Chi-squared distributed variable divided by a constant) will be approximately Chi-squared distributed and the ratio ESS/RSS (a ratio of two chi-squared distributed variables) will be F-distributed.

These differences mean that the ANOVA F-test makes two "errors" when being used as an approximation:

1. The ANOVA/F-test assumes a different distribution of the residuals, independent normal distributed residuals instead of (somewhat correlated) uniform distributed residuals. This effectively assumes longer tails, ie. higher values are expected to be more probable.
2. The ANOVA/F-test assumes that the denominator RSS is an independently distributed from ESS. This effectively assumes higher values are more improbable (in the rank test where RSS = TSS-ESS and TSS is constant, a high value of ESS is related to a low value of RSS making the ratio ESS/RSS more often high than one would expect when ESS and RSS are independent).

Interestingly these two "errors" sort of cancel-out each other and this makes that, when we use rank as dependent variable, then it actually works better with the linear model (and the ANOVA/F-test) than the Kruskall Wallis test. Or at least... that is when the Kruskall Wallis test uses the Chi-squared distribution (the R function kruskal.test does this) as approximation for the distribution of the H statistic (instead of using tables with exact values).

See the image with simulations below.

• The image compares simulations of distributions for the ratio ESS/RSS (which is used in the F-test) and ESS/TSS (which is used in the Kruskall Wallis test)
• Added to the image are lines that express approximations of these distributions based on a F-distribution (which is used in the F-test) and a Chi-squared distribution (which is used in approximations for the Kruskall Wallis test)
• Note the following:
1. You can see that the points ESS/TSS have lower values than the approximation with the chi-squared distribution (more on the left side of the graph). This is because the uniform distribution has shorter tails than the normal distribution.
2. You can see that while the points ESS/TSS are lower than the chi-squared distribution approximation the points ESS/RSS are actually higher. This is because of the direct coupling between ESS and RSS (they are fully correlated).
3. Everything becomes much the same when the sample size increases. The F-distribution will approach the Chi-distribution (since the denominator becomes less variable, ie more like a constant). The simulations will approach the Chi-distribution, because the group means, a sum of many uniform distributed variables, will approach a normal distribution.

Thus when comparing the the function lm and anova versus the function kruskal.test (the latter uses the chi-squared approximation) then the anova F-test actually performs much better.

• But note that this is only because there are two errors that happen to cancel each other out. I would consider it doubtful to use the method 'just because it works', while underlying are quite some wrong assumptions (the same idea is why the Chi-squared test often works, it is not exact as the fisher exact test, but also the assumed normal distribution is not correct, ie. in reality there may be more outliers, and as a result in practice these two cancel each other and it works out pretty well) .

• I would be skeptical about the computations from the blog. They are only made for a single situation. You should model a lot of values instead and not just compare the p-values of the two methods in a single case. Also, it is very much unclear why the cutoff value is at N>11 and how this value is being influenced by the number of groups.

Based on such considerations I see not why you would use the linear model in the first place. In the end you would still have to verify the validity of the approximation and the best method will still be to compute the distribution of the statistic. Then it does not matter whether you use the statistic $$H$$ or the statistic $$F$$ since they are basically the same and only transformations of each other and when you know the distribution of the one then you know the distribution of the other.

Linear model with more independent variables

We see that with the one-way ANOVA the use of a F-test instead of a Kruskall Wallis test is not too bad. Especially when you compare it with a Kruskall Wallis test that uses the chi-squared distribution as approximation (When you wonder about acceptance by journals; I guess that the chi-squared distribution approximation, e.g. by using R, is worse than the use of the F-test with a linear model, but probably a lot of articles/research that use it have been accepted).

It might be very well possible that you could extend this to linear models with more variables. The assumption of normality is also not so much of importance for ANOVA to work. A good approach would be, I believe, to use simulations and have some sort of Monte Carlo approach to compute significance. For instance, in the case of the one-way Anova why not use the curves that I computed/simulated above in the graph instead of the approximated value?

Something similar could be done for larger linear models, however the difficulty becomes mostly that your deterministic model needs to be a good representation of the mean of the rank. And because rank is such an artificial variable I imagine that it might be likely that rank is not well modeled by a linear model. The result is that your model may not capture the true group means of the ranks, and consequently your model will have less power (introducing a bad model will not reduce the residuals that much).

The use of rank is very arbitrary.

The advantage is in the fact that it allows to develop universal methods that are independent on the distribution of residuals and can be applied to all cases (since for any sample of variables, no matter how they are distributed, the rank turns anything into a uniform distribution).

However any function of the rank will have the same property. E.g. why not use the square of the rank or the logarithm etc?

You might be better of by actually modeling the distribution that you have in your case. There are ways to use (non)-linear models even when the error distribution is skewed (for instance generalized linear model). There are many other possible options here, for instance also note that the distribution of your outcome variable does not need to be normal distributed or non-skewed.

Anyway, it seems that, while you ask about the transformation of your dependent variable to a rank, your real question is how you should deal with your data (which may have more, and different, solutions), so maybe you should provide more information about that data (possibly in a different question since it would/might dramatically change the question here).

While using the rank-transformation allows your data to be used with standard methods, it is also a very crude method. Instead of selecting a ready-made method and hoping that it will work, it would always be better to use knowledge about your data gathering process as the starting point. What do we already know about the way that the data is created. Can we model it? Then possibly you might find out an alternative transformation that will allow you to model the data more precisely and this gives you more power in recognizing significant effects.

code for image and comparing p-values:

set.seed(1)
nsim <- 10^5

# for storing p-values
pkw <- matrix(rep(0,nsim*3),nsim)
pano <- matrix(rep(0,nsim*3),nsim)

# for storing statistic values
stf <- matrix(rep(0,nsim*3),nsim)
sth <- matrix(rep(0,nsim*3),nsim)

# for storing sum of squares
TSS <- matrix(rep(0,nsim*3),nsim)
ESS <- matrix(rep(0,nsim*3),nsim)

# simulate
for (Ngi in 1:3) {
Ng <- c(5,10,20)[Ngi]  # number per group
for (i in 1:nsim) {
x <- rep(1:Nc,each=Ng)
y <- sample(seq(1,Ng*Nc),Ng*Nc)

# compute statistic (and divide by degrees of freedom in groups)
kwt <- kruskal.test(y,x)
mod <- lm(y~as.factor(x))
ano <- anova(mod)

# compute and store sums of squares
means <- colMeans(matrix(y,Ng))
ym <- rep(means,each=Ng)

TSS[i,Ngi] <- sum((y-mean(y))^2)
ESS[i,Ngi] <- sum((ym-mean(y))^2)

#storage of p-values
pkw[i,Ngi] <- kwt$$p.value pano[i,Ngi] <- ano$$Pr(>F)[1]

# storage of H and f values
sth[i,Ngi] <- kwt$$statistic stf[i,Ngi] <- ano$$F value[1]
}
}

# plotting

layout(matrix(1:3,1))

for (Ngi in 1:3) {
Ng <- c(5,10,20)[Ngi]
plot(-1,-1, xlim = c(0,20), ylim = c(10^-3,2),log="y",type="l",
xlab = expression(chi^2), ylab = "p-value")
lines(chi,1-pchisq(chi,Nc-1),col=1)
lines(chi*(Nc-1),1-pf(chi,Nc-1,Nc*Ng-Nc),col=1,lty=2)