having constructed a glm we can pass them through an anova as such:

    anova(mnegbin1, mnegbin2, test = "Chisq")
    Likelihood ratio tests of Negative Binomial Models
    
    Response: count
                            Model    theta Resid. df    2 x log-lik.   Test    df      LR stat. Pr(Chi)
    1 origin + substrate + sample 18.46502      1232       -2459.886                                   
    2          substrate + sample 18.46502      1232       -2459.886 1 vs 2     0 -1.500666e-11       1

This gives us a comparison of the log likelihoods between the models, but can something similar be done for between residual deviances?:

    summary(mnegbin1)
        Null deviance: 5686.16  on 1304  degrees of freedom
        Residual deviance:  890.24  on 1232  degrees of freedom

    summary(mnegbin2)
        Null deviance: 5686.16  on 1304  degrees of freedom
        Residual deviance:  890.24  on 1232  degrees of freedom

In this instance they are of course identical but does that actually tell us anything?