# # weil_QL.q # # -------------------------------------------------- # # PURPOSE: Use quasilikelihood for Weil data # # AUTHOR: P. Heagerty # # DATE: 00/01/26 # # -------------------------------------------------- # data <- matrix(scan("weil.data"), ncol = 3, byrow = T) weil <- list( y = data[ , 1], n = data[ , 2], tx = data[ , 3], x = cbind(1, data[ , 3]), z = matrix(1, nrow(data), 1), id = c(1:nrow(data)) ) # ##### ##### scale variance model ##### # glmfit <- glm( cbind(y,n-y) ~ tx, family=binomial, data = weil ) summary( glmfit, cor=F ) # qlfit1 <- glm( cbind(y,n-y) ~ tx, family=quasi(link="logit", variance="mu(1-mu)"), data = weil ) summary( qlfit1, cor=F ) # qlfit2 <- glm( cbind(y,n-y) ~ 1, family=quasi(link=logit, variance="mu(1-mu)"), data = weil ) summary( qlfit2 ) # anova( qlfit1, qlfit2, test="F" ) # ##### ##### beta-binomial variance model ##### # ##### Iteration required! # fit <- glmfit resids <- residuals( glmfit, type="pearson" ) alpha <- mean( (resids^2 - 1)/(weil$n-1) ) cat( paste("alpha =", alpha,"\n") ) # tolerance <- 1e-6 converged <- F maxiter <- 10 # count <- 0 while( !converged && (count