# Kruskal-Wallis and Negative binomial regression do not agree

I am comparing the number of broods made by dung beetles (Brood_Number) across three temperature treatments (Temp_Offset, a 3 level factor (+0 deg, +2 deg and +4 deg)). I have run the following negative binomial model:

    m=glm.nb(Brood_Number~Temp_offset, data=Trial2)


Which gives me:

    Call:
glm.nb(formula = Brood_Number ~ Temp_offset, data = Trial2, init.theta = 7.152501499,
link = log)

Deviance Residuals:
Min       1Q   Median       3Q      Max
-2.6469  -0.9302   0.1293   0.5589   1.4051

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)    3.4247     0.1569  21.825  < 2e-16 ***
Temp_offset2  -0.6138     0.2228  -2.756  0.00586 **
Temp_offset3  -0.2210     0.2172  -1.017  0.30900
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(7.1525) family taken to be 1)

Null deviance: 32.663  on 22  degrees of freedom
Residual deviance: 24.939  on 20  degrees of freedom
AIC: 177.33

Number of Fisher Scoring iterations: 1

Theta:  7.15
Std. Err.:  2.92

2 x log-likelihood:  -169.326
>


Comparing the model with and without the Temp_Offset term improves model fit. I have interpreted this as saying that Temp_Offset significantly affects brood number, but that there is only a significant difference between the base level (+0 deg) and the +2 deg treatment. There is no difference between the base level and the +4 deg treatment.

I then performed a Kruskal-Wallis test:

    kruskal.test(Brood_Number~Temp_offset, data=Trial2)


Which gives me:

            Kruskal-Wallis rank sum test

data:  Brood_Number by Temp_offset
Kruskal-Wallis chi-squared = 5.1272, df = 2, p-value = 0.07703


So the Kruskal-Wallis and the NB regression give different answers. I feel that the NB regression is the one to be trusted because NB specifically models overdispersed count data (which I have). Am I correct in trusting the significant effect of temperature from the negative binomial regression instead of the non-significant Kruskal-Wallis? I can see no advantage in the KW test. (It was suggested that I try it by a reviewer ...)

Raw data are as follows. The response is the number of broods. The independent factor is the temperature offset. Brood numbers for the different temperature treatments are (its a small data set):

Brood numbers:

    +0 deg: 32, 47, 19, 23, 16, 45, 33

+2 deg: 21, 12, 3, 22, 20, 6, 29, 20

+4 deg: 16, 31, 14, 24, 30, 26, 40, 16


Thanks

• No mystery here: Assuming a negative binomial model provides information. The KW test does not use that information. – BruceET Jul 8 '19 at 4:41
• @BruceETCould you please clarify your comment? Do you mean that a NB uses additional info that the KW does not? – JeanDrayton Jul 8 '19 at 5:11
• any chance you can post your original data? It would make for a more complete/satisfying answer ... – Ben Bolker Jul 8 '19 at 6:02
• possible duplicate of stats.stackexchange.com/questions/264646/… ? – Ben Bolker Jul 8 '19 at 6:02
• @BenBolker It did look at that question, but its a little different to mine because that question was comparing a NB with a factor plus covariates, with just a KW (with one factor). My KW and NB both only have one factor. I will post the data :) – JeanDrayton Jul 8 '19 at 6:08

## 1 Answer

### take-home messages

• The negative binomial uses more information from the data, so it's expected that it would be slightly more powerful than the rank-based Kruskal-Wallis test. In general, you'd use K-W if you were concerned that the distributional assumptions of your model were badly violated.
• The difference between $$p=0.021$$ (overall significance of the NB model; see below) and $$p=0.077$$ is not enormous; they do lie on opposite sides of the "magic" $$p=0.05$$ line, but we wouldn't necessarily say the tests disagree strongly about the strength of evidence against the null hypothesis.
• The diagnostic plots for the NB suggest that the negative binomial model is entirely defensible in this case.
• You should probably focus on the overall test of significance (see drop1() below) rather than stating the significance of the two differences (+0 vs +2 and +0 vs +4) separately; report the overall $$p$$-value, then state the actual group values (e.g. 30.7, 16.6, 24.6, from emmeans::emmeans(m, ~Temp_offset, type="response"))
• try to avoid statements like "there is no difference between the base level and the +4 deg treatment" - I know you meant "no significant difference", but it's a slippery slope to fooling yourself that the two outcomes are actually identical, rather than just that you can't see a clear difference between them ...

### details

Data:

Trial2 <- data.frame(Brood_Number=c(
32, 47, 19, 23, 16, 45, 33, ## +0 deg
21, 12, 3,  22, 20, 6,  29, 20, ## +2 deg
16, 31, 14, 24, 30, 26, 40, 16), ## +4 deg
Temp_offset=rep(c("+0 deg","+2 deg", "+4 deg"), c(7,8,8)))
## generate rank information
Trial2$$Brood_rank <- rank(Trial2$$Brood_Number)


Compare raw values to ranked values. The Kruskal-Wallis is similar to an ANOVA on the ranks (in this case anova(lm(Brood_rank ~Temp_offset, data=Trial2)) gives $$p=0.070$$, pretty close to the K-W value.)

library(ggplot2); theme_set(theme_bw())
Trial2B <- tidyr::gather(Trial2,var_type,value,-Temp_offset)
ggplot(Trial2B, aes(Temp_offset, value))+geom_point()+
facet_wrap(~var_type,scale="free")


You can see that taking ranks pulls in the extreme values slightly, bringing the groups closer together and decreasing the strength of evidence for among-group differences.

Fit the model:

m <- MASS::glm.nb(Brood_Number~Temp_offset, data=Trial2)


Run tests:

drop1(m, test="Chisq")  ## p=0.021
kruskal.test(Brood_Number~Temp_offset, data=Trial2) ## p=0.077


Diagnostics:

(1) base R

png("nbplot1.png")
op <- par(mfrow=c(2,2), mar=c(2,3,2,1),las=1, mgp=c(1,0.5,0))
plot(m)
par(op)
dev.off()


(i) The residuals vs fitted plot looks OK, but it more or less has to (since there is no possibility for nonlinearity in a one-way categorical design); (ii) the Q-Q plot looks good, so little evidence for distributional differences; (iii) scale-location ditto - little evidence for heteroscedasticity; (iv) nothing that looks like an outlier or influential point.

An alternative set of diagnostics (using simulation-based residuals and among other things performing a significance test on the deviation of the distribution)

(2) DHARMa package

png("nbplot2.png")
DHARMa::simulateResiduals(m,plot=TRUE)
dev.off()