5
$\begingroup$

Update: I got the JAGS model running and this eliminates the distracting part of my question. It's really about the proper preparation of data for dinterval() and inits.

I can't find a concrete example of preparing right-censored survival data for use in JAGS when starting with a vector of times and a vector of censored status. The examples I've seen all have a constant censoring time.

My two questions are:

  1. What is the correct way to prepare dinterval(t[i], c[i]) when you are starting with a vector of event times and a matching vector of censored/failed status?
  2. What is the reason for supplying a initial value vector of event times where uncensored times are set to NA and censored times are set to a value > than the observed censoring time? I'm really struggling to understand that.

Update: I think I have the answer to (2). Apparently JAGS complains if you try to overwrite values in the init vector, but if the vector is initialized to NA it works.

Note that my understanding of these requirements comes from looking at existing online examples. If the requirements are unfounded I'd love to learn that.

I created a short test using the aml data set from the survival package which does a survival fit then tries to fit a Weibull distribution using JAGS and add those points for comparison. Originally I thought you had to pass zeros in for the unscensored times in c[i] for dinterval(t[i], c[i]), but that causes an error:

Error in node is.censored[1] Observed node inconsistent with observed parents at initialization. Try setting appropriate initial values.

Here is my sample R script. Some lines seem to be indented too much, not sure why. I converted all tabs to spaces and it looks clean in my edit window:

# Right-censoring test
require(survival)
require(runjags)
graphics.off()
rm(fit, S, results)

# SURVIVAL FIT AND PLOT
fit <- survfit(Surv(aml$time, aml$status) ~ 1, conf.type = "log-log")
plot(fit, xlab = "time", ylab = "survival", main = "AML survival")

# JAGS WEIBULL FIT WITH CENSORING
# Model uses three different versions of observed times.
# t.orig is the original set of observed times including censored times.
# t.mod has the censored times set to NA.
# t.cens is the set I have questions about. Currently the same as t.orig
# but is that the proper choice?
#
model <- "model {
    for(i in 1:N) {
        is.censored[i] ~ dinterval(t.mod[i], t.cens[i])
        t.mod[i] ~ dweib(a, b)
        S[i] <- 1 - pweib(t.orig[i], a, b)
    }
    a ~ dgamma(3.5, 2)
    b ~ dgamma(1.5, 100)
}"

censored <- !aml$status
    t.orig <- aml$time
t.mod <- aml$time
    t.mod[censored] <- NA
    t.cens <- aml$time
# Online examples seem to indicate special treatment of the c[i]
# passed to dinterval(t[i], c[i]), changing the uncensored times.
#t.cens[!censored] <- 0 # THIS CAUSES AN ERROR

data <- dump.format(list(
    t.orig = t.orig,
    t.mod = t.mod,
    N = length(t.mod),
    t.cens = t.cens,
    is.censored = as.numeric(censored)  
))

# Apparently you must provide initial values of the time vector
# with uncensored times of NA and censored times > the observed time.
# I don't understand this requirement at all.
t.init <- rep(NA, length(t.mod))
t.init[censored] <- t.cens[censored] + 1

inits <- c(
    dump.format(list(a = 1.1, b = 3.0e-4, t.mod = t.init)),
    dump.format(list(a = 2.3, b = 0.02, t.mod = t.init))
)
results = run.jags(model, monitor=c("S", "a", "b"), data=data, inits=inits)
S <- print(results, vars="S")
points(aml$time, S[,4], col="blue", pch=19)
    segments(aml$time, S[,1], aml$time, S[,3], col="blue", lwd=2)

Thanks for taking the time to look at this.

$\endgroup$
  • $\begingroup$ I'm having trouble following the model and this seems to be where the error lies... I find a hierarchical drawing showing "is distributed as" and "is equal to" with the priors can be helpful when visualizing these. Perhaps more comments in the model text also? $\endgroup$ – dardisco Sep 24 '13 at 5:15
  • $\begingroup$ I edited the code to do a fit without the censored data and prove that the model and priors work. Basically it's failure time distributed as F(t) = dweib(a, b), and survival probability S(t) = 1 - F(t). I think the error lies in how I am preparing the event data for use with dinterval. $\endgroup$ – Eldan Sep 24 '13 at 13:12
7
$\begingroup$

I was asked to re-post this answer here from my comment at http://doingbayesiandataanalysis.blogspot.com/2012/01/complete-example-of-right-censoring-in.html The specifics of this answer relate to the model in that comment, but the concepts apply to the topic here.


The core of the JAGS model for censored data is this:

isCensored[i] ~ dinterval( y[i] , censorLimitVec[i] )
y[i] ~ dnorm( mu , tau )

The key to understanding what JAGS is doing is that JAGS automatically imputes a random value for any variable that is not specified as a constant in the data. Thus, when y[i] is NA (i.e., a missing value, not a constant), then JAGS imputes a random value for it.

But what value should it generate?

The second line of the model, above, says that y[i] should be randomly generated from a normal distribution with mean mu and precision tau.

But the first line of the model, above, puts another constraint on the randomly generated value of y[i]. That line says that whatever value of y[i] is randomly generated, it must fall on the side of censorLimitVec[i] dictated by the value of isCensored[i].

To understand this part, let's unpack the dinterval() distribution. Suppose that censorLimitVec has 3 values in it, not just 1:

censorLimitVec = c(10,20,30)

Then randomly generated values from dinterval(y,c(10,20,30)) will be either 0, 1, 2, or 3 depending on whether $y<10$, $10 < y < 20$, $20<y<30$, or $30<y$. So, if $y=15$, dinterval(y,c(10,20,30)) has output of $1$ with 100% probability. The trick is this: We instead specify the output of dinterval, and impute a random value of y that could produce it. Thus, if we say

1 ~ dinterval(y,c(10,20,30))

then y is imputed as a random value between 10 and 20.

Putting the two model statements together,

1 ~ dinterval( y , censorLimit )

y ~ dnorm( mu , tau )

means that y comes from a normal density and y must fall above the censorLimit.

Hope that helps!!

$\endgroup$
  • $\begingroup$ Because I still had problems understanding the dinterval() distribution I searched a bit and found this presentation by Martyn Plummer. Slide 14 really helped me understand the concept. $\endgroup$ – elevendollar Jan 30 '15 at 14:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.