Extract the deviance from a list of GLM models using ldply()

Question

I want to extract the deviance from a list of GLM models using ldply()

Example data (from R base installation):

    library(reshape2)
    library(plyr)
    mtcars.1 <- mtcars[, c("am",  "qsec" , "drat")  ]
    mtcars.m <- melt(mtcars.1, id= 1 ) 

    glm.cars <- dlply( mtcars.m ,  .(variable) ,  
    glm,  formula=  am ~ value , family=binomial )

Got this far:

    ldply(  glm.cars  ,  summarise ,   "Null Deviance" = null.deviance , 
        "Residual Deviance" = deviance , "Deviance"= "??"    )

Which gives this:

      variable  Null Deviance     Residual Deviance    Deviance
1     qsec      43.22973          41.46512             ??
2     drat      43.22973          21.65003             ??

Deviance is missing! How do I extract it?

So how do I extract the deviance in the example above?

Off course i could just do null.deviance + deviance , but I just dont feel like doing it that way. I guess the reason I that I want to get better acquainted with the G statistic. I feel I go through the steops of extracting, subtracting and doing the chisqr, I will learn it better.

PS I was very confused to find that glm.model$devianc

Where do you see `Deviance` in the original data (that is, in `glm.cars`)? — David Robinson, Jan 25 '13 at 15:43
And what do you specifically mean by deviance (what is the difference between what you want and residual deviance?) — David Robinson, Jan 25 '13 at 15:47

IRTFM · Accepted Answer · 2013-01-25T16:24:49.213

As you say, you are confused. For each model you have two deviances. It's the difference in those two deviances (... not their sum) that is the interesting statistical measure. (I'm guessing you were making an analogy to the additive nature of residual sums of squares and model sums of squares, but if so, then you followed the wrong rabbit down the linguistic analogy hole.) You need to compare the difference to a 95% chi-square value with the same degrees of freedom as the difference in degrees of freedom between the null model and the "residual model". If you do str(.) on the models you can scroll down the list output to find that ther names are:

 str(glm(am~qsec, mtcars, family=binomial)  )
 .....
 $ deviance         : num 41.5
 $ aic              : num 45.5
 $ null.deviance    : num 43.2
 .....
 $ df.residual      : int 30
 $ df.null          : int 31
 .....

So your dlply code needs to extract those and and then you calculate null.deviance-deviance and df.null -df.residual and perhaps display qchisq(0.95, df.null-df.residual). If you want to see how it got packaged by the R-Core then look at:

 anova( glm(am~qsec, mtcars, family=binomial)  )

Ok, I think i got it now. The glm.cars$null.deviance is the -2*log(likelihood) of the model with the intercept only, and the glm.cars$deviance is the -2*log(likelihood) of the current model, in this case the model with the intercept plus a predictor? And i do null.deviance - deviance to get the G-statistic. Last question please: How do I get the p-value from the G-statistic pchisq(G , df.null - df ) wont do it :( — Rene Bern, Jan 28 '13 at 18:42
Exactly. P-value would just be `with(model_object, 1-pchisq(null.deviance-deviance , df.null -df.residual))` — IRTFM, Jan 28 '13 at 19:57

Extract the deviance from a list of GLM models using ldply()

1 Answers1