1
$\begingroup$

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

$\endgroup$
  • $\begingroup$ No mystery here: Assuming a negative binomial model provides information. The KW test does not use that information. $\endgroup$ – BruceET Jul 8 at 4:41
  • $\begingroup$ @BruceETCould you please clarify your comment? Do you mean that a NB uses additional info that the KW does not? $\endgroup$ – JeanDrayton Jul 8 at 5:11
  • $\begingroup$ any chance you can post your original data? It would make for a more complete/satisfying answer ... $\endgroup$ – Ben Bolker Jul 8 at 6:02
  • $\begingroup$ possible duplicate of stats.stackexchange.com/questions/264646/… ? $\endgroup$ – Ben Bolker Jul 8 at 6:02
  • 1
    $\begingroup$ @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 :) $\endgroup$ – JeanDrayton Jul 8 at 6:08
2
$\begingroup$

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")

enter image description here

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()

enter image description here

(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()

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.