# jh 2012.11.19

# Woburn Study Table 2 -- leukemia cases and peers

dx.year=1900+c(66,69,69,72,72, 73,74,75,75.5,76, 76,78,79,80,81,82,83)
birth.year.case=1900+c(59,57,64,65,68,70,65,64,75,63,72,63,69,66,68,79,74)
age.dx.case = dx.year - birth.year.case
n.in.riskset=c(218,290,265,182,183,170,213,239,115,219,132,219,164,199,187,154,84);
n.risksets = length(n.in.riskset); n.risksets

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);
controls = n.in.riskset - 1

require(survival)

# start of function definition


sample.and.plot.results = function(sizes) {
	
	n.ctls.from.each.set = sizes;
	SIMS=length(n.ctls.from.each.set)

   par(mfrow=c(1,SIMS), plt=c(0.1,0.99,.1,.99) );
   for (sim in 1:SIMS) {
	   n.exposed.ctls=c() ; n.unexposed.ctls=c(); 
      actual.n.from.each = c(); pe=c();
      for (i in 1:17) {
         all.ctls=c( rep(1,controls.exposed[i]), rep(0,controls.not.exposed[i]) )
         n.all=length(all.ctls); n.sampled = min(n.ctls.from.each.set[sim],n.all);
         exposed.ctls.in.sample = sum( sample( all.ctls, n.sampled ) );
         n.exposed.ctls  = c(n.exposed.ctls, exposed.ctls.in.sample);
         n.unexposed.ctls= c(n.unexposed.ctls, n.sampled-exposed.ctls.in.sample);
         actual.n.from.each = c(actual.n.from.each,n.sampled);
         pe = c(pe,exposed.ctls.in.sample/n.sampled)
      }
            
      # MH ratio
      
      n = actual.n.from.each+1;
      mhNum=sum( case.exposed * n.unexposed.ctls   / n ) ;
      mhDen=sum( case.not.exposed * n.exposed.ctls / n ) ;
      logHR.MH = log(mhNum/mhDen)
            
      # set up data for coxph
      
      age=age.dx.case;case=rep(1,17);
      exposed=case.exposed;freq=rep(1,17);riskset.no=1:17;
      ds.cases       = cbind(age,case,exposed,freq,riskset.no);
      age=age.dx.case;case=rep(0,17);exposed=rep(1,17);
      freq=n.exposed.ctls;riskset.no=1:17;
      ds.exp.ctls    = cbind(age,case,exposed,freq,riskset.no);
      age=age.dx.case;case=rep(0,17);exposed=rep(0,17);
      freq=n.unexposed.ctls;riskset.no=1:17;
      ds.nonexp.ctls = cbind(age,case,exposed,freq,riskset.no); 
      temp = rbind(ds.cases,ds.exp.ctls,ds.nonexp.ctls);
      ds=data.frame(temp[,1],temp[,2],temp[,3],temp[,4],temp[,5])
      names(ds)=c("age","case","exposed","freq","riskset.no")
   
      # fit ph model matching each riskset of age.at.dx and year of birth
      
      # remove observations with 0 frequency, then fit model
      ds = ds[ds$freq>0,];
      fit = coxph(Surv(age,case) ~ exposed + strata(riskset.no), 
        weights=freq,data=ds)
      alpha.hat = fit$coefficients ; me = 1.96*sqrt(fit$var)
      HR.hat=exp(alpha.hat) ; 
      alpha.lo = alpha.hat-me; alpha.hi = alpha.hat+me; ;
         
      plot(c(1964.5,1986),c(0,20),type="n",xlab="Year",ylab="Age",ylim=c(0,16)) 
      dy=2; dx=0.4;
      for (i in 1:17) {
	      segments(birth.year.case[i],0 ,
                 dx.year[i], dx.year[i]-birth.year.case[i],lty=3,col="grey")
        #rect(dx.year[i],dx.year[i]-birth.year.case[i],
        #     dx.year[i]+dy,dx+dx.year[i]-birth.year.case[i],
        #     col=c("green","red")[case.exposed[i]+1],border=NA )
        if(case.exposed[i]==0) points(
             c(dx.year[i]+1.00*dy),c(dx/3+dx.year[i]-birth.year.case[i]),
             col="green",pch=19,cex=0.75)
             
        if(case.exposed[i]==1) points(
             c(dx.year[i]-0.00*dy),c(dx/3+dx.year[i]-birth.year.case[i]),
             col="red",pch=19,cex=0.75)
             
        if(pe[i]>0)
             rect(dx.year[i],-0.05+dx.year[i]-birth.year.case[i],
             dx.year[i]+dy*pe[i],-0.05-dx+dx.year[i]-birth.year.case[i],
             col="red",border=NA );
        if(pe[i]<1)
           rect(dx.year[i]+dy*pe[i],-0.05+dx.year[i]-birth.year.case[i],
                dx.year[i]+dy      ,-0.05-dx+dx.year[i]-birth.year.case[i],
                col="green",border=NA )     
        text(dx.year[i]+1.65*dy, -0.05-dx+dx.year[i]-birth.year.case[i],
             toString(actual.n.from.each[i]),cex=0.8,adj=c(1,0))
        text(dx.year[i]+1.65*dy,dx.year[i]-birth.year.case[i],"1",
             cex=0.8,adj=c(1,0))
      }
      y.0=12; dy=0.8;
      text(1966,y.0+ 4.4*dy,"log HR",adj=c(1,0),cex=0.7);
      text(1967,y.0+ 4.4*dy,"Cox",adj=c(0.5,0),cex=0.6)
      text(1968.25,y.0+ 4.4*dy,"MH",adj=c(0.5,0),cex=0.6)
      
      for(alpha in seq(0,4,1)) {
         segments(1966,y.0+alpha*dy, 1966.5,y.0+alpha*dy);
   	    text(1966-0.2,y.0+ alpha*dy,toString(alpha),adj=c(1,0.5),cex=0.6)
      }
      points(c(1967),y.0+c(alpha.hat)*dy,pch=20)
      segments(1967,y.0+ alpha.lo*dy, 1967,y.0+ alpha.hi*dy)
      points(c(1968),y.0+c(logHR.MH)*dy,pch=20)
   }

}  # end of function definition


# n.ctls.from.each.set = c(1,4,10);

# n.ctls.from.each.set = c(1,4,10,300);

n.ctls.from.each.set = c(10,300); # will use min(these no.'s , available no.'s)
    
# CALL the function

 sample.and.plot.results(n.ctls.from.each.set)