From groemping at bht-berlin.de  Mon Apr 26 20:48:15 2010
From: groemping at bht-berlin.de (=?ISO-8859-15?Q?Ulrike_Gr=F6mping?=)
Date: Mon, 26 Apr 2010 20:48:15 +0200
Subject: [Rd] Bug in calculation of overdispersion for quasibinomial grouped
 data ?
Message-ID: <4BD5DFEF.2030505@bht-berlin.de>

Dear list,

in preparing a lecture, I have created an example for the same data in 
grouped and ungrouped form (the well-known Titanic data). The code 
included below shows that the overdispersion estimates are strongly 
different for the two approaches.
If the overdispersion parameter FI is seen as a multiplyer for the 
variance, the estimate should the same from grouped and ungrouped data, 
as far as I can see. With the typical pearson residual approach to 
estimating FI and a mean response modeled with number of trials given in 
weights, this is not the case. The function overdisp.weighted below 
shows how estimation of overdispersion would have to be modified 
(intended for demonstrating the formula, not the programming) in a 
quasibinomial model with grouped data in order to yield the same 
estimate (up to a slight difference, presumably for numerical reasons) 
as for the ungrouped model.

Best regards,
Ulrike

## the standard way of estimating FI
overdisp <- function(fit) sum(residuals(fit, 
type="pearson")^2)/fit$df.residual
## the proposed way of estimating FI for grouped data
overdisp.weighted <- function(fit){
       yi.mean <- fit$model[,1]
       fiti <- fitted(fit)
       ni <- weights(fit)
       dfr <- sum(ni)-(fit$df.null-fit$df.residual+1)
       (sum(residuals(fit, type="pearson")^2) + 
sum(ni*yi.mean*(1-yi.mean)/(fiti*(1-fiti))))/dfr
}

### grouped data
require(alr3)
titanic

### Individual data
require(epitools)
titanic.lang <- expand.table(Titanic)
head(titanic.lang)

### quasi-binomial logistic model for grouped data
gruppmod <- glm(Surv/N~Class+Age+Sex, family=quasibinomial,
    data=titanic, weights=N)
summary(gruppmod)
overdisp(gruppmod)
overdisp.weighted(gruppmod)

### quasi-binomial logistic model for individual data
indivmod <- glm(Survived~Class+Age+Sex, family=quasibinomial,
    data=titanic.lang)
summary(indivmod)
overdisp(indivmod)

-- 
***********************************************
* Ulrike Groemping                            *
* BHT Berlin - University of Applied Sciences *
***********************************************
* +49 (30) 39404863 (Home Office)             *
* +49 (30) 4504 5127 (BHT)                    *
***********************************************
* http://prof.tfh-berlin.de/groemping         *
* groemping at bht-berlin.de                     *



[Rd] Bug in calculation of overdispersion for quasibinomial grouped data ?
=?ISO-8859-15?Q?Ulrike_Gr=F6mping?=