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[j + 1] - obs.t[i] - eps) * fail[i] } } # Model for(j in 1:T) { for(i in 1:N) { dN[i, j] ~ dpois(Idt[i, j]) # Likelihood Idt[i, j] <- Y[i, j] * exp(beta * Z[i]) * dL0[j] # Intensity } 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)); } c <- 0.001 r <- 0.1 for (j in 1 : T) { dL0.star[j] <- r * (t[j + 1] - t[j]) } beta ~ dnorm(0.0,0.000001) }