########################################################################
# Survival models:
log(time[i])=b[1]
+ b[2]*treat[i] + b[3]*age[i]+sigma*eps[i], eps[i]~Gumbel
########################################################################
model
{
for(i in 1:n)
{
cens[i] <- (1-status[i])*time[i] # If censored, status[i] = 0.
time[i] ~ dweib(alpha,
lambda[i])I(cens[i],)
lambda[i] <- exp( -mu[i]/sigma )
mu[i] <- b[1] + b[2]*treat[i] + b[3]*age[i]
}
alpha <-
1/sigma
b[1] ~ dnorm(0, 0.000001)
b[2] ~ dnorm(0, 0.000001)
b[3] ~ dnorm(0, 0.000001)
sigma ~ dgamma(0.0001, 0.0001)
#
Estimate Survival function at t=10,20,...,2000.
for(i in 1:200)
{
t[i] <- 10*i
z[i] <- (log(t[i]) - b[1] - b[2] -b[3]*60 )/sigma
S[i] <- exp(-exp( z[i])
)
}
}
# Initial values
list( b=c(0, 0, 0), sigma = 1)
# Data 1
list(n = 26)
# Data 2
time[] status[] treat[] age[]
156
1
1 66
1040
0
1 38
59
1
1 72
421
0
2 53
329
1 1
43
769
0
2 59
365
1
2 64
770
0
2 57
1227
0
2 59
268
1
1 74
475
1
2 59
1129
0
2 53
464
1
2 56
1206
0
2 44
638
1
1 56
563
1
2 55
1106
0
1 44
431
1
1 50
855
0
1 43
803
0
1 39
115
1
1 74
744 0 2 50
477
0
1 64
448
0
1 56
353
1
2 63
377
0
2 58
END