Шаг за шагом анализа выживаемости правой цензурированные в JAGS

Question

Это своего рода наблюдения на ранее пост на SE: ​​https://stats.stackexchange.com/questions/70858/right-censored-survival-fit-with-jags

Но здесь, я хотел бы увидеть полный R сценарий (от начала до конца), на котором анализируется выживаемость данных с правой цензурой в JAGS. Все сайты, которые я нашел, требуют очень высокого уровня владения JAGS, поэтому мне трудно понять, как перейти от одной строки кода к другой. Я знаю, что это много, чтобы спросить ...

В любом случае, вот некоторые примеры данных о выживании. Группы - t1, t2, t3. НС относятся к цензурным данным (обрезание цензора = 3).

t1 <- c(1.73, NA, NA, NA, NA,0.77, NA, NA, NA, NA, NA, NA,0.5,1.06, NA, NA, NA, NA, NA,0.46, NA) 
t2 <- c(1.42, NA, NA, NA, NA, NA,0.69,1.84, NA, NA, NA,1.47,1.6,1.8, NA, NA, NA, NA, NA,0.73, NA,1.28,3,2.97) 
t3 <- c(0.88, NA, NA,1.65,1.75, NA, NA,1.01,1.46,1.95, NA, NA,2.93, NA,0.78,1.05,1.52, NA) 


#Specify model in BUGS language 


sink("model.txt") 
cat(" 
model 
{ 




} 
",fill = TRUE) 
sink() 

#Bundle data 
data<- list() 

#Parameters monitored 
parameters<-c() 

#Initial values 
inits <- list(

# MCMC settings 
ni <- 
nt <- 
nb <- 
nc <- 


fit <- jags(data, inits, parameters, "model.txt", n.iter=ni, n.thin=nt, n.burnin=nb, n.chains=nc, working.directory = getwd()) 

Я знаю, что это много, чтобы спросить, но я потратил дни, пытаясь собрать что-то вместе, и я сбиваюсь/запутаться. Я знаю, что теперь есть пакеты для запуска такого анализа, но я действительно хочу научиться строить это с нуля и самостоятельно! Спасибо, читатели!

Answer 1

Я не занимаюсь анализом выживаемости (и вы не указываете, какой дистрибутив вы хотите использовать для этой части - есть разные варианты), но этот код должен начать работу с частичной цензурой :

library("runjags") 

    # Your data 
    t1 <- c(1.73, NA, NA, NA, NA,0.77, NA, NA, NA, NA, NA, NA,0.5,1.06, NA, NA, NA, NA, NA,0.46, NA) 
    t2 <- c(1.42, NA, NA, NA, NA, NA,0.69,1.84, NA, NA, NA,1.47,1.6,1.8, NA, NA, NA, NA, NA,0.73, NA,1.28,3,2.97) 
    t3 <- c(0.88, NA, NA,1.65,1.75, NA, NA,1.01,1.46,1.95, NA, NA,2.93, NA,0.78,1.05,1.52, NA) 

    # Combine these into a single vector to make the code cleaner 
    alldata <- rbind(cbind(t1, 1), cbind(t2, 2), cbind(t3, 3)) 
    T.obs <- alldata[,1] 
    Group <- alldata[,2] 
    N <- length(T.obs) 

    # The censoring setup - in this case 0 for T.obs < 3 and 1 for T.obs > 3 
    Censoring <- as.numeric(is.na(T.obs)) 
    Breakpoint <- 3 

    # A simple JAGS model: 
    m <- " 
    model{ 
     for(i in 1:N){ 
      # The censoring part: 
      Censoring[i] ~ dinterval(T.obs[i], Breakpoint) 
      # The regression part - you may well want to change dexp to something else: 
      T.obs[i] ~ dexp(rate[Group[i]]) 
     } 
     rate[1] ~ dgamma(0.01, 0.01) 
     rate[2] ~ dgamma(0.01, 0.01) 
     rate[3] ~ dgamma(0.01, 0.01) 

     #data# N, Censoring, Breakpoint, T.obs, Group 
     #monitor# rate, T.obs 
    } 
    " 

    # One of the things we need to do is help JAGS initialise T.obs: 
    T.obs.init <- ifelse(is.na(T.obs), 4, NA) 

    # The function call: 
    results <- run.jags(m, n.chains=2, inits=list(T.obs=T.obs.init)) 

    # Look at results: 
    results 

использует пакет runjags, который делает некоторые автоматизированные диагностики и т.д. конвергенции и позволяет сокращенное использование # данных # и # монитор # внутри кода модели, а не код R - для получения дополнительной информации об этом пакете см http://runjags.sourceforge.net/quickjags.html

[редактировать: На самом деле нет необходимости контролировать T.obs, но это демонстрирует, что неверное п. п., оцениваются как> 3, а наблюдаемые значения нестационарны, как ожидалось.