R help - May 2008 - Estimating QAIC using glm with the quasibinomial family

Hello R-list.  I am a "long time listener - first time caller" who has
been using R in research and graduate teaching for over 5 years.  I  
hope that my question is simple but not too foolish.  I've looked  
through the FAQ and searched the R site mail list with some close hits  
but no direct answers, so...
I would like to estimate QAIC (and QAICc) for a glm fit using the  
quasibinomial family.  I found a general reference suggesting a simple  
solution:

"we calculated QAICc adjusting for overdispersion by dividing the  
residual deviance (i.e. -2 loglikelihood) with the overdispersion  
parameter calculated from the most complex model as the sum of squares  
Pearson residuals divided by the number of degrees of freedom (Burnham  
& Anderson, 2002). "
- Mystrud et al. 2007. Animal Conservation. 10:77-87.


My question is: Will this calculation be valid with the residual  
deviance returned by the glm() function using the quasibinomial family  
as reported in R?


I thought I should ask to be certain that there is no dispersion  
correction applied to the reported deviance, as encouraged by Burnham  
and Anderson, 2nd ed., 2002 on p.69:

"When data are overdispersed and c > 1, the proper likelihood is
log(L)/c".

Regards,  Darren Gillis

Department of Biological Sciences
Faculty of Science
University of Manitoba
Winnipeg, MB
Canada, R3T 2N2

R-users
E-mail: r-help@r-project.org
>My question is: Will this calculation be valid with the residual deviance
>returned by the glm() function using the quasibinomial family as
>reported in R?
   Let me show you a simple example, assuming c=2.5:

function ()
{
   xx <- c(1,2,3,4,5,6,7,8,9,10)
   yy <- c(1,0,1,0,0,1,0,1,1,1)
   data1 <- data.frame(x=xx, y=yy)
   out1 <- glm(y~x, data=data1, family=binomial)
   print(out1)
   aic0 <- out1$aic
   print("aic0")
   print(aic0)

   dev1 <- out1$deviance
   aic1 <-  dev1+ 2*2
   print("aic1")
   print(aic1)

   c1 <- 2.5
   qaic1 <-  dev1/c1+ 2*2
   print("qaic1")
   print(qaic1)
}

The result is:

Call:  glm(formula = y ~ x, family = binomial, data = data1)

Coefficients:
(Intercept)            x
    -0.7300       0.2131

Degrees of Freedom: 9 Total (i.e. Null);  8 Residual
Null Deviance:      13.46
Residual Deviance: 12.63        AIC: 16.63
[1] "aic0"
[1] 16.63054
[1] "aic1"
[1] 16.63054
[1] "qaic1"
[1] 9.052216

   I hope that this R program will be of some help to you.

K. Takezawa

	[[alternative HTML version deleted]]

Thank you for the rapid response.  I have no trouble in following your
example, but I was specificly interested in calculating QAIC and QAICc from
a glm fitted with the "family=quasibinomial" option.  This provides an
estimate of dispersion, but does not contain an explicit value for AIC or
log-likelihood.  Do you have any insights into how to deal with the output
from this specific option?  Or are you suggesting not to use
family=quasibinomial when faced with potential overdispersion?  Regardless,
thank you again for your response.

Cheers, Darren


-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
On
Behalf Of Kunio takezawa
Sent: Tuesday, May 06, 2008 8:41 PM
To: r-help at r-project.org
Subject: Re: [R] Estimating QAIC using glm with the quasibinomial family

R-users
E-mail: r-help at r-project.org
>My question is: Will this calculation be valid with the residual 
>deviance returned by the glm() function using the quasibinomial family 
>as reported in R?
   Let me show you a simple example, assuming c=2.5:

function ()
{
   xx <- c(1,2,3,4,5,6,7,8,9,10)
   yy <- c(1,0,1,0,0,1,0,1,1,1)
   data1 <- data.frame(x=xx, y=yy)
   out1 <- glm(y~x, data=data1, family=binomial)
   print(out1)
   aic0 <- out1$aic
   print("aic0")
   print(aic0)

   dev1 <- out1$deviance
   aic1 <-  dev1+ 2*2
   print("aic1")
   print(aic1)

   c1 <- 2.5
   qaic1 <-  dev1/c1+ 2*2
   print("qaic1")
   print(qaic1)
}

The result is:

Call:  glm(formula = y ~ x, family = binomial, data = data1)

Coefficients:
(Intercept)            x
    -0.7300       0.2131

Degrees of Freedom: 9 Total (i.e. Null);  8 Residual
Null Deviance:      13.46
Residual Deviance: 12.63        AIC: 16.63
[1] "aic0"
[1] 16.63054
[1] "aic1"
[1] 16.63054
[1] "qaic1"
[1] 9.052216

   I hope that this R program will be of some help to you.

K. Takezawa

	[[alternative HTML version deleted]]

______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.