# buyseph0204
# 02/01/04
# fits Cox proportional hazards model from 
#        Cowles, 2004, "Bayesian methods for evaluating surrogate markers,"
#        Technical report, Department of Statistics and Actuarial Science,
#        The University of Iowa
# uses data augmentation
# See note below data listing regarding the calculation of the variables "newobs.t" and "new.cent"
# N -- number of subjects
# T -- number of distinct failure times

model 
{

logsqrt2pi <- log( 2 * 3.14159265) / 2
log2 <- log(2.0)

for(i in 1:N) {
	for (j in 1:n[i]) {
		lambda[i,j ] ~ dexp(1)
		logmu[i,j] <-  beta0[j] + beta[1] * trt[i] + beta[2] * strat[i]
		mu[i,j] <-  exp( logmu[i,j])
	
		z[i,j] <- newobs.t[i,j] / (mu[i,j] * sqrt(2 * lambda[i,j])  )						# z is standard half normal
		logjacobian[i,j] <-  - logmu[i,j] - log( 2 * lambda[i,j]) / 2
		llik1[i,j] <- logjacobian[i,j] + log2 - logsqrt2pi  - pow( z[i,j], 2) / 2
		zeroes[i,j] <- 0
		psi[i,j] <- -llik1[i,j]  + 1 # minus log likelihood
		zeroes[i,j] ~ dpois( psi[i,j])
		newobs.t[i,j] ~ dunif(newcen.t[i,j], 5)
	}	
	s[i] ~ dnorm( mumark[i], deltastar1 )
	mumark[i] <- betas[1] + betas[2] * trt[i] + betas[3] * strat[i] + deltastar2 * sum(z[i,1:n[i]])	# E(s|z)

}	


for (k in 1:T) {
	beta0[k] ~ dflat()
}
for(k in 1:3) {
	betas[k] ~ dflat()
}

for(k in 1:2) {
	beta[k] ~ dflat()
}

	# The next lines are equivalent to Tau[1:2,1:2] ~ dwish( R[1:2,1:2], nu)

deltastar3 <- 1 
as <- nu / 2
bs1 <- R[1,1] - pow(R[1,2], 2) / R[2,2]
bs <- bs1 / 2
deltastar1 ~ dgamma(as, bs)
priorprec <- R[2,2] * deltastar1
priormean <- R[1,2] / R[2,2]
deltastar2 ~ dnorm(priormean, priorprec)

# transformations to quantities of interest
Sigma[2,2] <- 1 / deltastar3
Sigma[1,2] <- deltastar2 * Sigma[2,2]
Sigma[1,1] <- 1/deltastar1 + pow(Sigma[1,2],2) / Sigma[2,2]

gammaz <- Sigma[1,2] / sqrt(Sigma[1,1] * Sigma[2,2])

RE3 <-beta[1] / (betas[2] / sqrt(Sigma[1,1])  )      
}


#inits
list(betas = c(0.2,0.4,0.2), beta = c(-0.90,-1.3),  beta0 = c(-5,-5,-5,-5,-5,-5,-5,-5,-5,-5), deltastar1 = 1.414, deltastar2 = 0.5),
lambda = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1)
) 	


#data

list(  R = structure(.Data = c(1,.5,.5,1), .Dim = c(2,2)), nu = 4)

#Also requires data file C:\my documents\actg320\buysephrna.txt
# The data file contains the following variables:
#      s[]      -- values of surrogate marker; one value per subject; missing values not permitted
#      trt[]     --  0/1 indicator variable for treatment group
#      strat[] --  0/1 indicator variable for treatment group
#      newobs.t[,] -- N  by T matrix of values
#                newobs.t[i,j] = NA:  subject i was still in risk set at distinct failure time j
#                newobs.t[i,j] = 1:     subject i failed at distinct failure time j
#                newobs.t[i,j] = 0:     subject i is no longer in risk set at distinct failure time j
#     newcent.t[,] -- N by T matrix of values
#                newcent.t[i,j] = 1:     subject i is still in risk set at distinct failure time j
#                newcent.t[i,j] = 0:     subject i failed at distinct failure time j or earlier