Validity of using a combination of Wilcoxon tests and Spearman rho as an alternative to GLM?

Question

I recently read a paper where the authors used either (a) a generalised linear model or (b) a combination of Wilcoxon rank sum tests and Spearman’s rank correlation, depending on whether the data followed a normal distribution or not. They also used a Bonferroni correction when applicable, to account for multiple comparisons.

The responses were two continuous variables, each analysed separately. They had 3 categorical variables and 1 continuous variable as independent variables. (a) For the first response variable, data followed a normal distribution, so a GLM was fitted assuming a normal distribution with identity link. Inference was (apparently) made by significance of model terms, i.e. whether or not the independent variable was significant in the model or not. (b) For the second response variable, the data did not follow a normal distribution, so a Wilcoxon rank sum tests and a Spearman’s rank correlation test were used instead.

My questions are:

1) Are these two approaches appropriate and equally valid alternatives, depending on (non)normality?

2) Even in case of non-normality, wouldn't a different GLM type be better than a combination of Wilcoxon rank sum tests and a Spearman’s rank correlation? (Not sure which type of GLM would be appropriate though, as the response is still continuous...gamma perhaps?)

3) In a GLM framework, aren't we worried about non-normality of residuals, rather than raw data (which is not always the same thing)?

4) Doesn't a GLM (generalized linear model) with "normal distribution and identity link" simply mean a general linear model (i.e. an ordinary linear model)?

5) Is this approach of using p-values from a full model (including all variables) appropriate, or would a variable selection (using AIC or similar) be more appropriate?

Answer 1

To clarify terminology, a GLM with Gaussian distribution and identity link is identical to the general linear model or linear regression. So I will simply refer to it as linear regression going forward.

Additionally, the normality assumption in regression is on the error term. Since the residuals are estimated errors, we verify this assumption using the residuals or variants of the residuals, like studentized deleted residuals. Claiming the distributional assumption applies to the error term works in linear regression. However, for other GLMs like logistic regression, we do not assume the error term is binomial. A more general assumption for linear regression is that the data are normally distributed with a mean that depends on the predictors and a variance that is constant. Since the mean of the distribution depends on the predictors, to verify normality, we extract these means, hence giving all of the data the same center/location. This is what residuals are and then we can plot these residuals to check normality. Beyond linear regression, such empirical verification becomes more complicated.

It happens that the normality assumption in linear regression is one of the assumptions we do not need to care much about. The least squares coefficients do not depend on normality. However, classical inference using t-tests relies on the normality assumption. But from central limit theorem, we know this assumption does not matter too much once the sample sizes are large enough.

Rather than moving away from linear regression because of non-normality, I might move to another GLM if I could hypothesize a theoretical reason for why another distribution would be more appropriate. For example, if my data were count so non-negative and they were largely low counts with a few high values, I might start to consider a Poisson distribution. If the data were binary, then Binomial distribution is reasonable. The normal distribution is a reasonable choice if all we are willing to assume is that the data have finite variance and their range is $(-\infty, \infty)$. Usually, we can make more assumptions about the data than these two assumptions.

Additionally, the linear regression can handle multiple predictors predicting a single outcome. Wilcoxon/Spearman are for single predictor situations, so they are not comparable to linear regression in this way. They are directly comparable to specialized linear regressions like the Pearson correlation and the t test. However, they ask different questions of the data. Wilcoxon tests are stochastic dominance tests comparing two groups, while Spearman test is a measure of monotonic relation (like an ordinal linear correlation). So even when they are similar to linear regression, they ask different questions of the data.

Finally, with regard to your final question, modeling is hard. Model selection is a difficult problem. A reasonable approach is to attempt different models then observe their varying implications. If their implications are similar, that makes life easy. If they differ, that can be interesting for research purposes. With Bayesian modeling, we have a few more options in that we do not simply have to select a model, we can combine different models to account for our uncertainty in our final model.

Edit

An error on my part from the first paragraph: general linear model also includes MANOVA; this differentiates general linear model from the standard linear regression. All other members of the general linear model family are linear regression models or generalized linear models with gaussian distribution and identity link function.

Answer 2

The model that is a generalization of Wilcoxon, Kruskal-Wallis, and, to some extent, Spearman's $\rho$, is the proportional odds ordinal logistic semiparametric model. For short we call it the proportional odds (PO) model. The PO model handles the full generality of the linear model plus more, because it is invariant to how $Y$ is transformed. This circumvents the need to guess about whether residuals have a normal distribution. The PO model is robust with regard to $Y$ and is 0.95 as efficient as the linear model if normality holds. Otherwise it can be more efficient than the linear model.

Software needs to make it fast to handle the case where there is a large number $k$ of distinct $Y$ values, as the number of intercepts in the PO model is $k-1$. The R rms package orm function is fast for $k$ up to about 6000.

See my detailed case study in my RMS course notes for modeling continuous $Y$ using ordinal regression.

Once you are comfortable with ordinal response models, you no longer even need the special cases such as Wilcoxon. And ironically the PO model handles extreme ties better than the Wilcoxon test.

Answer 3

I think Jim's answer addresses all your points adequately, but I'll try to boil it down to the per-question points.

1) Are these two approaches appropriate and equally valid alternatives, depending on (non)normality?

Definitely not. Regression, in general, seeks to describe the relationship between typical values of a random variable given known information. In their comment, whuber describes the (very different) purpose of the other two tests.

2) Even in case of non-normality, wouldn't a different GLM type be better than a combination of Wilcoxon rank sum tests and a Spearman’s rank correlation? (Not sure which type of GLM would be appropriate though, as the response is still continuous...gamma perhaps?)

Definitely so. We would need more information about the variable to provide specific guidance. Residual plots are a good way to convey the presumed lack of normality.

3) In a GLM framework, aren't we worried about non-normality of residuals, rather than raw data (which is not always the same thing)?

Correct. This is almost never the same thing. Residual plots are a good tool to explore this in linear regression. As Jim mentions, deviations from normality are not in and of themselves a great concern when you have a moderately sized sample and care only about the regression coefficients or mean predicted values.

On the other hand, a more important condition that residual plots also allow you to assess is a mean value of zero for the residuals. If the residuals stray too far at times, it suggests the model may be misspecified.

4) Doesn't a GLM (generalized linear model) with "normal distribution and identity link" simply mean a general linear model (i.e. an ordinary linear model)?

Not quite. A "general linear model" refers to a, well, general version of linear models. This generality also allows for multiple outcomes to be analyzed simultaneously. Generalized linear models only handle univariate outcomes. The normal-identity GLM could be described as "multiple linear regression". That makes it both a particular case of a GLM and of the general linear model.

5) Is this approach of using p-values from a full model (including all variables) appropriate, or would a variable selection (using AIC or similar) be more appropriate?

A p-value only allows you to say whether a statistically significant difference was observed, not that it wasn't, but it fails at communicating anything beyond that. I suggest using the confidence intervals to describe the findings provided by the model, as they show the range of plausible values for the estimated parameters. Even if a variable is significant, a wide confidence interval can suggest that little useful information about that particular relationship was gleaned. Conversely, a non-significant variable with a narrow confidence interval suggests that the presumed relationship is either absent or small in magnitude.

Using AIC for variable selection addresses a very different question. Selecting the AIC-best model will, on average, give you the model that makes the best predictions while also being as small as possible.