我正在尝试在 JAGS 中编写自定义的威布尔 AFT 模型。但我定制的 weibull aft 模型的输出与使用
~dweib()
的 weibull aft 模型有显着不同。
我尝试了三种不同的方式:
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)
}
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)
}
flexsurv
使用威布尔分布对常客参数 AFT 模型进行建模,以确认方法 1 和 2 的结果。在拟合模型之前,我将审查指标更改为死亡指标。flexsurvreg(formula = Surv(survT, dead) ~ I(trt),
data = df, dist = 'weibull')
拟合所有三个模型后,我注意到模型 2 与模型 1 和模型 3 的结果显着不同(处理系数是 11 而不是 0.2),所以显然我定制的威布尔模型有问题。
我仍在尝试弄清楚如何上传数据集。同时,我如何检查模型 2 的代码?一旦我弄清楚了,我就会上传数据集。
我分享的是威布尔模型的 JAGS 模型示例,其中形状参数为伽马先验分布,截距和斜率参数为正态先验分布:
model {
# Likelihood function
for (i in 1:n) {
z[i] <- x[i] > 0
L[i] <- pow(S[t[i]], exp(-lambda * x[i])) * z[i] + (1 - z[i]) * (1 - S[t[i]])
y[i] ~ dinterval(t[i], T)
log(S[t[i]]) <- -exp(beta[1] + beta[2] * x[i])
}
# Prior distributions
lambda ~ dgamma(0.001, 0.001)
beta[1] ~ dnorm(0, 0.001)
beta[2] ~ dnorm(0, 0.001)
}
JAGS 中的 dweib 函数主要可用于指定第 i 个观测值的形状参数 b 和尺度参数 λ[i] 的威布尔分布。
model {
# Likelihood function
for (i in 1:n) {
y[i] ~ dweib(b, lambda[i])
}
# Prior distributions
b ~ dgamma(0.001, 0.001)
for (i in 1:n) {
lambda[i] ~ dgamma(0.001, 0.001)
}
}
To use the zero or one trick in JAGS
model {
# Likelihood function
for (i in 1:n) {
z[i] <- x[i] > 0
y[i] ~ dinterval(t[i], T)
log(S[t[i]]) <- -exp(beta[1] + beta[2] * x[i])
L[i] <- pow(S[t[i]], exp(-lambda * x[i])) * z[i] + (1 - z[i]) * (1 - S[t[i]])
}
# Prior distributions
lambda ~ dgamma(0.001, 0.001)
beta[1] ~ dnorm(0, 0.001)
beta[2] ~ dnorm(0, 0.001)
}