az.fit<-glm(count~sex+dose,family=poisson(link=log),
             data=azin.df)
 expected<-predict(az.fit,type="response")

 power.diverg<-function(observed,expected,NR=2,NC=2,lambda=1) {
 # Power-Divergence Statistic
 # NR (NC) = # of rows (columns) in the table
 # lambda = 1 (default to Pearson X2)
   c2.lambda<- 2/(lambda*(lambda+1)) *
         sum(observed*((observed/expected)^lambda - 1))
   p.value<- 1 - pchisq(c2.lambda,df=(NR-1)*(NC-1))
   data.frame(C2.lambda=c2.lambda,lambda=lambda,p.value=p.value)
 }

 # Pearson X2 ==> lambda = 1
 power.diverg(azin.df$count,expected,
              NR=length(unique(azin.df$sex)),
              NC=length(unique(azin.df$dose)))


    C2.lambda   lambda      p.value
 1  22.60906      1       1.231699e-05


 # Cressie-Read ==> lambda = 2/3
 power.diverg(azin.df$count,expected,lambda=2/3,
              NR=length(unique(azin.df$sex)),
              NC=length(unique(azin.df$dose)))


    C2.lambda    lambda     p.value
 1  22.38332   0.6666667  1.37887e-05


 # Power-divergence ==> lambda = 1/2
 power.diverg(azin.df$count,expected,lambda=1/2,
              NR=length(unique(azin.df$sex)),
              NC=length(unique(azin.df$dose)))


    C2.lambda    lambda     p.value
 1  22.28192      0.5     1.450583e-05