394

時変共変量を可能にする JAGS で生存モデルを構築しようとしています。たとえば、生存率がワイブル分布に従うと仮定します (ただし、ハザードを変化させたいので、指数関数は単純すぎます)。したがって、これは基本的にパッケージで実行できることのベイジアン バージョンでありflexsurv、パラメトリック モデルで時変共変量を可能にします。

したがって、「カウントプロセス」形式でデータを入力できるようにしたいと考えています。各被験者には複数の行があり、それぞれが共変量が一定のままである時間間隔に対応しています (この pdfまたはここで説明されています。これはまたはパッケージが許可する配合(start, stop]。survivalflexurv

残念ながら、JAGS で生存分析を実行する方法に関するすべての説明は、被験者ごとに 1 つの行を想定しているようです。

この単純なアプローチを採用して、カウント プロセス形式に拡張しようとしましたが、モデルは分布を正しく推定しません。

失敗した試み:

例を次に示します。まず、いくつかのデータを生成します。

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

次に、データを被験者ごとに 2 つの観測値に分割します。各被験者を時間 = 300 で分割しているだけです (時間 = 300 に達しなかった場合を除き、観察結果は 1 つしか得られません)。

## 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 とはあまり関係がなく、問題の定式化に関連している場合は、提案を歓迎します。時変共変量がパラメトリック生存モデルで使用できないことに絶望しているかもしれません (また、一定のハザードを仮定し、実際には基になる分布を推定しない Cox モデルのようなモデルでのみ使用できます)。ただし、上で述べたように、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 で打ち切られていない観測が見られた場合、この人物がは、ワイブル パラメータに関する推論に 2 つのデータ ポイントを提供していますが、実際には 1 つのポイントだけを提供する必要があります。したがって、私は各人にランダム効果を組み込んでみました。ただし、これは完全に失敗し、nuパラメーターの推定値が非常に大きくなりました (50 ～ 90 の範囲)。それがなぜなのかはわかりませんが、おそらくそれは別の投稿の質問です。私は問題が関連しているかどうかわからないので、そのモデルの JAGS コードを含む、この記事全体のコードをここで見つけることができます。

4

1 に答える 1