JAGS에서 '계수 프로세스' 형태로 파라메트릭 생존 모델 표현
저는 JAGS에서 시간에 따라 변하는 공변량을 허용하는 생존 모델을 만들려고 합니다.예를 들어 생존이 Weibull 분포를 따른다고 가정하면 모수 모형이 되었으면 합니다(그러나 위험이 변동될 수 있으므로 지수는 너무 단순합니다).그래서, 이것은 본질적으로 베이지안 버전에서 할 수 있는 것입니다.flexsurv패키지 - 모수 모형에서 시간 간격 공변량을 허용합니다.
따라서 각 주체가 공변량이 일정하게 유지된 시간 간격(이 PDF 또는 여기에 설명된 대로)에 해당하는 여러 행을 갖는 '계산 과정' 형식으로 데이터를 입력할 수 있기를 원합니다.여기가 바로(start, stop]공식적으로 그것은.survival또는flexurv패키지가 허용합니다.
안타깝게도 JAGS에서 생존 분석을 수행하는 방법에 대한 모든 설명은 피험자당 하나의 행을 가정하는 것처럼 보입니다.
이 간단한 방법을 사용하여 계수 프로세스 형식으로 확장하려고 했지만 모형이 분포를 정확하게 추정하지 못합니다.
실패한 시도:
여기 예가 있어요.먼저 몇 가지 데이터를 생성합니다.
library('dplyr')
library('survival')
## Make the Data: -----
set.seed(3)
n_sub <- 1000
current_date <- 365*2
true_shape <- 2
true_scale <- 365
dat <- data_frame(person = 1:n_sub,
true_duration = rweibull(n = n_sub, shape = true_shape, scale = true_scale),
person_start_time = runif(n_sub, min= 0, max= true_scale*2),
person_censored = (person_start_time + true_duration) > current_date,
person_duration = ifelse(person_censored, current_date - person_start_time, true_duration)
)
person person_start_time person_censored person_duration
(int) (dbl) (lgl) (dbl)
1 1 11.81416 FALSE 487.4553
2 2 114.20900 FALSE 168.7674
3 3 75.34220 FALSE 356.6298
4 4 339.98225 FALSE 385.5119
5 5 389.23357 FALSE 259.9791
6 6 253.71067 FALSE 259.0032
7 7 419.52305 TRUE 310.4770
그런 다음 데이터를 대상별로 두 개의 관측치로 나눕니다.저는 단지 각 주제를 시간 = 300으로 나누고 있습니다. (그들이 시간 = 300에 도달하지 못한 경우, 그들은 단지 하나의 관측치를 얻습니다.)
## Split into multiple observations per person: --------
cens_point <- 300 # <----- try changing to 0 for no split; if so, model correctly estimates
dat_split <- dat %>%
group_by(person) %>%
do(data_frame(
split = ifelse(.$person_duration > cens_point, cens_point, .$person_duration),
START = c(0, split[1]),
END = c(split[1], .$person_duration),
TINTERVAL = c(split[1], .$person_duration - split[1]),
CENS = c(ifelse(.$person_duration > cens_point, 1, .$person_censored), .$person_censored), # <— edited original post here due to bug; but problem still present when fixing bug
TINTERVAL_CENS = ifelse(CENS, NA, TINTERVAL),
END_CENS = ifelse(CENS, NA, END)
)) %>%
filter(TINTERVAL != 0)
person split START END TINTERVAL CENS TINTERVAL_CENS
(int) (dbl) (dbl) (dbl) (dbl) (dbl) (dbl)
1 1 300.0000 0 300.0000 300.00000 1 NA
2 1 300.0000 300 487.4553 187.45530 0 187.45530
3 2 168.7674 0 168.7674 168.76738 1 NA
4 3 300.0000 0 300.0000 300.00000 1 NA
5 3 300.0000 300 356.6298 56.62979 0 56.62979
6 4 300.0000 0 300.0000 300.00000 1 NA
이제 JAGS 모델을 설정할 수 있습니다.
## Set-Up JAGS Model -------
dat_jags <- as.list(dat_split)
dat_jags$N <- length(dat_jags$TINTERVAL)
inits <- replicate(n = 2, simplify = FALSE, expr = {
list(TINTERVAL_CENS = with(dat_jags, ifelse(CENS, TINTERVAL + 1, NA)),
END_CENS = with(dat_jags, ifelse(CENS, END + 1, NA)) )
})
model_string <-
"
model {
# set priors on reparameterized version, as suggested
# here: https://sourceforge.net/p/mcmc-jags/discussion/610036/thread/d5249e71/?limit=25#8c3b
log_a ~ dnorm(0, .001)
log(a) <- log_a
log_b ~ dnorm(0, .001)
log(b) <- log_b
nu <- a
lambda <- (1/b)^a
for (i in 1:N) {
# Estimate Subject-Durations:
CENS[i] ~ dinterval(TINTERVAL_CENS[i], TINTERVAL[i])
TINTERVAL_CENS[i] ~ dweibull( nu, lambda )
}
}
"
library('runjags')
param_monitors <- c('a', 'b', 'nu', 'lambda')
fit_jags <- run.jags(model = model_string,
burnin = 1000, sample = 1000,
monitor = param_monitors,
n.chains = 2, data = dat_jags, inits = inits)
# estimates:
fit_jags
# actual:
c(a=true_shape, b=true_scale)
분할점의 위치에 따라 모형은 기본 분포에 대해 매우 다른 모수를 추정합니다.데이터가 계수 프로세스 양식으로 분할되지 않은 경우에만 매개 변수를 올바르게 가져옵니다.이런 종류의 문제에 대한 데이터를 포맷하는 방식은 아닌 것 같습니다.
만약 제가 가정을 놓치고 제 문제가 JAGS와 덜 관련이 있고 문제를 공식화하는 방법과 더 관련이 있다면, 제안은 매우 환영할 것입니다.저는 시간 변동 공변량이 파라메트릭 생존 모델에서 사용될 수 없다는 것에 실망할 수 있습니다(그리고 일정한 위험을 가정하고 실제로 기초 분포를 추정하지 않는 콕스 모델과 같은 모델에서만 사용될 수 있습니다). 그러나 위에서 언급했듯이,flexsurvregR의 패키지는 다음을 수용합니다.(start, stop]파라메트릭 모델의 공식화.
이와 같은 모델을 다른 언어(예: JAGS 대신 STAN)로 구축하는 방법을 아시는 분도 계시다면 감사하겠습니다.
편집:
Chris Jackson은 이메일을 통해 몇 가지 유용한 조언을 제공합니다.
JAGS에서의 절단을 위한 T() 구조가 여기에 필요하다고 생각합니다.기본적으로 사람이 살아 있지만 공변량이 일정한 각 기간(t[i], t[i+1])에 대해 생존 시간은 해당 기간의 시작 부분에서 왼쪽으로 잘리고 마지막 부분에서도 오른쪽으로 관측 중단됩니다.그래서 당신은 다음과 같은 것을 쓸 것입니다.
y[i] ~ dweib(shape, scale[i])T(t[i], )
저는 이 제안을 다음과 같이 구현하려고 했습니다.
model {
# same as before
log_a ~ dnorm(0, .01)
log(a) <- log_a
log_b ~ dnorm(0, .01)
log(b) <- log_b
nu <- a
lambda <- (1/b)^a
for (i in 1:N) {
# modified to include left-truncation
CENS[i] ~ dinterval(END_CENS[i], END[i])
END_CENS[i] ~ dweibull( nu, lambda )T(START[i],)
}
}
불행하게도 이것은 그다지 효과가 없습니다.이전 코드에서는 모델이 대부분 척도 매개변수를 정확하게 파악하고 있었지만 형상 매개변수는 매우 잘못 처리했습니다.이 새 코드를 사용하면 올바른 형상 모수에 매우 가깝지만 척도 모수를 지속적으로 과대평가합니다.과도한 추정의 정도는 분할 지점이 얼마나 늦게 오는지와 상관이 있다는 것을 알게 되었습니다.분할점이 이른 경우(cens_point = 50), 실제로 과대평가는 없습니다; 만약 늦는다면 (cens_point = 350), 많이 있습니다.
저는 이 문제가 관측치를 '이중 관측치'로 바꾸는 것과 관련이 있을 수 있다고 생각했습니다. 만약 우리가 t=300에서 관측치를 관측한다면, 같은 사람으로부터 관측치가 t=400에서 관측치를 관측할 수 있습니다.이 사람이 Weibull 매개 변수에 대한 우리의 추론에 두 개의 데이터 포인트를 기여하는 것은 직관적으로 보이는데, 실제로는 그들이 한 포인트만 기여해야 합니다.따라서, 저는 각 사람에 대한 무작위 효과를 통합하려고 노력했지만, 이것은 완전히 실패했고, 50-90 범위의 엄청난 추정치를 가지고 있었습니다.nu매개 변수왜 그런지는 모르겠지만, 아마 그것은 별도의 게시물에 대한 질문일 것입니다.문제가 관련된 것인지 아닌지는 알 수 없으므로, 해당 모델의 JAGS 코드를 포함하여 전체 게시물에 대한 코드를 여기에서 찾을 수 있습니다.
사용할 수 있습니다.rstanarm패키지는 STANSTAN 입니다.표준 R 공식 표기법을 사용하여 생존 모형을 설명할 수 있습니다. stan_surv함수는 "임의 프로세스" 형식의 인수를 허용합니다.Weibull을 포함한 다양한 기저 위험 함수를 사용하여 모형을 적합시킬 수 있습니다.
의 생존 .rstanarm-stan_surv이 기능은 여전히 CRAN에서 사용할 수 없으므로 mc-stan.org 에서 직접 패키지를 설치해야 합니다.
install.packages("rstanarm", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))
아래 코드를 참조하십시오.
library(dplyr)
library(survival)
library(rstanarm)
## Make the Data: -----
set.seed(3)
n_sub <- 1000
current_date <- 365*2
true_shape <- 2
true_scale <- 365
dat <- data_frame(person = 1:n_sub,
true_duration = rweibull(n = n_sub, shape = true_shape, scale = true_scale),
person_start_time = runif(n_sub, min= 0, max= true_scale*2),
person_censored = (person_start_time + true_duration) > current_date,
person_duration = ifelse(person_censored, current_date - person_start_time, true_duration)
)
## Split into multiple observations per person: --------
cens_point <- 300 # <----- try changing to 0 for no split; if so, model correctly estimates
dat_split <- dat %>%
group_by(person) %>%
do(data_frame(
split = ifelse(.$person_duration > cens_point, cens_point, .$person_duration),
START = c(0, split[1]),
END = c(split[1], .$person_duration),
TINTERVAL = c(split[1], .$person_duration - split[1]),
CENS = c(ifelse(.$person_duration > cens_point, 1, .$person_censored), .$person_censored), # <— edited original post here due to bug; but problem still present when fixing bug
TINTERVAL_CENS = ifelse(CENS, NA, TINTERVAL),
END_CENS = ifelse(CENS, NA, END)
)) %>%
filter(TINTERVAL != 0)
dat_split$CENS <- as.integer(!(dat_split$CENS))
# Fit STAN survival model
mod_tvc <- stan_surv(
formula = Surv(START, END, CENS) ~ 1,
data = dat_split,
iter = 1000,
chains = 2,
basehaz = "weibull-aft")
# Print fit coefficients
mod_tvc$coefficients[2]
unname(exp(mod_tvc$coefficients[1]))
값과 실제 값일출력는하치과출(력,()true_shape <- 2; true_scale <- 365):
> mod_tvc$coefficients[2]
weibull-shape
1.943157
> unname(exp(mod_tvc$coefficients[1]))
[1] 360.6058
는 "STAN"을 사용하여 .rstan::get_stanmodel(mod_tvc$stanfit)STAN에서 시도한 것과 합니다.
언급URL : https://stackoverflow.com/questions/36733638/representing-parametric-survival-model-in-counting-process-form-in-jags
'programing' 카테고리의 다른 글
| R에서 언어 설정을 변경하는 방법 (0) | 2023.06.08 |
|---|---|
| Laravel 마이그레이션을 통해 mysql/mariaDB 함수를 사용할 수 없습니다. (0) | 2023.06.08 |
| 종속성 속성.등록() 또는 .첨부() 등록 (0) | 2023.06.03 |
| xact_abort가 켜져 있을 때 raise 오류가 발생한 후 SQL Server가 계속 실행되는 이유는 무엇입니까? (0) | 2023.06.03 |
| CSS 선택기 "~"는 무엇을 의미합니까? (0) | 2023.06.03 |