# # seizure3.q # # ---------------------------------------------------------------------- # # PURPOSE: example analysis for B571 # # AUTHOR: P. Heagerty # # DATE: 22 Jan 2003 # # ---------------------------------------------------------------------- # seizure <- read.table("seizure.dat") seizure[1:10,] names(seizure) # options( contrasts="contr.treatment" ) # seizure$logA <- log( seizure$age ) seizure$logB <- log( seizure$base + 0.5 ) model.matrix(~ age + base + trt, data=seizure)[1:10,] # ##### ##### fit a poisson regression ##### # glmfit <- glm( y4 ~ logA + logB + trt, family = "poisson", data = seizure) summary( glmfit, cor=F) # ##### ##### what about scale parameter? ##### # qfit <- glm(y4 ~ logA + logB + trt, family = quasi(link=log,variance=mu), data = seizure ) summary( qfit, cor=F) # ##### ##### negative binomial ML ##### # source("negbin.q") # nbfit <- glm.nb( y4 ~ logA + logB + trt, data = seizure ) summary( nbfit, cor=F) # mmm <- fitted( nbfit ) ddd <- rep( nbfit$theta , length(mmm) ) y <- seizure$y4 logL <- sum( lgamma( ddd + y ) - lgamma( ddd ) - lgamma( y+1 ) + ddd * log( ddd/(ddd+mmm) ) + y * log( mmm/(ddd+mmm) ) ) logLb <- sum( lgamma( ddd + y ) - lgamma( ddd ) + ddd * log( ddd/(ddd+mmm) ) + y * log( mmm/(ddd+mmm) ) ) # print( nbfit$twologlik ) print( 2*logL ) print( 2*logLb ) # # end-of-file