How to write a Weibull (Parametric AFT) model in JAGS using Zero or One Trick

135 views Asked by At

I am trying to write a customized weibull AFT model in JAGS. But the output from my customized weibull aft model is significantly different from the weibull aft model using ~dweib().

I tried in three different ways:

  1. Using time[i] ~ dweib(b, λ[i]).
model{
  for(i in 1 : N) {
    is.censored[i] ~ dinterval(time[i], cen[i])
    time[i] ~ dweib(b, lambda[i])
    lambda[i] <- exp(-mu[i] * b)
    mu[i] <- beta0 + beta1 * trt[i]
  }
  
  ##priors for betas
  beta0 ~ dnorm(0, 0.001)
  beta1 ~ dnorm(0, 0.001)
  
  ##prior for b
  b ~ dgamma(0.001, 0.001)
  sigma <- pow(alpha, -1)
}
  1. Using zero or one tricks to generate a customized likelihood distribution. I use the formula of hazard and survival function from the jags manual. I am using the one trick in the code below.
data{
  for(z in 1:N){
    ones[z] <- 1
  }
  C <- 1000000
}

model{
  for(i in 1 : N) {
    is.censored[i] ~ dinterval(time[i], cen[i])
    ones[i] ~ dbern(ones.mean[i])
    
    ##one trick
    ones.mean[i] <- L[i] / C
    
    ##customize log-likelihood of the weibull distribution
    L[i] <- ifelse(is.censored[i],
            S[i],
            h[i]*S[i])
    
    h[i] <- b * lambda[i] * pow(time[i], b - 1)
    S[i] <- exp(-lambda[i] * pow(time[i], b))
    lambda[i] <- exp(-(mu + beta_formula[i]) * b)
    beta_formula[i] <- beta1 * trt[i]
  }
  
  beta1 ~ dnorm(0, 0.001)
  
  ##prior for mu and sigma of the weibull distribution
  mu ~ dnorm(0, 0.001)
  b ~ dgamma(0.001, 0.001)
  sigma <- pow(b, -1)
}
  1. I use the package flexsurv to model a frequentist's parametric AFT model with weibull distribution to confirm the results from method 1 and 2. I changed the censoring indicator into a death indicator before fitting the model.
flexsurvreg(formula = Surv(survT, dead) ~ I(trt), 
                           data = df, dist = 'weibull')

After fitting all three models, I notice that model 2 has significantly different result from model 1 and model 3 (the coefficient for treatment is 11 instead of 0.2), so clearly something is wrong with my customized weibull model.

I am still trying to figure out how to upload the dataset. In the meantime, how can I check my code for model 2? I will upload the dataset as soon as I figure it out.

0

There are 0 answers