6
$\begingroup$

I am trying to perform some logistic regressions (and I am a neophyte user of R). Initially I used "glm" to compute coefficients, AIC and p-values; this worked great until I ran across a data set suffering from complete separation. In [1], Gelman et alia suggest using an (informative) prior to address this problem; the corresponding algorithm is implemented in R as "bayesglm" (in the ARM package).

Here is my problem. Previously, with "glm", I would compute p-values as follows:

mylogit <- bayesglm(a ~ b+c+d+e+f+g+h, data = mydata, family="binomial")
with(mylogit, pchisq(null.deviance - deviance, df.null - df.residual, lower.tail = FALSE))

There are 53-48=5 degrees of freedom:

Null deviance: 71.188  on 53  degrees of freedom
Residual deviance: 37.862  on 48  degrees of freedom

However, if I use "bayesglm" instead of "glm", the resulting degrees of freedom are a bit surprising to me:

Null deviance: 22.279  on 53  degrees of freedom
Residual deviance: 39.030  on 54  degrees of freedom

If I plug in the preceding formula for a p-value, I have -1 degrees of freedom! Can someone help me get a more sensible answer (or help me interpret this)?

By the way, the documentation on the "bayesglm" command includes the following ominous comment:

We include all the glm() arguments but we haven’t tested that all the options (e.g., offests, contrasts, deviance for the null model) all work.

[1] Gelman, Andrew, et al. "A weakly informative default prior distribution for logistic and other regression models." The Annals of Applied Statistics (2008): 1360-1383.

$\endgroup$
6
$\begingroup$

I'm not sure how you got 5 degrees of freedom with seven independent variables for your glm-based model, but I'll assume that's just a typo somewhere or that I'm missing something minor.

Anyway, counting degrees of freedom with models that are constrained by a prior can be tricky, and there isn't necessarily a "correct" way to do it in many cases. Perhaps the authors of arm used -1 degrees of freedom as a way to keep people from blindly misinterpreting the results.

Although we can't easily calculate the number of degrees of freedom for most regularized models, we can at least put an upper bound on it: the number of degrees of freedom must be less than or equal to the degrees of freedom for the corresponding un-regularized model.

So (assuming the 5 degrees of freedom you reported above is correct), you can plug in 5 and be confident that the true P-value will be no larger than what your Chi-square test predicts. Thus, if it's significant with 5 degrees of freedom, the true value will also be significant.

If you want something more exact, you might want to look into using the lasso or ridge regression for regularization instead: statisticians have invested a lot of effort into counting degrees of freedom for these models, and have even developed some significance tests for them. Andrew Gelman talks about one recent advance on his blog here.

Edited to add: If you do stick with bayesglm but don't trust the null deviance estimates, you can find it yourself by running a model with no predictors except the intercept. The formula syntax for this would be a ~ 1.

$\endgroup$
  • 1
    $\begingroup$ Thanks for the advice, that's very helpful! And I appreciate your suggestion for double-checking the null deviance estimates. As for the degrees of freedom, that is not a typo-- although there are 7 (probabilistically) independent variables, the rank of the matrix of observations is only 5 (i.e. the observations are linearly dependent). Thanks! $\endgroup$ – Bill Bradley Oct 16 '13 at 14:29

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.