########################################################################

# Section 6.10 Stanford Heart Transplant Data (using 0's trick)

########################################################################

model

{

            for(i in 1:n)

            {

                       # Likelihood

                       llk_nt_0[i] <- log(p)+p*log(lambda)-(p+1)*log(lambda+survtime[i])

                       llk_nt_1[i] <- p*log(lambda)-p*log(lambda+survtime[i])

                       llk_t_0[i]  <- log(p)+log(tau)+p*log(lambda)

                                      -(p+1)*log(lambda+timetotransplant[i]+tau*survtime[i])

                       llk_t_1[i]  <- p*log(lambda)

                                      -p*log(lambda+timetotransplant[i]+tau*survtime[i])

                       loglk[i]    <- (1-transplant[i])*(1-state[i])*llk_nt_0[i]

                                     +(1-transplant[i])* state[i]   *llk_nt_1[i]

                                     +   transplant[i] *(1-state[i])*llk_t_0[i]

                                     +   transplant[i] * state[i]   *llk_t_1[i]

                       mu[i] <-  -loglk[i]+10000   # Add some large constant to make lam[i] > 0

                       zeros[i] <- 0

                       zeros[i] ~ dpois(mu[i])  

                                                         

            }

           

tau    <- exp(ltau)

            lambda <- exp(llambda)

            p      <- exp(lp)

 

            ltau    ~dnorm(0, 0.000001)

            llambda ~dnorm(0, 0.000001)

            lp      ~dnorm(0, 0.000001)

                                  

            # Estimate Pr(Tau<= 1|y): Lower risk with a transplant operation?

            tau.le.1 <- step(1 - tau) 

     

     # Prediction of survival function

      for(t in 1:240)

            {

                       S[t] <- exp( p*log(lambda) - p*log(lambda + t) )

                       time[t] <- t

            }

}                                                                                

# Initial values. 

            list( ltau = 0, llambda=0, lp=0 )

# Data

            list(n = 82)

# Data

survtime[]       transplant[]     timetotransplant[]      state[]

49        0          0          0

5          0          0          0

17        0          0          0

2          0          0          0

39        0          0          0

84        0          0          0

7          0          0          0

0          0          0          0

35        0          0          0

36        0          0          0

1400    0          0          1

5          0          0          0

34        0          0          0

15        0          0          0

11        0          0          0

2          0          0          0

1          0          0          0

39        0          0          0

8          0          0          0

101      0          0          0

2          0          0          0

148      0          0          0

1          0          0          0

68        0          0          0

31        0          0          0

1          0          0          0

20        0          0          0

118      0          0          1

91        0          0          1

427      0          0          1

15        1          0          0

3          1          35        0

624      1          50        0

46        1          11        0

127      1          25        0

61        1          16        0

1350    1          36        0

312      1          27        0

24        1          19        0

10        1          17        0

1024    1          7          0

39        1          11        0

730      1          2          0

136      1          82        0

1379    1          24        1

1          1          70        0

836      1          15        0

60        1          16        0

1140    1          50        1

1153    1          22        1

54        1          45        0

47        1          18        0

0          1          4          0

43        1          1          0

971      1          40        1

868      1          57        1

44        1          0          0

780      1          1          1

51        1          20        0

710      1          35        1

663      1          82        1

253      1          31        0

147      1          40        0

51        1          9          0

479      1          66        1

322      1          20        0

442      1          77        1

65        1          2          0

419      1          26        1

362      1          32        1

64        1          13        0

228      1          56        0

65        1          2          0

264      1          9          1

25        1          4          0

193      1          30        1

196      1          3          1

63        1          26        0

12        1          4          0

103      1          45        1

60        1          25        1

43        1          5          1

END