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 *
* groemp...@bht-berlin.de *
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
