Ulrike Grömping
2010-Apr-26 18:48 UTC
[Rd] Bug in calculation of overdispersion for quasibinomial grouped data ?
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 *