# Woburn data from Table 2 of Lagakos Wesson and Zelen JASA 1986
# jh 2008.02.20

dx=1970+c(-4,-1,-1,2,2,3,4,5,5,6,6,8,9,10,11,12,13)
born=1900+c(59,57,64,65,68,70,65,64,75,63,72,63,69,66,68,79,74)

case.cum.exposure=c(1.26,0,.75,4.3,2.76,.94,0,0,0,.37,0,7.88,2.41,0,0,.39,0)
n=c(218,290,265,182,183,170,213,239,115,219,132,219,164,199,187,154,84)
e.cum.exposure=c(.31,.34,.17,.90,.58,.20,.56,.99,.09,1.18,.24,1.41,.73,1.38,1.14,.08,.25)
v.cum.exposure=c(.26,.36,.10,2.23,.88,.20,1.04,2.78,.03,3.87,.32,6.23,2.56,6,4.2,.02,.45)
p.riskset.exposed=c(.33,.26,.25,.36,.32,.19,.29,.38,.25,.40,.18,.40,.31,.39,.35,.23,.23)

case.exposed = case.cum.exposure > 0

n.exposed = round(p.riskset.exposed * n,0)
n.not.exposed = n - n.exposed

#####  ML fitting

## directly... by eye

log.lik = function(hr){
 	       sum( log( case.exposed*hr + (1-case.exposed) ) - 
 	            log( n.exposed*hr    +  n.not.exposed*1  )
 	          ) 
 	}

hr=2^seq(-0.75,3.5,0.25);
plot(hr,sapply(hr,log.lik))

log.hr=seq(-1,3,0.1);
plot(log.hr,sapply(exp(log.hr),log.lik))

log.lik(1)

## estimating equation, and trial and error

U = function(beta){ 
	            hr= exp(beta);
	            return(
	             sum(case.exposed) - 
 	             sum( (1*hr*n.exposed + 0*1*n.not.exposed)/(hr*n.exposed+1*n.not.exposed)  )
 	                  )
 	                     }

I = function(beta){ 
	            hr= exp(beta);
	            return(
	             sum( (1*1*hr*n.exposed + 0*0*1*n.not.exposed) /
	                  (    hr*n.exposed +     1*n.not.exposed) 
	                  -
	                  ((1*hr*n.exposed + 0*1*n.not.exposed)/(hr*n.exposed+1*n.not.exposed))^2 
	                )
	                   )
 	                     }

beta.range=seq(-1,3,.1);

plot(beta.range,sapply(beta.range,U))

## newton-raphson iteration

beta.prev = 0;
for (iteration in 1:6){
	if (iteration==1) print("  beta       hr      log.lik       U          I   ");
	print(round(c(beta.prev, exp(beta.prev), log.lik(exp(beta.prev)), U(beta.prev), I(beta.prev) ),5)
	             );
	beta.prev= beta.prev + U(beta.prev)/I(beta.prev);
	}
	
	
# X^2 statistic

U(0)^2 / I(0)