n.in.riskset=c(218,290,265,182,183,170,213,239,115,219,132,219,164,199,187,154,84);
risksets = length(n.in.riskset);
prop.exposed=c(.33,.26,.25,.36,.32,.19,.29,.38,.25,.40,.18,.40,.31,.39,.35,.23,.23);

n.exposed     = round(n.in.riskset*prop.exposed); 
n.not.exposed = n.in.riskset - n.exposed ;

case.exposure    = c(1.26,0,.75,4.3,2.76,.94,0,0,0,.37,0,7.88,2.41,0,0,.39,0);
case.exposed     = (case.exposure > 0);
case.not.exposed = 1 - case.exposed;

controls.exposed     = n.exposed     - case.exposed; 
controls.not.exposed = n.not.exposed - (1-case.exposed);

mhNum=sum( case.exposed     * ctls.not.expposed / n.in.riskset ) ;
mhDen=sum( case.not.exposed * ctls.expposed    / n.in.riskset  ) ;
c("mhNum",round(mhNum,2),"mhDen",round(mhDen,2),"mhNum/mhDenmhDen", round(mhNum/mhDen,2))

require(survival)

datasets=201;

sizes = c(4,9,16,25,36); n.sizes=length(sizes);
Q1.beta.hat  = c();
median.beta.hat = c();
Q3.beta.hat  = c();
for (ctls.per.case in sizes  )
{

   beta.hat=c();

   for (dataset in 1:datasets)
   {
      e=c(); case=c(); which.set=c()
      for (set in 1:risksets) {
         which.set=c(which.set,set);
         case = c(case,1); 
         e = c(e,case.exposed[set]);
         which.set = c(which.set,rep(set,ctls.per.case));
         case = c(case,rep(0,ctls.per.case));
         e = c(e, sample( c( rep(1,controls.exposed[set]),
                             rep(0,controls.not.exposed[set]) ), ctls.per.case) );
      };
   
      ds=data.frame(e=e,case=case,which.set=which.set);
      beta.hat = c(beta.hat,clogit(case~e+strata(which.set) )$coefficients[1])
   };
   Q1.beta.hat = c(Q1.beta.hat, summary(beta.hat)[2] ) ; 
   median.beta.hat = c(median.beta.hat, median(beta.hat) ) ; 
   Q3.beta.hat = c(Q3.beta.hat, summary(beta.hat)[5] ) ;
}

min.beta.hat = .5; max.beta.hat = 1.5;

plot( seq(sizes[1],sizes[n.sizes],(sizes[n.sizes]-sizes[1])/(n.sizes-1) ), 
         seq( min.beta.hat, max.beta.hat ,( max.beta.hat - min.beta.hat)/(n.sizes-1) ), 
     type = "n") 
     
points(sizes,  Q1.beta.hat)
points(sizes, median.beta.hat,col="red")
points(sizes, Q3.beta.hat)