7

So, I'm trying to compare two models, fit1 and fit2.

Initially, I was just doing anova(fit1,fit2), and this yielded output that I understood (including a p-value).

However, when I switched my models from lm()-based models to glm()-based models, anova(fit1,fit2) now yielded Residual Degrees of Freedom, Residuals Deviances, and Df Deviances, which I am having trouble interpreting (resources explaining these metrics seem scarce). I was hoping to extract a p-value for the comparison between the two models, but for some reason anova(fit1,fit2, test='Chisq') isn't working. Any suggestions?

I realize that, depending on the link function in my glms, Chi-squared may not be the most appropriate test, but I have used 'F' in appropriate contexts as well with similar disappointment.

Is this problem familiar to anybody else? Suggestions? Many thanks!

Example:

make_and_compare_models <- function(fitness_trait_name, data_frame_name, vector_for_multiple_regression, predictor_for_single_regression, fam){
        fit1<-glm(formula=as.formula(paste(fitness_trait_name,"~", paste(vector_for_multiple_regression, sep="+"))), family=fam, data=data_frame_name)
        print ("summary fit 1")
        print(summary(fit1))
        fit2<- glm(data=data_frame_name, formula=as.formula(paste(fitness_trait_name,"~",predictor_for_single_regression)), family=fam)

        print("summary fit 2")
        print(summary(fit2))
        print("model comparison stats:")
        mod_test<-anova(fit2,fit1)

        ##suggestion #1
        print(anova(fit2,fit1, test="Chisq"))

        #suggestion #2
        print ("significance:")
    print (1-pchisq( abs(mod_test$Deviance[2]),df=abs(mod_test$Df[2])))

        }


data<-structure(list(ID = c(1L, 2L, 4L, 7L, 9L, 10L, 12L, 13L, 14L, 
15L, 16L, 17L, 18L, 20L, 21L, 22L, 23L, 24L, 25L, 27L, 28L, 29L, 
31L, 34L, 37L, 38L, 39L, 40L, 41L, 43L, 44L, 45L, 46L, 47L, 48L, 
49L, 52L, 55L, 56L, 59L, 60L, 61L, 62L, 63L, 65L, 66L, 67L, 68L, 
69L, 71L), QnWeight_initial = c(158L, 165L, 137L, 150L, 153L, 
137L, 158L, 163L, 159L, 151L, 145L, 144L, 157L, 144L, 133L, 148L, 
151L, 151L, 147L, 158L, 178L, 164L, 134L, 151L, 148L, 142L, 127L, 
179L, 162L, 150L, 151L, 153L, 163L, 155L, 163L, 170L, 149L, 165L, 
128L, 134L, 145L, 147L, 148L, 160L, 131L, 155L, 169L, 143L, 123L, 
151L), Survived_eclosion = c(0L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 
1L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), Days_wrkr_eclosion_minus20 = c(NA, 
1L, NA, 3L, 0L, 2L, 0L, 1L, 0L, 0L, 0L, 1L, NA, 0L, 7L, 1L, 0L, 
1L, 0L, 1L, 2L, 2L, NA, 2L, 3L, 2L, 2L, NA, 0L, 1L, NA, NA, 0L, 
0L, 0L, 0L, 3L, 3L, 3L, 1L, 0L, 2L, NA, 1L, 0L, 1L, 1L, 3L, 1L, 
2L), MLH = c(0.5, 0.666666667, 0.555555556, 0.25, 1, 0.5, 0.333333333, 
0.7, 0.5, 0.7, 0.5, 0.666666667, 0.375, 0.4, 0.5, 0.333333333, 
0.4, 0.375, 0.3, 0.5, 0.3, 0.2, 0.4, 0.875, 0.6, 0.4, 0.222222222, 
0.222222222, 0.6, 0.6, 0.3, 0.4, 0.714285714, 0.4, 0.3, 0.6, 
0.4, 0.7, 0.625, 0.555555556, 0.25, 0.5, 0.5, 0.6, 0.25, 0.428571429, 
0.3, 0.25, 0.375, 0.555555556), Acon5 = c(0.35387674, 0.35387674, 
0.35387674, 0.35387674, 0.35387674, 0.35387674, 0.35387674, 0, 
0, 1, 0, 1, 0.35387674, 0, 0, 0.35387674, 1, 1, 0, 0, 0, 1, 0, 
0.35387674, 0, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 
0, 0, 1, 0, 0, 0, 1, 0, 0.35387674), Baez = c(1, 1, 1, 0.467836257, 
1, 1, 0, 0, 1, 1, 0, 0.467836257, 1, 0, 0, 0, 0, 1, 0, 0, 0, 
0, 0, 0.467836257, 1, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 1, 
1, 0, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1), C294 = c(0, 1, 0, 0, 1, 
0.582542694, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 
0, 1, 1, 0, 0, 0.582542694, 1, 1, 0, 0, 1, 0, 1, 1, 0, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1), C316 = c(1, 1, 0, 0, 0.519685039, 
0.519685039, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0.519685039, 0, 
1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0.519685039, 1, 0, 1, 
1, 0, 0.519685039, 1, 0.519685039, 1, 1, 1, 0.519685039, 0.519685039, 
0, 0.519685039, 0.519685039, 0), i_120_PigTail = c(1, 1, 0, 1, 
0.631236443, 0.631236443, 1, 1, 1, 1, 1, 0, 0.631236443, 1, 1, 
1, 0, 0.631236443, 1, 1, 1, 0, 0, 1, 1, 1, 0.631236443, 0, 1, 
1, 0, 1, 0.631236443, 1, 0, 1, 0, 0, 1, 0.631236443, 0.631236443, 
0, 1, 0, 0.631236443, 0.631236443, 1, 0.631236443, 0.631236443, 
1), i129 = c(0L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 0L, 1L, 0L, 0L, 
1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 
0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L), Jackstraw_PigTail = c(0L, 1L, 1L, 0L, 
1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 
1L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 
0L, 1L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), Neil_Young = c(0.529636711, 
0, 1, 0, 0.529636711, 0.529636711, 1, 1, 0, 1, 1, 1, 0, 0, 1, 
1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 
1, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1), Ramble = c(0, 0, 0, 
0, 0.215163934, 0.215163934, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 
0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0.215163934, 0, 
0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0.215163934, 0, 0, 0, 0), Sol_18 = c(1, 
0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 
0, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0.404669261, 
1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1)), .Names = c("ID", "QnWeight_initial", 
"Survived_eclosion", "Days_wrkr_eclosion_minus20", "MLH", "Acon5", 
"Baez", "C294", "C316", "i_120_PigTail", "i129", "Jackstraw_PigTail", 
"Neil_Young", "Ramble", "Sol_18"), class = "data.frame", row.names = c(NA, 
-50L))

