.ND
.pl 100i
.po 3
.na
.SH


.SH
Usage:
.LP
.nf
.in +0.5i

data(Growth)
data(budworm)

.in -0.5i
.fi

.SH
Details:
.IP "" 0.5i
These are examples of counting data analysis.
.in -0.5i

.SH
References:
.IP "" 0.5i
Venables and Ripley (1999) Modern Applied Statistics
with S-PLUS, 3rd ed., ch. 7.
.in -0.5i

.SH
Examples:
.LP
.nf
.in +0.5i

# attach library; get data
#library( pda )
data( Growth )

Growth <- Growth[,c("trt","advplt")]
# drop the added BHTA control ( trt = 20 )
Growth <- Growth[ Growth$trt != 20, ]
 
Growth$BA <- c(rep(0,4),rep(0.44,4),rep(4.4,3),rep(NA,9),0)[1+Growth$trt]
Growth$TDZ <- c(rep(c(0,0.2,2,20),3)[1:11],rep(NA,9),0)[1+Growth$trt]
Growth$code <- c(0:9,"a",rep(NA,9),"C")[1+Growth$trt]

Growth$trt <- factor(Growth$trt)
Growth$BA <- factor(Growth$BA)
Growth$TDZ <- factor(Growth$TDZ)
Growth <- Growth[!is.na(Growth$advplt),]

print( sample.size( Growth$BA, Growth$TDZ ))
hist( Growth$advplt, nclass = 25 )
print( table( Growth$advplt ))
Growth$nonzero <- as.numeric( Growth$advplt > 0 )

Growth.bin <- glm( nonzero ~ BA * TDZ, data = Growth,
   family = binomial )
print( anova( Growth.bin, test = "Chisq" ))
Growth.bin <- glm( nonzero ~ TDZ *  BA, data = Growth,
   family = binomial )
print( anova( Growth.bin, test = "Chisq" ))
print( drop1( Growth.bin, formula( Growth.bin),
   test = "Chisq" ))

Growth.des <- data.frame(
   BA = factor( rep( levels( Growth$BA ), 4 )),
   TDZ = factor( rep( levels( Growth$TDZ ), rep(3,4) )))
Growth.pred <- 1 / ( 1 + exp( - predict(
   Growth.bin, Growth.des )))
plot( c(1,21),range( Growth.pred ), type = "n", log = "x",
   xlab = "TDZ", ylab = "fraction nonzero" )
for( i in levels( Growth.des$BA ))
{
   Growth.BA <- i == Growth.des$BA
   text( 1 + unfactor( Growth.des$TDZ[Growth.BA] ),
      Growth.pred[Growth.BA], i )
}
rm( Growth.BA )

# library( vr )
data( budworm )
budworm$dead <- budworm$dead / 20
budworm$total <- rep( 20, nrow( budworm ))
budworm.fit <- glm( dead ~ sex * ldose, data = budworm,
   weight = total, family = binomial )
summary( budworm.fit )

# Plot of Logit by Sex

plot( c(1,32), c(0,1), type = "n", xlab = "dose",
   ylab = "prob", log = "x" )
text( 2^budworm$ldose, budworm$dead, as.character( budworm$sex ))
budworm.ld <- seq( 0, 5, 0.1 )
for( i in levels( budworm$sex ))
   lines( 2^budworm.ld, predict( budworm.fit,
      data.frame( ldose = budworm.ld,
         sex = factor( rep (i, length( budworm.ld )),
         levels = levels( budworm$sex ))),
      type = "response" ))
rm( budworm.ld )

# Parameters for Lines

print( summary( glm( dead ~ sex + ldose - 1, data = budworm,
   weight = total, family = binomial ))$coefficients )

# Four-way Contingency Table

brandx <- cbind( expand.grid( Brand = c("X","M"),
   Temp = c("Low","High"), M.user = c("N","Y"),
   Soft = c("Hard","Medium","Soft")),
   Fr = c(68,42,42,30,37,52,24,43,
        66,50,33,23,47,55,23,47,
        63,53,29,27,57,49,19,29))
brandx$Soft <- ordered( brandx$Soft,
 levels = c("Soft","Medium","Hard"))
contrasts( brandx$Brand ) <- contr.poly
contrasts( brandx$Temp ) <- contr.poly
contrasts( brandx$M.user ) <- contr.poly
contrasts( brandx$Soft ) <- contr.poly

# fit poisson model

brandx.fit <- glm( Fr ~ M.user*Temp*Soft+Brand,
   family = poisson, data = brandx )
print( anova( brandx.fit, test = "Chisq" ))
print( drop1( brandx.fit, formula( brandx.fit), test = "Chisq" ))

brandx.step <- step( brandx.fit,
   list( lower = formula( brandx.fit ), upper = ~.^3 ),
   scale = 1, trace = F )
print( brandx.step$anova )
print( anova( brandx.step, test = "Chisq" ))
brandx.mod <- glm( terms( Fr ~ M.user*Temp*Soft +
   Brand*M.user*Temp, keep.order = T ),
   family = poisson, data = brandx )
summary(brandx.mod, correlation = F, test = "Chisq" )

.in -0.5i
.fi