#jointaft.bug # exponential survival model # random effects + random walk model for rna and cd4 # ability to compute AUC for rna # computes approximations to betaR1 and propcap # Kate Cowles # nint -- number of intervals for numeric integration (32) model{ for (i in 1:N) { trt2[i] <- trt[i] + 1 trt0[i] <- 1 - trt[i] # trt[i] = 1 if trt grp 1, 0 otherwise # trt0 is the reverse strat0[i] <- 1 - strat[i] grp00[i] <- trt0[i] * strat0[i] grp01[i] <- trt0[i] * strat[i] grp10[i] <- trt[i] * strat0[i] grp11[i] <- trt[i] * strat[i] for (j in 1:(nint+1)) { # compute endpoints of intervals for numeric integration t[i,j] <- obs.t[i] * (j - 1) / nint } } # The next part is for computing PTE n00 <- inprod(grp00[], obs.t[]) n10 <- inprod(grp10[], obs.t[]) n01 <- inprod(grp01[], obs.t[]) n11 <- inprod(grp11[], obs.t[]) mean00 <- inprod(grp00[1:N], lambda2[1:N]) / n00 mean01 <- inprod(grp01[1:N], lambda2[1:N]) / n01 mean10 <- inprod(grp10[1:N], lambda2[1:N]) / n10 mean11 <- inprod(grp11[1:N], lambda2[1:N]) / n11 betaR1 <- (log(mean10) - log(mean00) + log(mean11) - log(mean01) ) / 2 propcap <- 1 - beta1 / betaR1 for (i in 1:N) { # Use "zeroes trick" to fit AFT model zeroes[i] <- 0 zeroes[i] ~ dexp(lambda[i]) lambda2[i] <- exp( beta0 + beta1 * trt[i] +beta2 * strat[i] + log(t0[i]) ) lambda[i] <- exp( fail[i] * (beta0 + beta1 * trt[i] + beta2 * strat[i] + beta3 * predrna[i, (nint+1)] + beta4 * predcd4[i,(nint+1)] ) - lambda2[i] ) for (j in 1:(nint + 1)) { # use the following one for instantaneous value of rna # predrna[i, j] <- alpha[i,1] + alpha[i,2] * t[i,j] + # ((t[i, j] - x[ leftend[i,j] ]) * w[i, leftend[i, j] + 1,1] + # (x[ leftend[i, j] + 1] - t[i, j]) * w[i, leftend[i, j] ,1]) / # (x[ leftend[i, j] + 1 ] - x[leftend[i, j] ] ) # use the following one instead for area under curve of rna # Following is integrated area under random walk for (k in 1:(leftend[i,j] - 1) ) { summand[i,j,k] <- (w[i,k+1,1] + w[i, k, 1]) * (x[k+1] - x[k]) / 2 } summand[i,j,leftend[i,j]] <- ( t[i,j] - x[leftend[i,j]]) * (w[i, leftend[i,j], 1] + ( t[i,j] - x[leftend[i,j] ]) * (w[i, leftend[i,j]+1, 1] - w[ i, leftend[i,j], 1]) / (2 * (x[leftend[i,j]+1] - x[leftend[i,j]])) ) randarea[i,j] <- sum(summand[i,j,1:leftend[i,j]]) predrna[i,j] <- alpha[i,1] * t[i,j] + mu[trt2[i],3] * pow( t[i,j], 2) / 2 + randarea[i,j] predcd4[i, j] <- alpha[i,2] + mu[trt2[i],4] * t[i, j] + ((t[i, j] - x[ leftend[i, j] ]) * w[i, leftend[i, j] + 1,2] + (x[ leftend[i, j] + 1] - t[i,j]) * w[i, leftend[i, j] ,2]) / (x[ leftend[i, j] + 1 ] - x[leftend[i, j] ] ) func[i,j] <- exp(beta3 * predrna[i,j] + beta4 * predcd4[i,j]) ; } t0[i] <- obs.t[i] * inprod( func[i,], mult[1:(nint+1)]) / (3 * nint) for (j in 1:times) { mu1[i, j] <- alpha[i,1] + mu[trt2[i],3] * x[j] + w[i,j,1] rna[i, j] ~ dnorm(mu1[i,j], tau1) mu2[i,j] <- alpha[i,2] + mu[trt2[i],4] * x[j] + w[i,j,2] cd4[i,j] ~ dnorm(mu2[i,j], tau2) } # set up w's as random walk w[i,1,1] <- 0 w[i,1,2] <- 0 for (j in 2:times) { for (k in 1:2) { w[i,j,k] <- w[i,j-1,k] + sqrt(x[j] - x[j-1]) * incr[i,j,k] } incr[i,j,1:2] ~ dmnorm(zero2[1:2], precw[1:2,1:2]) } for (k in 1:2) { # needed in computing predrna and predcd4 w[i,times + 1, k] <- w[i,times,k] } alpha[i,1:2] ~ dmnorm(mu[ trt2[i], 1:2], precalph[1:2,1:2]) } # Priors mu[1,1:4] ~dmnorm(zero4[1:4], precmu[1:4,1:4]) mu[2,1:4] ~dmnorm(zero4[1:4], precmu[1:4,1:4]) precw[1:2,1:2] ~ dwish( omegaw[1:2,1:2], 4 ) precalph[1:2,1:2] ~ dwish( omegalph[1:2,1:2], 4) tau1 ~ dgamma(0.1, 0.1) tau2 ~ dgamma(0.1, 0.1) beta0 ~ dnorm(0.0, 0.0001) beta1 ~ dnorm(0.0, 0.0001) beta2 ~ dnorm(0.0, 0.0001) beta3 ~ dnorm(0.0, 0.0001) beta4 ~ dnorm(0.0, 0.0001) # Deterministic transformations sigw[1,1] <- inverse(precw[,],1,1) sigw[1,2] <- inverse(precw[,],1,2) sigw[2,2] <- inverse(precw[,],2,2) sigalph[1,1] <- inverse(precalph[,],1,1) sigalph[1,2] <- inverse(precalph[,],1,2) sigalph[2,2] <- inverse(precalph[,],2,2) sig1 <- 1/tau1 sig2 <- 1/tau2 } # Data files # data for full dataset with nint = 32 # requires two additional files as indicated below list(N=197, x = c(0,4,8,24,40,55), times = 5, zero4 = c(0,0,0,0), precmu = structure(.Data = c(0.001,0,0,0, 0, 0.001,0,0, 0,0,0.001,0, 0,0,0,0.001), .Dim = c(4,4)), omegalph = structure( .Data = c(1,0,0,1), .Dim = c(2,2)), nint=32, mult = c(1, 4,2,4,2,4,2,4,2,4,2,4,2,4,2,4,2,4,2,4,2,4,2,4,2,4,2,4,2,4,2,4,1 ), zero2 = c(0,0), omegaw = structure( .Data = c(0.25,0,0,0.25), .Dim = c(2,2))) # Note: x[times + 1] must be greater than t[T + 1] # requires additional data file in tabular format trt[] strat[] rna[,1] rna[,2] rna[,3] rna[,4] rna[,5] cd4[,1] cd4[,2] c d4[,3] cd4[,4] cd4[,5] obs.t[] fail[] 0 1 4.2479 3.27323 4.0566 3.9829 NA 3.705 3.511 3.653 3.488 NA 32.8571 0 0 0 5.56951 5.10036 4.96781 5.41695 5.01041 2.088 2.561 2.432 1.627 NA 43.8571 0 0 0 4.96314 5.4852 4.60326 5.29003 NA 1.968 1.968 2.088 1.968 NA 38.8571 0 . . . # requires third data file of leftend list(leftend=structure(.Data=c(1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3 , 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, . . . 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5), .Dim = c(197,33))) # Initial values list(mu = structure(.Data = c(4.8, 2.7,-0.4, 0.0, 4.8, 2.8, -0.09, 0.03), .Dim = c( 2,4)), tau1 = 5, tau2 = 50, precalph = structure(.Data = c(4,0,0,1.5), .Dim = c(2,2)), beta0 = -7.5, beta1 = 0.5, beta2 = 0.9, beta3 = 0.1, beta4 = -1.3, precw = structure( .Data = c(5,0 ,0,30), .Dim = c(2,2))) list(mu = structure(.Data = c(5,3,0.0, 0.0, 5, 3,0.0, 0.0), .Dim = c( 2,4)), tau1 = 10, tau2 = 10, precalph = structure( .Data = c(10,0,0,10), .Dim = c(2,2)), beta0 = -10.0, beta1 =-0.25, beta2 = 0, beta3 = 0.0, beta4 = -1.7, precw = structure( .Data = c(10,0 ,0,10), .Dim = c(2,2))) list(mu = structure(.Data = c(6,4,-0.01, -0.01, 6,4,-0.01, -0.01), .Dim = c( 2,4)), tau1 = 1, tau2 = 1, precalph = structure( .Data = c(1,0,0,1), .Dim = c(2,2)), beta0 = -6.0, beta1 = 1.25, beta2 = 1.7, beta3 = 0.2, beta4 = -0.9, precw = structure( .Data = c(1,0 ,0,1), .Dim = c(2,2)))