make_and_compare_models("QnWeight_initial", data, c("Acon5","Baez","C294","C316","i_120_PigTail","i129","Jackstraw_PigTail","Neil_Young","Ramble","Sol_18"), "MLH", "gaussian")
Atticus29
  • 4,190
  • 18
  • 47
  • 84
  • 3
    `anova(fit1,fit2,test="Chisq")` **should** work, unless the nested models happen to have identical fits. Can you provide more detail? – Ben Bolker Dec 20 '12 at 14:04
  • PS it's not the link function but the family that determines whether you should use Chi-square or F (specifically, whether the scale parameter is fixed [Poisson, binomial] or estimated [Gaussian, Gamma, quasi-likelihood fits] – Ben Bolker Dec 20 '12 at 14:06
  • @BenBolker thanks for the clarification. Just to be sure, it Chi-square for fixed scale parameters and F for estimated? Also, the output from anova(fit1,fit2, test="Chisq") yields a Pr( – Atticus29 Dec 20 '12 at 22:21
  • Also, is there a test="F" analogue? I can't find anything about test as a parameter for anova() in the manual... – Atticus29 Dec 20 '12 at 22:22
  • (1) Yes, chi-square for fixed and F for estimated. (2) Can you provide a reproducible example of the pattern you're describing? It would be **really** helpful One confusing "feature" is that if the deviance difference is less than or approximately zero (i.e. the models are essentially identical), then R doesn't print the p-value, and what you're really seeing is the deviance difference (e.g. see @DWin's example). (3) For "F", "Chisq", etc. see `?anova.glm`, which in turn references `?stat.anova`. – Ben Bolker Dec 21 '12 at 03:05
  • Ok @BenBolker, I have added an example. I hope that this provides you with the information you need... – Atticus29 Dec 21 '12 at 17:53
  • 1
    Your example shows that you are comparing *non-nested* models: the df difference (shown in the `Df` column) is zero! All of the `anova()` framework (as discussed in the answers below) is framed around *nested* models. If you want to compare goodness-of-fit of *non*-nested models, you can use AIC (with caution) or the Vuong test ... – Ben Bolker Dec 21 '12 at 17:55
  • I will post more of an answer later, but I also question your *statistical* approach -- this might be evolving into a http://stats.stackexchange.com question, but question are you trying to answer?? – Ben Bolker Dec 21 '12 at 18:11
  • I am taking the advice of statistical geneticist, who wrote about this method in his Evolution paper in 1997. "The appropriate procedure is to test whether a multiple regression incorporating specific effects for each locus explains more variance than a simple regression on MLH". I don't think that the statistical approach is problematic... – Atticus29 Dec 21 '12 at 18:18
  • Citation is David, P. 1997. Modeling the genetic basis of heterosis: tests of alternative hypotheses. Evolution 51:1049–1057. – Atticus29 Dec 21 '12 at 18:18
  • Per a personal discussion with David a few months ago, "You don’t need to do that because the models are actually nested even if it’s not apparent. MLH is the sum of all single-locus heterozygosities so the regression reads fitness = alpha MLH + mu = alpha (H1+H2+H3…)+ mu the multiple regression with all single locus heterozygosities reads fitness = alpha1 H1 + alpha2 H2 + alpha3 H3 …+mu therefore the MLH regression is simply a special case of the multiple regression in which we impose alpha1=alpha2=alpha3=…; so the models ARE already nested." – Atticus29 Dec 21 '12 at 18:23
  • let us [continue this discussion in chat](http://chat.stackoverflow.com/rooms/21535/discussion-between-ben-bolker-and-atticus29) – Ben Bolker Dec 21 '12 at 18:28

2 Answers2

8

The difference in deviance between a "larger" or more complex model and a nested or "reduced" model is distributed (asymptotically) as a chi-squared variate with the difference in degrees of freedom of the two models. So you would extract the deviance estimate and the difference in degrees of freedom and compare that to pchisq( deviance, diff(df) ). The "p-value" is just 1 minus that value.

> 1-pchisq(3.84,1)
[1] 0.05004352

If you run the first example in the glm help page and then add a reduced model without the "treatment" variable, you get:

glm.D93.o <- glm(counts ~ outcome, family=poisson())
 anova.res <-anova(glm.D93, glm.D93.o)
 anova.res
#------------
Analysis of Deviance Table

Model 1: counts ~ outcome + treatment
Model 2: counts ~ outcome
  Resid. Df Resid. Dev Df    Deviance
1         4     5.1291               
2         6     5.1291 -2 -2.6645e-15
#---------------
 str(anova.res)
Classes ‘anova’ and 'data.frame':   2 obs. of  4 variables:
 $ Resid. Df : num  4 6
 $ Resid. Dev: num  5.13 5.13
 $ Df        : num  NA -2
 $ Deviance  : num  NA -2.66e-15
 - attr(*, "heading")= chr  "Analysis of Deviance Table\n" "Model 1: counts ~ outcome + treatment\nModel 2: counts ~ outcome"

So after looking at how things were stored in the object itself, this give the p-value for "outcome":

 1-pchisq( abs(anova.res$Deviance[2]), abs(anova.res$Df[2]))
[1] 1

And this would be the corresponding procedure on the treatment+outcome model versus the treatment-only model:

> glm.D93.t <- glm(counts ~ treatment, family=poisson())
> anova.res2 <-anova(glm.D93, glm.D93.t)
> 1-pchisq( abs(anova.res2$Deviance[2]), abs(anova.res2$Df[2]))
[1] 0.06547071
IRTFM
  • 258,963
  • 21
  • 364
  • 487
  • Thanks, DWin! That answers my question! – Atticus29 Nov 05 '12 at 18:49
  • the 1-pchisq() can't be right. I've run simulations with completely scrambled data (i.e., there should be no significant difference between the two models, because neither model successfully predicts the response), and the reported p-value is consistently "0". Are you sure it's not just pchisq() in this case? – Atticus29 Dec 19 '12 at 20:47
  • I am quite sure that `1-pchisq(3.84,1)` returns 0.05. You need to make sure you are putting the absolute value of the correct deviance difference in the first argument and the correct degrees of freedom in the second. The order of the model arguments will flip the sign of the anova `$Deviance` results but the `abs()` should take care of that. – IRTFM Dec 19 '12 at 21:23
  • Point taken. Absolute values are there. Hmm... ok I just specifically designated the second argument as "df=model$DF[2]", and that cleared it up. Interesting... – Atticus29 Dec 19 '12 at 22:03
  • This is a somewhat pathological example! For reasons that I don't yet understand, the `treatment` variable is redundant (has exactly zero predictive power), so R doesn't print the chisq p-value even when it is requested. `glm.D93.i <- glm(counts~1,family=poisson); anova(glm.D93.i,glm.D93.o,test="Chisq")` is somewhat easier to understand. – Ben Bolker Dec 20 '12 at 14:03
1

If your 2 models are nested, then you can use the change in deviance of the 2 models to see if the model containing extra parameters yields an improved fit. If model 1 contains k parameters and model 2 contains those same k parameters plus an additional m parameters, then the change in deviance follows an (approximately) chi-square distribution with m degrees of freedom. You can use this test statistic to see if model 2 is an improvement on model 1.

If you are new to this area, I would strongly recommend reading an introductory text on GLMs

mathematician1975
  • 21,161
  • 6
  • 59
  • 101
  • that's perfect, except I'm not quite sure how to actually implement that. I.e., do you happen to know the R syntax for it? – Atticus29 Nov 05 '12 at 18:23
  • Unfortunately it is years since I used R. As far as I recall, the glm.summary output used to provide everything that was needed for this calculation. Hopefully you will get an `R` specific answer rather than just theoretical. – mathematician1975 Nov 05 '12 at 18:27