model
   {
   # Set up data
       for(i in 1 : N) {
          for(j in 1 : T) {
    # risk set = 1 if obs.t >= t
             Y[i, j] <- step(obs.t[i] - t[j] + eps) 
    # counting process jump = 1 if obs.t in [ t[j], t[j1] )
    # i.e. if t[j] <= obs.t < t[j1]
             dN[i, j] <- Y[i, j ] *step(t[j1] - obs.t[i] - eps)*fail[i] 
          }
      }
   # Model 
      for(j in 1 : T) {
         for(i in 1 : N) {
            dN[i, j] ~ dpois(Idt[i, j]) 
            Idt[i, j] <- Y[i, j] * exp(beta * Z[i]+b[pair[i]]) * dL0[j] 
         } 
         dL0[j] ~ dgamma(mu[j], c)
         mu[j] <- dL0.star[j] * c # prior mean hazard
   # Survivor function = exp(-Integral{l0(u)du})^exp(beta * z) 
         S.treat[j] <- pow(exp(-sum(dL0[1 : j])), exp(beta * -0.5))
         S.placebo[j] <- pow(exp(-sum(dL0[1 : j])), exp(beta * 0.5))   
      }
      for(k in 1 : Npairs) {
         b[k] ~ dnorm(0.0, tau);
      }
      tau ~ dgamma(0.001, 0.001)
      sigma <- sqrt(1 / tau)
      c <- 0.001 r <- 0.1 
      for (j in 1 : T) { 
         dL0.star[j] <- r * (t[j1]-t[j]) 
      } 
      beta ~ dnorm(0.0,0.000001) 
   }