# Interpreting residual plots for heteroscedasticity

I am trying to measure the effect of different antibiotics on the growth of E. coli and S.aureus. Since my data for both bacteria is not normal I have chosen to do a Kruskal–Wallis test. But one of the assumptions of the Kruskal–Wallis test is there should be no heteroscedasticity. To check for heteroscedasticity I made some residual plots but I am having a hard time interpreting them.

For my E. coli data, I made the following plot:

Does this show that my E. coli data is heteroscedastic?

Also, when I made a residual plot for my S.aureus data I got the following result. This one seems to suggest that the S.aureus data is not heteroscedastic. Am I correct about this?

Edit: My response variable is the diameter of growth inhibition zones (in mm). For example, the table below shows my data for E.coli.

I did a Welch's anova on my E.coli data but got no p-value. Can someone explain what I am doing wrong?

One-way analysis of means (not assuming equal variances)

data:  values and ind
F = NaN, num df = 17, denom df = NaN, p-value = NA
$$$$

• What is your response variable and how is it measured? Your plots give one the impression it might be discrete, such as a small count, suggesting that perhaps using an appropriate model (such as a GLM) is called for.
– whuber
Dec 30, 2019 at 22:18
• I made an edit and mentioned the response variable. Dec 31, 2019 at 14:47
• 1. I suspect the Welch's anova is failing because there is no variance in some groups. Yes... The function in R doesn't like that. (Although the function for Welch's t test works in this condition). If you fit a traditional anova model, the residuals aren't particularly normal or homoscedastic, so, yes, a nonparametric approach or generalized linear model (glm) may work. Dec 31, 2019 at 19:03
• 2. While your measured variable is in theory continuous (diameter in mm), the measurements have little precision (that is, many values have only one significant figure). I'm not sure what approach I'd take if I were doing it myself, but I'll add one potential solution to my answer. Dec 31, 2019 at 19:04

Actually, the Kruskal-Wallis test does not have an assumption of homoscedasticity. The assumptions for the test given by Conover † are the following.

1. All samples are random samples from their respective populations.

2. In addition to independence within each sample, there is mutual independence among the various samples.

3. The measurement scale is at least ordinal.

4. Either the k population distribution functions are identical, or else some of the populations tend to yield larger values than other populations do.

† Conover, Practical Nonparametric Statistics, 3rd.

Edit: Since OP has provided some data, I can add more to this answer.

As discussed some in the comments, the data that OP has is in nature continuous (diameter in mm), but has low precision, so acts more like a discrete variable. In this case, I'd probably be tempted to use a model for discrete values. (Perhaps negative binomial regression?) But because the data are unusual in that some groups have no variance, and some have considerably more, I might just feel better with the nonparametric approach.

The following conducts a Kruskal-Wallis test, and then a Dunn (1964) test for post-hoc. I made no p value corrections here for multiple tests, although you may want to.

Also, at the end, there's a function that reports the pairwise effect size statistics, including Vargha and Delaney's A. This effect size statistic is related to the probability of an observation in one group being larger than an observation in the other group. For example, comparing C30 to the control, vda=0, and in fact every observation in C30 is greater than every observation in the control, but these two treatments are not significantly different in the Dunn test. It is up to you how you use this information.

if(!require(tidyr)){install.packages("tidyr")}
if(!require(FSA)){install.packages("FSA")}
if(!require(PMCMRplus)){install.packages("PMCMRplus")}
if(!require(rcompanion)){install.packages("rcompanion")}

### Import data and translate it to long format

Control KF30 P10  APR15 x510 TE30   AK30  C30   TOB10  CN10   N3     efiran Bleach Soap ETOH70 ETOH96 Detergent Mouthwash
6       23  6    15    11    11     18    7     17     16     17     16     6      6     6     9     6          15
6       10  6    15    12    22     18    7     16     16     17     15     6      6     6     7     6          6
6       11  6    15    12    22     18    8     15     15     17     14     6      6     6     6     6          6
6       10  6    17    12    16     21    8     18     20     19     19     6      6     6     6     6          6
6       12  6    15    11    23     17    8     16     16     15     16     7      6     7     7     6          6
6       12  6    20    10    22     18    8     16     16     16     19     6      6     6     7     6          6
6       11  6    16    12    18     19    8     15     12     18     14     12     6     6     7     6          6
6       12  6    16    13    21     20    9     17     15     17     13     12     7     6     7     6          6
6       12  6    15    12    16     20    8     16     15     16     30     10     6     6     10    6          6
")

library(tidyr)

DataLong <- gather(Data, Treatment, Diameter, Control:Mouthwash, factor_key=TRUE)

str(DataLong)

### 'data.frame':    162 obs. of  2 variables:
###  $$Treatment: Factor w/ 18 levels "Control","KF30",..: 1 1 1 1 1 1 1 1 1 2 ... ###$$ Diameter : int  6 6 6 6 6 6 6 6 6 23 ...

### Plot and summarize data

library(FSA)

Summarize(Diameter ~ Treatment, data=DataLong)

###    Treatment n      mean        sd min Q1 median Q3 max
### 1    Control 9  6.000000 0.0000000   6  6      6  6   6
### 2       KF30 9 12.555556 4.0034707  10 11     12 12  23
### 3        P10 9  6.000000 0.0000000   6  6      6  6   6
### 4      APR15 9 16.000000 1.6583124  15 15     15 16  20
### 5       X510 9 11.666667 0.8660254  10 11     12 12  13
### 6       TE30 9 19.000000 4.0311289  11 16     21 22  23
### 7       AK30 9 18.777778 1.3017083  17 18     18 20  21
### 8        C30 9  7.888889 0.6009252   7  8      8  8   9
### 9      TOB10 9 16.222222 0.9718253  15 16     16 17  18
### 10      CN10 9 15.666667 2.0615528  12 15     16 16  20
### 11        N3 9 16.888889 1.1666667  15 16     17 17  19
### 12    efiran 9 17.333333 5.1961524  13 14     16 19  30
### 13    Bleach 9  7.888889 2.6666667   6  6      6 10  12
### 14      Soap 9  6.111111 0.3333333   6  6      6  6   7
### 15    ETOH70 9  6.111111 0.3333333   6  6      6  6   7
### 16    ETOH96 9  7.333333 1.3228757   6  7      7  7  10
### 17 Detergent 9  6.000000 0.0000000   6  6      6  6   6
### 18 Mouthwash 9  7.000000 3.0000000   6  6      6  6  15

plot(Diameter ~ Treatment, data=DataLong)

### (Plot)

### Kruskal-Wallis and Dunn tests

kruskal.test(Diameter ~ Treatment, data=DataLong)

### Kruskal-Wallis rank sum test
###
### Kruskal-Wallis chi-squared = 142.47, df = 17, p-value < 2.2e-16

### Order groups by median

Sum = Summarize(Diameter ~ Treatment, data=DataLong)

Sum2 = Sum[order(Sum$median),] DataLong$$Treatment = factor(DataLong$$Treatment, levels=Sum2$Treatment)

plot(Diameter ~ Treatment, data=DataLong)

### (Plot with groups ordered by medians)

library(PMCMRplus)

Dunn = kwAllPairsDunnTest(Diameter ~ Treatment, data = DataLong, p.adjust.method = "none")

Dunn

### (Large output table)

summaryGroup(Dunn)

###           median Q25 Q75 n Sig. group
### Control        6   6   6 9          a
### P10            6   6   6 9          a
### Bleach         6   6  10 9         ab
### Soap           6   6   6 9          a
### ETOH70         6   6   6 9          a
### Detergent      6   6   6 9          a
### Mouthwash      6   6   6 9          a
### ETOH96         7   7   7 9        abc
### C30            8   8   8 9        abc
### KF30          12  11  12 9         cd
### x510          12  11  12 9        bcd
### APR15         15  15  16 9         de
### TOB10         16  16  17 9         de
### CN10          16  15  16 9         de
### efiran        16  14  19 9         de
### N3            17  16  17 9         de
### AK30          18  18  20 9          e
### TE30          21  16  22 9          e

### Report Vargha and Delaney's A for pairwise comparisons

multiVDA(Diameter ~ Treatment, data=DataLong)\$pairs

### (Large table)

### Note with vda, 0.50 = no effect; 0 or 1 = stochastic dominance


Seems like in both of your plots, suffer from heteroscedasticity, although the latter is in better shape. Well the kruskal wallis test indeed does not assume homoscedasticity, however still is better not to use it, because it assumes that the distributions between the populations are identical, so having heteroscedasticity doesn't help. You should first try Bartlett's test in order to test your samples for heteroscedasticity and then Welch's anova as an alternative. Try studying here as a resource material. Some really good papers should be:

The last source I provided is about how welch's test is type error I robust which is the main problem heteroscedasticity causes.

• because it assumes that the distributions between the populations are identical`... What could this mean? If a statistical test assumed the populations were identical, what would a statistical test be testing for? Dec 31, 2019 at 1:02
• The whole idea behind kruskal wallis is that if the distributions are identically shaped between the populations, then their medians are equal - that is the null Hypothesis (apart from the other assumptions). Having equally distributed samples does not mean they have the same median. Dec 31, 2019 at 1:16
• I tried to do a Welch's anova but I am getting no p-value. Can you please check my edit above? Dec 31, 2019 at 15:34
• Kruskal-Wallis doesn't assume the distributions in the groups have the same shape and spread. If you make that assumption, then the test becomes a test of medians. But that's much less fun than using as a test of stochastic dominance. Dec 31, 2019 at 19:07