right censored data

Dave
2012-01-10
2012-09-01
  • Dave

    Dave - 2012-01-10

    Hello:

    I am having a problem fitting a Weibull model with right-censored data. I use
    dinterval in the BUGS code, but when I run the (simple) program on simulated
    data (where I know the true coefficients), the results are not even close.
    When I code the log likelihood myself and use the 'Poisson trick' described in
    Ntzoufras' text (page 276) everything works fine.

    I am not sure of the reliability of the dinterval distribution. Perhaps I am
    not using it correctly, but the documentation for its use is sparse.

    Here is my code that works, followed by my original code with dinterval which
    does not work:

    model

    {C<-100

    for(i in 1:N){

    zeros_~dpois(zeros.mean_)

    zeros.mean_<- -l_+C

    l_<-y.cens_(log(beta/lambda_)+(beta-1)log(y_/lambda_)-(y_/lambda_)^beta)+(1-
    y.cens_)*(-(y_/lambda_)^beta)

    log(lambda_)<-alpha0+alpha1/temp_

    }

    beta~dgamma(0.1,0.1)

    alpha0~dnorm(0,0.01)

    alpha1~dnorm(0,0.000001)

    }

    Original code:

    model

    {for(i in 1:N){

    is.censored~dinterval(y_,y.cens_)

    y_~dweib(beta,lambda_)

    log(eta_)<- alpha0+alpha1/temp_

    lambda_<-pow(eta_,-beta)

    }

    beta~dgamma(0.1,0.1)

    alpha0~dnorm(0.0,0.01)

    alpha1~dnorm(0.0,0.000001)

    }

    Any suggestions would be appreciated.______

     
  • Martyn Plummer

    Martyn Plummer - 2012-01-10

    You have two problems. Firstly, there are two distinct parameterizations of
    the Weibull distribution and you are using the wrong one. BUGS uses the
    proportional hazards parameterization, whereas you are using the shape/scale
    parameterization with the zeros trick.

    Secondly you need to change the data representation. In a standard
    representation of survival analysis, we start with theoretical failure time T
    and censoring time C. But we do not work with T and C. Instead we observe

    Y = min(T, C)
    

    and status indicator

    Y.cens  = I(T > C)
    

    In JAGS, you work with the original T and C. C is known for all observations
    (at least theoretically, C is the end of the observation period). T is missing
    if the observation is censored. Then the model is described as

    T[j] ~dweib(lambda[j], mu)
    Y.cens[j] ~ dinterval(T[j], C[j])
    

    The reason why JAGS uses this counter-intuitive parameterization is because a
    JAGS model should describe the way that the data are generated. Concretely,
    you should be able to generate a replicate data set from a JAGS model.

     
  • Dave

    Dave - 2012-01-10

    Thank you for the quick response. I do understand the nature of censored time-
    to-event data, but I don't see how the specification of dinterval generates
    the correct likelihood. Here is the code I used, again with incorrect results:

    model

    {for(i in 1:N){

    y_~dweib(beta,lambda_)

    is.censored_~dinterval(y_,y.cens_)

    log(eta_)<-alpha0+alpha1/temp_

    lambda_<-pow(eta_,-beta)

    }

    beta~dgamma(0.1,0.1)

    alpha0~dnorm(0,0.01)

    alpha1~dnorm(0,0.000001)

    }

    The true values of beta, alpha0 and alpha1 are 2, 2.5 and 800, respectively.
    This code produces estimates far from these; for example, the estimated alpha1
    is 76. I am trying to replicate example 7.10 on page 239 of Hamada, et al's
    text on Bayesian Reliability. Can you see what I am doing wrong?___

     
  • Dave

    Dave - 2012-01-11

    To elaborate on my previous message. When I run the following WinBUGS program,
    I get the correct results:

    model

    {for(i in 1:40){

    y_~dweib(beta,lambda_)I(y.cens_,)

    log(eta_)<- alpha0+alpha1/temp_

    lambda_<-pow(eta_,-beta)

    }

    beta~dgamma(0.1,0.1)

    alpha0~dnorm(0.0,0.01)

    alpha1~dnorm(0.0,0.000001)

    }

    Data

    list( y=c(80.7 , NA, NA , NA , NA ,29.1, NA , NA , NA, NA, 47.5, 73.7, NA , NA
    ,86.2 , NA , NA , NA , NA ,71.8, 29.5 , NA ,52.0, 63.5 , NA, 99.5, 56.3, 92.5,
    NA, NA ,80.9 ,76.6, 53.4,NA ,47.5, 26.1 ,77.6 , NA ,61.8 ,56.1),y.cens=c(0,
    100 ,100 ,100 ,100 , 0 ,100 ,100, 100, 100, 0, 0, 100, 100, 0, 100, 100, 100,
    100, 0, 0, 100, 0 , 0 ,100 , 0 , 0 , 0 ,100 ,100 , 0 , 0 , 0 ,100 , 0 , 0 , 0
    ,100 , 0 , 0),temp=c(300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 350,
    350, 350, 350, 350, 350, 350, 350, 350, 350, 400, 400, 400, 400, 400, 400,
    400, 400, 400, 400, 500, 500, 500, 500, 500, 500, 500, 500, 500, 500))

    Inits

    list(beta=1,alpha0=0,alpha1=0)

    Output:

    node mean sd MC error 2.5% median 97.5% start sample

    alpha0 3.097 0.5645 0.04656 2.052 3.075 4.251 4001 16000

    alpha1 646.6 227.4 18.9 198.7 647.6 1080.0 4001 16000

    beta 2.295 0.4602 0.01487 1.485 2.268 3.289 4001 16000

    However, when I run the code in my previous message with the same data, I get
    the following output:

    Mean SD Naive SE Time-series SE

    beta 3.280 0.6283 0.003628 0.009807

    alpha0 4.063 0.4552 0.002628 0.030315

    alpha1 79.322 183.1883 1.057638 12.252298

    Why the discrepancy?


     
  • Martyn Plummer

    Martyn Plummer - 2012-01-11

    You've set up the model so that your outcome variable y is always higher than
    y.cens. For the non-censored observations, y.cens=0 so this is non-
    informative. To translate this into JAGS, you need to define an indicator
    variable - the outcome of the dinterval distribution, which is always 1. We
    can define such a variable in a data block:

    data {
       for(j in 1:40){
          one[j] <- 1
       }
    }
    model {
       for(j in 1:40){
          one[j] ~ dinterval(y[j], y.cens[j])      
          y[j]~dweib(beta,lambda[j])
          log(eta[j])<- alpha0 + alpha1/temp[j]
          lambda[j]<-pow(eta[j],-beta)
       }
       beta~dgamma(0.1,0.1)
       alpha0~dnorm(0.0,0.01)
       alpha1~dnorm(0.0,0.000001)
    }
    

    Then to run the model you just need to ensure that the initial conditions y >
    y.cens are fulfilled

    weib.data <- list( y=c(80.7 , NA,   NA ,  NA ,  NA ,29.1,   NA ,  NA ,  NA,   NA, 47.5, 73.7,
    NA ,  NA ,86.2 ,  NA ,  NA ,  NA ,  NA ,71.8, 29.5 ,  NA ,52.0, 63.5 ,  NA,
    99.5, 56.3, 92.5,   NA,   NA ,80.9 ,76.6, 53.4,NA ,47.5, 26.1 ,77.6 ,  NA ,61.8
    ,56.1),y.cens=c(0, 100 ,100 ,100 ,100 ,  0 ,100 ,100, 100, 100,   0,   0, 100,
    100,   0, 100, 100, 100, 100,   0,   0, 100,   0 ,  0 ,100 ,  0 ,  0 ,  0 ,100
    ,100 ,  0 ,  0 ,  0 ,100  , 0  , 0  , 0 ,100 ,  0 ,  0),temp=c(300, 300, 300,
    300, 300, 300, 300, 300, 300, 300, 350, 350, 350, 350, 350, 350, 350, 350, 350,
    350, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 500, 500, 500, 500, 500,
    500, 500, 500, 500, 500))
    
    weib.inits <- list(beta=1,alpha0=0,alpha1=0)
    weib.inits$y <- with(weib.data, ifelse(is.na(y), y.cens+1, NA))
    
    library(rjags)
    m <- jags.model("weib.bug", data=weib.data, inits=weib.inits)
    s <- coda.samples(m, c("alpha0", "alpha1", "beta"), n.iter=10000, thin=10)
    
     
  • Martyn Plummer

    Martyn Plummer - 2012-01-11

    ... You can do it that way, but it would be better to say we observe failure
    times from time 0 to time C=100 and define an indicator variable that is 1 if
    the failure time is censored (y > C) and 0 otherwise.

    weib.data <- list( y=c(80.7 , NA,   NA ,  NA ,  NA ,29.1,   NA ,  NA ,  NA,   NA, 47.5, 73.7,
    NA ,  NA ,86.2 ,  NA ,  NA ,  NA ,  NA ,71.8, 29.5 ,  NA ,52.0, 63.5 ,  NA,
    99.5, 56.3, 92.5,   NA,   NA ,80.9 ,76.6, 53.4,NA ,47.5, 26.1 ,77.6 ,  NA ,61.8
    ,56.1), [b]C=100[/b],temp=c(300, 300, 300,
    300, 300, 300, 300, 300, 300, 300, 350, 350, 350, 350, 350, 350, 350, 350, 350,
    350, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 500, 500, 500, 500, 500,
    500, 500, 500, 500, 500))
    [b]weib.data$is.censored <- is.na(weib.data$y)[/b]
    

    Then the model is

    model {
       for(j in 1:40){
          [b]is.censored[j] ~ dinterval(y[j], C)[/b]      
          y[j]~dweib(beta,lambda[j])
          log(eta[j])<- alpha0 + alpha1/temp[j]
          lambda[j]<-pow(eta[j],-beta)
       }
       beta~dgamma(0.1,0.1)
       alpha0~dnorm(0.0,0.01)
       alpha1~dnorm(0.0,0.000001)
    }
    

    And to run it

    weib.inits <- list(beta=1,alpha0=0,alpha1=0)
    [b]weib.inits$y <- with(weib.data, ifelse(is.censored, C+1, NA))[/b]
    
    library(rjags)
    m <- jags.model("weib2.bug", data=weib.data, inits=weib.inits)
    s <- coda.samples(m, c("alpha0", "alpha1", "beta"), n.iter=10000, thin=10)
    
     
  • Dave

    Dave - 2012-01-12

    Many thanks. Both forms run and produce the correct results. I see the
    importance of initial values for censored data. Again, thank you. I appreciate
    your patience.

     

Log in to post a comment.