CHAPTER 12:  Time to Event Data
One and two group methods

Exercise 12.4

model{
 lambda ~ dgamma(a,b)
 alpha ~ dlnorm(c,d)
 # alpha ~ dgamma(c,d)
 # alpha ~ dunif(0,u)
 med <- pow(0.69/lambda,1/alpha)
 surv <- exp(-lambda*pow(20,alpha))
}
list(a=1 , b=1, c=1, d=1) # Insert the values from part (a) here
Note:  The value 1 for a,b,c,d is NOT what you get when you do Exercise 12.4. 
           I used them just to make sure the code ran.

	 node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	alpha	4.425	5.549	0.08169	0.3636	2.706	19.34	1	5000
	lambda	1.002	1.024	0.01349	0.02404	0.6789	3.691	1	5000
	med	4.37E+15	3.09E+17	4.383E+15	0.2066	1.004	25.9	1	5000
	surv	0.0319	0.1323	0.001793	0.0	0.0	0.5555	1	5000

Exercise 12.7

model{
 for(i in 1:21){
   n[i] ~ dbin(p[i],r[i])
   r[i] <- n[i]
   p[i] <- phi((log(u[i]) - mu)/sigma) - phi((log(l[i]) - mu)/sigma)
 }
 mu ~ dnorm(0,0.001)
 tau ~ dgamma(0.001,0.001)I(0.001,)
 sigma <- 1/sqrt(tau)
# Inference
 med <- exp(mu)
 ninetypct <- med*exp(1.28*sigma)
 surv <- 1-phi((log(16) - mu)/sigma)
 junk <- age[1]
}
list(mu=0, tau=1)

DATA

age[ ] n[ ] l[ ] u[ ]
10 4 10 10.99
11 12 11 11.99
12 19 12 12.99
12 2 12.5 100
13 24 13 13.99
13 15 13.5 100
13 1 0 13.5
14 20 14 14.99
14 24 14.5 100
14 2 0 14.5
15 13 13 13.99
15 18 15.5 100
15 3 0 15.5
16 3 16 16.99
16 14 16.5 100
16 2 0 16.5
17 1 17 17.5
17 6 17.5 100
17 3 0 17.5
18 1 0 18.5
18 4 18.5 100
END

LEUKEMIA DATA.  EXPONENTIAL SAMPLING MODEL
# Code for Example 12.3.3

model
{
 for(i in 1:n[1]){t1[i] ~ dexp(theta1)}
 for(i in 1:n[2]){t2[i] ~ dexp(theta2)}
 theta1 ~ dgamma(a1,b1)
 theta2 ~ dgamma(a2,b2)
 median1 <- log(2)/theta1
 median2 <- log(2)/theta2
 relmedian <- theta2/theta1
 S[1] <- exp(-24*theta1)
 S[2] <- exp(-24*theta2)
 Sdiff <- S[1]-S[2]
}
list(n=c(17,16),
a1=1.53,b1=26.4,a2=2.31,b2=37.95,
t1=c(65,156,100,134,16,108,121,4,39,
          143,56,26,22,1,1,5,65),
t2=c(56,65,17,7,16,22,3,4,2,3,8,4,3,30,4,43))
list(theta1=0.05,theta2=0.02)

RESULTS

	 node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	S[1]	0.6675	0.06226	2.913E-4	0.5406	0.6697	0.7831	1	50000
	S[2]	0.2716	0.08202	3.599E-4	0.129	0.2656	0.4469	1	50000
	Sdiff	0.3959	0.1026	4.67E-4	0.1834	0.4003	0.5846	1	50000
	median1	43.03	10.56	0.04936	27.04	41.49	68.06	1	50000
	median2	13.03	3.234	0.01427	8.124	12.55	20.65	1	50000
	relmedian	3.495	1.202	0.005393	1.718	3.303	6.38	1	50000
	theta1	0.01703	0.003966	1.851E-5	0.01019	0.01671	0.02563	1	50000
	theta2	0.05632	0.01325	5.768E-5	0.03357	0.05524	0.08533	1	50000

Sensitivity analysis
model
{
 for(i in 1:n[1]){t1[i] ~ dexp(theta1)}
 for(i in 1:n[2]){t2[i] ~ dexp(theta2)}
 #theta1 ~ dgamma(0.001,0.001)
 #theta2 ~ dgamma(0.001,0.001)
 theta1 ~ dunif(0,1000)
 theta2 ~ dunif(0,1000)
 median1 <- log(2)/theta1
 median2 <- log(2)/theta2
 relmedian <- theta2/theta1
 S[1] <- exp(-24*theta1)
 S[2] <- exp(-24*theta2)
 Sdiff <- S[1]-S[2]
}
list(n=c(17,16),t1=c(65,156,100,134,16,108,121,4,39,143,56,26,22,1,1,5,65),
t2=c(56,65,17,7,16,22,3,4,2,3,8,4,3,30,4,43))
list(theta1=0.05,theta2=0.02)

