Simulation

### Functions
test_stat <- function(x){
  ts <- numeric()
  
  X <- matrix(x,r)
  rs <- rowSums(X)  # ni.
  cs <- colSums(X)  # n.j
  
  # E(n)  
  M <- outer(rs,cs)/sum(X)  
    
  # 1st test statistic (LRT)
  ts[1]<-1-pchisq(2*sum(X*log(X/M),na.rm=T),(r-1)*(c-1));

  # 2nd test statistic (Pearson's Chi-square)
  ts[2]<- 1-pchisq(sum((X-M)^2/M),(r-1)*(c-1));
  
  ts[3:12]<-exact(X,rs,cs);
  
  return(ts)
}
exact <- function(x,rs,cs){
  T <- r2dtable(iter,rs,cs)
  obs <- unlist(mmh(x))
  exp <- sapply(T,mmh)

  return(apply(exp>obs,1,mean))
}
mh<-function(x,i,j){
  N<-x[1,1]+x[1,j]+x[i,1]+x[i,j];
  return((x[i,j]*x[1,1]/N-x[i,1]*x[1,j]/N)/sqrt((x[1,1]+x[i,1])/N*(x[1,1]+x[1,j])/N*(x[i,j]+x[i,1])/(N-1)*(x[i,j]+x[1,j])))
}
mmh<-function(x){
  Z1<-mh(x,2,2);  Z2<-mh(x,2,3);  Z3<-mh(x,3,2);  Z4<-mh(x,3,3)
  Z_one<-max(Z1,Z2,Z3,Z4)
  Z_two<-max(abs(Z1),abs(Z2),abs(Z3),abs(Z4))
  return(list(Z_one,Z1,Z2,Z3,Z4,
              Z_two,abs(Z1),abs(Z2),abs(Z3),abs(Z4)))
}


### sample size and number of iteration
n <- 100; r<-3; c<-3; s <- 1000; iter <- 1000; 

### test statistics
num_ts<-12; names_ts <- c('QL', 'Qp',
                          'One-sided', 'MH22', 'MH23', 'MH32','MH33',
                          'Two-sided', 'MH22', 'MH23', 'MH32','MH33')

### parameters
alpha <- c(1,0.2,0.2); beta <- c(1,0.2,0.1); gamma <- matrix(c(1,1,1,1,3.,1.,1,1.,1.),r)

### s 2*2 random tables with n from multinomial distribution
tables <- rmultinom(s, n, outer(alpha,beta)*gamma)

p <- apply(tables,2,test_stat)
rownames(p)<-names_ts

### empirical tail probability
percent <- cbind(apply((p<=0.05),1,mean,na.rm=TRUE))
rownames(percent)<- names_ts

write.csv(rbind(alpha,beta,gamma),"E:/documents/maxor/out/sim.csv",append=TRUE)
write.csv(percent,"E:/documents/maxor/out/sim.csv",append=TRUE)