1

我在从 rjags 上的 winbugs 重新实现模型时遇到了麻烦。我收到Invalid parent values错误，这是您在未正确设置审查时遇到的错误，但我看不到我的错误。

这是 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)
}

设置可重现的示例：

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)

这是 rjags 部分：

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))

输出：

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
4

0 回答 0