model cox1;

  const 
   N = 42, # number of patients 
   T = 18,  # number of unique failure times
   eps = 0.000001; # used to guard against numerical imprecision in step function


  var
   beta1,   # regression coefficient 
   D, MAR, dsq[N],
  failtime[N],   # observed remission   or censoring time for each patient
  t[T+1],     # unique failure  times + maximum censoring time
  dN[N,T],    # counting process increment
  Y[N,T],     # 1 = if subject  observed; 0 = not observed
  Idt[N,T],   # intensity process
  trt[N],     # treatment covariate
  dL0[T],     # increment in unknown cumulative baseline hazard function
  beta0[T],   # log(increment in unknown cumulative baseline hazard function)
  dL0.star[T], # prior specification for cumulative baseline hazard function
  c0,           # degree of confidence in prior specification for dL0
  mu[T],        # location parmater for gamma process prior ( = c0 * dL0.star)  
  r0,           # prior specification for failure rate
  failcens[N],      # failure  = 1, censored  = 0
  S.trt[N],    # survivor function for treatment group
  S.placebo[N],   # survivor function for placebo  group
  M[N],           # martingale residual
  d[N],           # deviance residual
  pm.post[2], likes[2], llike[N,T,2], sumlike, Idtm[N,T,2], g[2];
  

  data failtime, failcens, trt in "leuk.dat", 
  t in "failtimes.dat";
  inits in "leuk.in";


  {
  # priors for regression coefficients

  beta1 ~ dnorm(0,.0000001);

  # set up data

  for(i in 1:N) {
   for(j in 1:T) {

   # risk set = 1 if failtime >= t
   Y[i,j] <- step(failtime[i] - t[j] + eps);
   
   # counting process jump = 1 if failtime in [ t[j], t[j+1] )
   #  i.e., if t[j] <= failtime < t[j+1] 

   dN[i,j] <- Y[i,j]*step(t[j+1] - failtime[i] - eps)*failcens[i];
   }
   }

 # Model


   for(j in 1:T) {

  # beta0[j] ~ dnorm(0,.001);  # include this when using the Poisson trick

   for (i in 1:N) {

    dN[i,j] ~ dpois(Idt[i,j]);  # likelihood
   Idt[i,j] <- Y[i,j]*exp(beta1*trt[i])*dL0[j]; 
      # Intensity without Poisson trick 

 # Try Poisson Trick - independent log-normal hazard increments 
 #                     enables dL0, c0, r0, mu to be dropped from the model
 # Idt[i,j] <- Y[i,j]*exp(beta1*trt[i]+ beta0[j]);
                               # intensity function based on Poisson trick
  
 }

   dL0[j] ~ dgamma(mu[j],c0)   # gamma process prior on the cumulative baseline hazard rate
   mu[j] <- dL0.star[j] * c0;  # prior mean cumulative hazard

 # Survivor function = exp(- Integral{l0(u) du})^{exp(x'beta)}

S.trt[j] <- pow(exp(-sum(dL0[1:j])), exp(beta1* -0.5));
S.placebo[j] <- pow(exp(-sum(dL0[1:j])), exp(beta1*0.5))
  

 }

c0 <- .001;
r0 <- 0.1;
for(j in 1:T) {
 dL0.star[j] <- r0 * (t[j+1] -t[j])
 }

for (i in 1:N) {
M[i] <- sum(dN[i,]) - sum(Idt[i,]);
d[i] <- (2*step(M[i]) - 1)*pow(2*(-M[i] - sum(dN[i,])*log(sum(dN[i,]) - M[i])),.5);
dsq[i] <- pow(d[i],2)

}

D <- sum(dsq[]);
MAR <- sum(M[]);


g[1] <- 0;
g[2] <- 1;


for(i in 1:N) {
   for(j in 1:T) {
    for (mod in 1:2) {
    Idtm[i,j,mod] <-  Y[i,j]*exp(g[mod]*beta1*trt[i])*dL0[j];
    llike[i,j,mod] <- dN[i,j]*log(Idtm[i,j,mod]) - Idtm[i,j,mod];
}
}
}

for (mod in 1:2){  
 likes[mod]<- exp(sum(llike[,,mod]));
 }
 sumlike<- sum(likes[]);
 for ( mod in 1:2){
 pm.post[mod]<- likes[mod]/sumlike;
 }


 }