Results

 theta1 ~ dgamma(0.001,0.001)
 theta2 ~ dgamma(0.001,0.001)
 	node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	S[1]	0.6842	0.06248	2.716E-4	0.5553	0.6865	0.8003	1	50000
	S[2]	0.2766	0.08798	4.229E-4	0.1248	0.27	0.4644	1	50000
	Sdiff	0.4076	0.1081	5.168E-4	0.1828	0.4127	0.6035	1	50000
	median1	46.06	11.91	0.05241	28.28	44.23	74.67	1	50000
	median2	13.26	3.549	0.01706	7.995	12.7	21.69	1	50000
	relmedian	3.708	1.362	0.006365	1.731	3.48	6.954	1	50000
	theta1	0.01599	0.003883	1.683E-5	0.009285	0.01567	0.02451	1	50000
	theta2	0.0558	0.01407	6.81E-5	0.03196	0.05456	0.0867	1	50000

	
	theta1 ~ dunif(0,1000)
   theta2 ~ dunif(0,1000)
	 node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	S[1]	0.6842	0.06248	2.723E-4	0.5554	0.6866	0.8003	1	50000
	S[2]	0.2766	0.08798	4.23E-4	0.1248	0.27	0.4644	1	50000
	Sdiff	0.4076	0.1081	5.175E-4	0.1828	0.4127	0.6035	1	50000
	median1	46.06	11.91	0.05253	28.29	44.24	74.68	1	50000
	median2	13.26	3.55	0.01707	7.994	12.7	21.69	1	50000
	relmedian	3.709	1.362	0.006373	1.731	3.48	6.957	1	50000
	theta1	0.01599	0.003882	1.687E-5	0.009282	0.01567	0.0245	1	50000
	theta2	0.05579	0.01407	6.809E-5	0.03196	0.05456	0.08672	1	50000

LEUKEMIA DATA.  WEIBULL SAMPLING MODEL
#  Section 12.3.2.  Gamma priors on the alpha's

model{
 for(i in 1:n[1]){ t1[i] ~ dweib(alpha[1],lambda[1]) }
 for(i in 1:n[2]){ t2[i] ~ dweib(alpha[2],lambda[2])  }
 lambda[1] ~ dgamma(a[1],b[1]) 
 alpha[1] ~ dgamma(c[1],d[1])
 lambda[2] ~ dgamma(a[2],b[2]) 
 alpha[2] ~ dgamma(c[2],d[2])
 med[1] <- pow(log(2)/lambda[1],1/alpha[1])
 med[2] <- pow(log(2)/lambda[2],1/alpha[2])
 relmedian <- med[1]/med[2]
 S[1] <- exp(-pow(24,alpha[1])*lambda[1])
 S[2] <- exp(-pow(24,alpha[2])*lambda[2])
}
list(n=c(17,16), a= c(0.001,0.001), b=c(0,001,0,001), c=c(0.001,0.001),d=c(0.001,0.001),
t1=c(65,156,100,134,16,108,121,4,39,143,56,26,22,1,1,5,65),
t2=c(56,65,17,7,16,22,3,4,2,3,8,4,3,30,4,43))
list(alpha=c(1,1), lambda=c(1,1))

RESULTS

 node	 mean	 sd	 MC error	2.5%	median	97.5%	start	sample
	S[1]	0.639	0.09415	8.326E-4	0.4431	0.6436	0.8082	10000	200001
	S[2]	0.2774	0.08829	2.495E-4	0.1265	0.2704	0.4679	10000	200001
	alpha[1]	0.8811	0.1819	0.002033	0.5574	0.8689	1.267	10000	200001
	alpha[2]	0.9238	0.1745	0.001326	0.6073	0.9156	1.288	10000	200001
	lambda[1]	0.03664	0.03117	3.083E-4	0.004304	0.02793	0.1202	10000	200001
	lambda[2]	0.08294	0.05187	3.604E-4	0.01794	0.07112	0.2148	10000	200001
	med[1]	41.78	14.04	0.08917	18.76	40.32	73.36	10000	200001
	med[2]	12.5	4.116	0.01954	5.858	12.04	21.84	10000	200001
	relmedian	3.737	1.935	0.0109	1.278	3.339	8.555	10000	200001

LEUKEMIA DATA.  WEIBULL SAMPLING MODEL
#  Section 12.3.2.  Log normal priors on the alpha's
# This what it is in the book

