# Trend tests under extra-binomial variability

 Z.GEE.bin<-function(conc,yy,rr)   {
 # formula (6.12)
 # conc = vector of concentrations
 # yy = vector of observed events
 # rr = vector of replicate counts

          ngrp <- length(unique(conc))
      conc.lst <- split(conc,conc)
        yy.lst <- split(yy,conc)
        rr.lst <- split(rr,conc)

          xbar <- sum(rr*conc)/sum(rr)
          pbar <- sum(yy)/sum(rr)

    top <- sum( (conc-xbar)*yy )
    bot <- 0
    for (ii in 1:ngrp) {
       bot <- bot + sum( ((conc.lst[[ii]]-xbar)^2)*
              ((yy.lst[[ii]] - rr.lst[[ii]]*pbar)^2))
  }

    cat("xbar = ",signif(xbar),"\n")
    cat("pbar = ",signif(pbar),"\n")
    cat("Z.GEE= ",signif(top/sqrt(bot)),"\n")
    top/sqrt(bot)
 }

 # place the teratogenicity data into a data frame

 dead   <-c(dead.0,dead.1,dead.2,dead.4)
 implant<-c(implant.0,implant.1,implant.2,implant.4)
 conc   <-rep(c(0,.1,.2,.4),c(27,27,27,26))

 terat.df<-data.frame(conc=conc,dead=dead,implant=implant)

 Z.GEE<-Z.GEE.bin(terat.df$conc,terat.df$dead,terat.df$implant)

# output follows ...

 xbar =   0.172 
 pbar =   0.102 
 Z.GEE=   2.427

 1-pnorm(Z.GEE)
 [1] 0.007611753