Vfrdtyky

I'm having troubles reimplementing a model from winbugs on rjags. I'm getting the Invalid parent values error which is the error you get when censoring was not correctly setup, but I can't see my mistake.

This is the original model on WinBugs:

model {

    for(i in 1 : N) {                          

        times[i] ~ dweib(v, lambda[i]) T(censor[i],)

        lambda[i] <- exp(beta0 + beta1*type[i])

        S[i]  <- exp(-lambda[i]*pow(times[i],v));

        f[i]  <- lambda[i]*v*pow(times[i],v-1)*S[i]

        h[i] <- f[i]/S[i]

    }

    beta0 ~ dnorm(0.0, 0.0001)

    beta1 ~ dnorm(0.0, 0.0001)

    v ~ dexp(0.001)



    median0 <-  pow(log(2) * exp(-beta0), 1/v)

    median1 <-  pow(log(2) * exp(-beta0-beta1), 1/v)

}

Setting up a reproducible example:

type <- as.factor(c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,

                    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0))

censor <- c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,882,892,1031,

            1033,1306,1335,0,1452,1472,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,381,0,0,0,0,0,0,0,0,0,529,0,

            0,0,0,0,0,0,0,0,945,0,0,1180,0,0,1277,1397,1512,1519)

times <-c (17,42,44,48,60,72,74,95,103,108,122,144,167,170,183,185,193,195,197,208,234,235,254,307,315,401,

           445,464,484,528,542,567,577,580,795,855,NA,NA,NA,NA,NA,NA,1366,NA,NA,1,63,105,129,182,216,250,262,

           301,301,342,354,356,358,380,NA,383,383,388,394,408,460,489,499,524,NA,535,562,675,676,748,748,778,

           786,797,NA,955,968,NA,1245,1271,NA,NA,NA,NA)

df <- tibble(type = type, censor = censor, time = times) %>%

  mutate(censor_limit = replace(censor, censor == 0, max(times, na.rm = TRUE))) %>%

  mutate(is_censored = ifelse(is.na(time), 1, 0)) %>%

  mutate(time_init = ifelse(is_censored == 1, censor_limit + 1, NA))

df$censor <- NULL

head(df)

And this is the rjags part:

m <- textConnection("model {

    for(i in 1 : N) {

        isCensored[i] ~ dinterval(times[i], censorLimit[i])

        times[i] ~ dweib(v, lambda[i])

        lambda[i] <- exp(beta0 + beta1*type[i])

        S[i]  <- exp(-lambda[i]*pow(times[i],v));

        f[i]  <- lambda[i]*v*pow(times[i],v-1)*S[i]

        h[i] <- f[i]/S[i]

    }

    beta0 ~ dnorm(0.0, 0.0001)

    beta1 ~ dnorm(0.0, 0.0001)

    v ~ dexp(0.001)



    # Median survival time 

    median0 <-  pow(log(2) * exp(-beta0), 1/v)

    median1 <-  pow(log(2) * exp(-beta0-beta1), 1/v)

}")



d <- list(N = nrow(df), times = df$time, type = df$type, isCensored = df$is_censored,

          censorLimit = df$censor_limit)

inits1 = function() {

    inits = list(v = 1, beta0 = 0, beta1=0, times = df$time_init)

}

mod <- jags.model(m, data = d, inits = inits1, n.chains = 3)

update(mod, 1e3)

mod_sim <- coda.samples(model = mod, variable.names = c("lambda", "median0", "median1"), n.iter = 5e3)

mod_csim <- as.mcmc(do.call(rbind, mod_sim))

Output:

Compiling model graph

   Resolving undeclared variables

   Allocating nodes

Graph information:

   Observed stochastic nodes: 164

   Unobserved stochastic nodes: 19

   Total graph size: 910



Initializing model

Deleting model



Error in jags.model(m, data = d, inits = inits1, n.chains = 3): Error in node h[35]

Invalid parent values