zd<-function(X2,ctable,observed,expected) {
 #       X2 = value of Pearson's X2
 #   ctable = cross-classification table (as constructed by "crosstabs")
 # observed = observed counts
 # expected = expected counts (under independence)
 #
   NR <- nrow(ctable)
   NC <- ncol(ctable)
   row.sums <- apply(ctable,1,sum)
   col.sums <- apply(ctable,2,sum)
   total.sum<- sum(row.sums)

   Dstat <- X2 - sum(observed/expected)

   muD <- (total.sum*(NR-1)*(NC-1)/(total.sum-1)) - NR*NC

   part1 <-  2*total.sum/(total.sum - 3)
   part2 <- (NR-1)*(total.sum-NR)/(total.sum-1) -
            (total.sum*sum(1/row.sums) - NR^2)/(total.sum-2)
   part3 <- (NC-1)*(total.sum-NC)/(total.sum-1) -
            (total.sum*sum(1/col.sums) - NC^2)/(total.sum-2)
   part4 <- 4/(total.sum-2)*
            (total.sum*sum(1/row.sums) - NR^2)/(total.sum-2)*
            (total.sum*sum(1/col.sums) - NC^2)/(total.sum-2)

   sigma2D <- part1*part2*part3 + part4

   ZD <- (Dstat - muD)/sqrt(sigma2D)

  data.frame(D=Dstat,mu=muD,sigma=sqrt(sigma2D),Z=ZD,p=2*(1-pnorm(ZD)))
 }