model{
 for(i in 1:n[1]){ t1[i] ~ dweib(alpha[1],lambda[1]) }
 for(i in 1:n[2]){ t2[i] ~ dweib(alpha[2],lambda[2])  }
 lambda[1] ~ dgamma(a[1],b[1]) 
 alpha[1] ~ dlnorm(c[1],d[1])
 lambda[2] ~ dgamma(a[2],b[2]) 
 alpha[2] ~ dlnorm(c[2],d[2])
 med[1] <- pow(log(2)/lambda[1],1/alpha[1])
 med[2] <- pow(log(2)/lambda[2],1/alpha[2])
 relmedian <- med[1]/med[2]
 S[1] <- exp(-pow(24,alpha[1])*lambda[1])
 S[2] <- exp(-pow(24,alpha[2])*lambda[2])
}
list(n=c(17,16), a= c(0.001,0.001), b=c(0,001,0,001), c=c(0,0),d=c(0.001,0.001),
t1=c(65,156,100,134,16,108,121,4,39,143,56,26,22,1,1,5,65),
t2=c(56,65,17,7,16,22,3,4,2,3,8,4,3,30,4,43))
list(alpha=c(1,1), lambda=c(1,1))

Exercise 12.12

model{
# Likelihood
 for(i in 1:24){ t[i] ~ dlnorm(mu[h[i]],tau[h[i]])I(c[i],) }
 for(i in 1:21){ tt[i] ~ dlnorm(mu[hh[i]], tau[hh[i]])I(cc[i],) }
# Prior
 for(i in 1:2){ 
   mu[i] ~ dnorm(0,.001)
   tau[i] ~ dgamma(0.001,0.001)
# Inference
 med[i] <- exp(mu[i])
 surv[i] <- 1- phi((log(120) - mu[i])/sigma[i])
 sigma[i] <- 1/sqrt(tau[i])
 } 
 relmed <- med[1]/med[2]
 survdiff <- surv[1] - surv[2]
}
list(mu=c(0,0), tau=c(1,1), tt=c(122,123,130,130,133,134,136,
141,143,148,151,152,153,154,156,162,164,165,182,189,144))

DATA

t[]	c[]	h[]
19	0	1
25	0	1
30	0	1
34	0	1
37	0	1
46	0	1
47	0	1
51	0	1
56	0	1
57	0	1
61	0	1
66	0	1
67	0	1
74	0	1
78	0	1
86	0	1
22	0	2
23	0	2
38	0	2
42	0	2
73	0	2
77	0	2
89	0	2
115	0	2
END

tt[]	cc[]	hh[]
NA	122	1
NA	123	1
NA	130	1
NA	130	1
NA	133	1
NA	134	1
NA	136	1
NA	141	1
NA	143	1
NA	148	1
NA	151	1
NA	152	1
NA	153	1
NA	154	1
NA	156	1
NA	162	1
NA	164	1
NA	165	1
NA	182	1
NA	189	1
NA	144	2
END

model{
# Likelihood
 for(i in 1:24){ t[i] ~  dlnorm(mu[h[i]],   tau)I(c[i],) }
 for(i in 1:21){ tt[i] ~ dlnorm(mu[hh[i]], tau)I(cc[i],) }
# Prior
 for(i in 1:2){ mu[i] ~ dnorm(0,.001) }
 tau ~ dgamma(0.001,0.001)
 sigma <- 1/sqrt(tau)
# Inference
 for(i in 1:2){
 med[i] <- exp(mu[i])
 surv[i] <- 1- phi((log(120) - mu[i])/sigma)
 }
 relmed <- med[1]/med[2]
 survdiff <- surv[1] - surv[2]
}
list(mu=c(0,0), tau=1, tt=c(122,123,130,130,133,134,136,
141,143,148,151,152,153,154,156,162,164,165,182,189,144))

DATA

t[]	c[]	h[]
19	0	1
25	0	1
30	0	1
34	0	1
37	0	1
46	0	1
47	0	1
51	0	1
56	0	1
57	0	1
61	0	1
66	0	1
67	0	1
74	0	1
78	0	1
86	0	1
22	0	2
23	0	2
38	0	2
42	0	2
73	0	2
77	0	2
89	0	2
115	0	2
END

tt[]	cc[]	hh[]
NA	122	1
NA	123	1
NA	130	1
NA	130	1
NA	133	1
NA	134	1
NA	136	1
NA	141	1
NA	143	1
NA	148	1
NA	151	1
NA	152	1
NA	153	1
NA	154	1
NA	156	1
NA	162	1
NA	164	1
NA	165	1
NA	182	1
NA	189	1
NA	144	2
END

Induced prior on beta based on independent log normal priors on medians mtilde

model{ 
  for(i in 1:r){
    mtilde[i] ~ dlnorm(mu0[i],tautilde0[i]) # mtilde[i] ~ LN
    mu0[i] <- log(mtilde0[i])  # Prior guess for tilde m_i is mtilde0[i]
    tautilde0[i] <- pow(1/sigmatilde0[i],2)
    lmtilde[i] <- log(mtilde[i])    
    beta[i] <- inprod(G[i,1:r] , lmtilde[1:r])  
  }
}
list(r=*, mtilde0 = c(*,*), sigmatilde0 = c(*,*))
# The matrix G = solve(Xtilde) must be obtained in R