生存時間解析の一般化ガンマ分布のモデルの推定の学習記録

【学習動機】

一般化ガンマ分布について、前回色々調べたが、

一般化ガンマ分布の数式の学習記録 - 霞と側杖を食らう

実際に、生存時間解析の、一般化ガンマ分布のパラメトリックモデルで推定すると、上手く推定できないケースがちょこちょこある。一般化ガンマ分布でtime-to-event dataの人工データを生成して、それを一般化ガンマ分布のモデルで推定することにする。

【学習内容】

パッケージの読み込みと乱数の種の設定。

library(tidyverse)
library(flexsurv)
set.seed(2021)

シミュレーションのコードはこのように。100回シミュレーション。サンプルサイズ100でデータ生成。\((\mu,\sigma,Q)\)はテキトーに乱数生成。一般化ガンマ分布でモデル指定して、BFGS法で推定。

nsim <- 100    # シミュレーション回数の設定
# 推定値を保存する
m_par_true <- matrix(NA_real_, nrow=nsim, ncol=3)
m_par_est_BFGS <- matrix(NA_real_, nrow=nsim, ncol=3)
m_est_opt <- matrix(NA_real_, nrow=nsim, ncol=2)

for(i in 1:nsim){
  n <- 100    # サンプルサイズ
  # 真のパラメータ生成
  mu <- rnorm(1,mean = 0, sd=10)
  sigma <- rgamma(1, shape = 2, rate = 1)
  Q <- rnorm(1, mean = 0, sd=2)
  #Q <- 4 # Qが4あたりより大きいとn=100じゃ推定が厳しいように見える(後ほど実行)
  #Q <- -4 # Qが負値がある程度以上だとBFGS法では収束しないことが増えるように見える(後ほど実行)

  par_true <- c(mu, sigma, Q)
  names(par_true) <- c("mu", "sigma", "Q")
  
  # データ生成
  df <- data.frame(
    time = rgengamma(n = n, mu = mu, sigma = sigma, Q = Q),
    event = sample(c(0,1), size = n, replace = TRUE, prob = c(0.2, 0.8))  # 2割がランダムで打ち切り
  )
  m_par_true[i,] <- par_true
  
  tryCatch(
    {
      fit_BFGS <- flexsurvreg(Surv(time, event)~1, data = df,
                              dist = "gengamma", method = "BFGS")
      m_par_est_BFGS[i,] <- fit_BFGS$res[,1]
    },
    error=function(e){
      m_par_est_BFGS[i,] <- numeric(3)
    },
    warning=function(e){
      m_par_est_BFGS[i,] <- numeric(3)
    },
    finally = {},
    silent = TRUE
  )
}

df_res <- 
  cbind(m_par_true, m_par_est_BFGS) %>% 
    as.data.frame() %>% 
    set_names(c("mu_true","sigma_true","Q_true",
                "mu_est","sigma_est","Q_est"))

df_res2 <- df_res %>% 
  mutate(
    mu_diff = abs(mu_true-mu_est),
    sigma_diff = abs(sigma_true-sigma_est),
    Q_diff = abs(Q_true-Q_est)
  )

結果を可視化してどうなっているか確認。赤色は真値を表し、青色は推定値を表す。青色がない部分は推定ができていない部分。赤色と青色が近くにあればその回の推定は上手くいっている一方で、離れていれば上手くいっていない。軸タイトルはちゃんと変更していない。

p_mu <- 
ggplot(df_res2) +
  theme_bw() +
  geom_line(aes(x=1:nsim,y=mu_true),colour="red") + 
  geom_line(aes(x=1:nsim,y=mu_est),colour="blue")

p_sigma <- 
ggplot(df_res2) +
  theme_bw() +
  geom_line(aes(x=1:nsim,y=sigma_true),colour="red") + 
  geom_line(aes(x=1:nsim,y=sigma_est),colour="blue")

p_Q <-
ggplot(df_res2) +
  theme_bw() +
  geom_line(aes(x=1:nsim,y=Q_true),colour="red") + 
  geom_line(aes(x=1:nsim,y=Q_est),colour="blue")

gridExtra::grid.arrange(p_mu, p_sigma, p_Q, ncol = 1)

Qが真値と推定値が大きく離れているケース確認。Qが大きいと起こりやすい？

# Qが真値と推定値が大きく離れているケース確認
df_res2 %>% filter(Q_diff>1)

##      mu_true sigma_true    Q_true     mu_est  sigma_est     Q_est    mu_diff
## 1  -7.256602  2.3907504  3.309848  -5.269766 0.22579896 33.316234 1.98683683
## 2  11.170970  3.5850739  5.138105  13.496568 0.40728452 32.798580 2.32559811
## 3   5.498213  2.6399826  4.675271   6.547039 0.30846360 32.267718 1.04882608
## 4 -14.714114  0.6803831 -4.738894 -14.805919 0.49633897 -7.011095 0.09180474
## 5   2.708861  1.5552280  2.783587   3.619433 1.14671749  3.842243 0.91057210
## 6  -3.969919  1.2063532  2.131750  -2.580868 0.05956122 46.632574 1.38905148
## 7 -19.410441  1.6320494  1.888442 -18.352138 0.79317497  3.818769 1.05830329
##   sigma_diff    Q_diff
## 1  2.1649514 30.006386
## 2  3.1777894 27.660475
## 3  2.3315190 27.592447
## 4  0.1840442  2.272201
## 5  0.4085105  1.058656
## 6  1.1467919 44.500824
## 7  0.8388745  1.930327

Qが推定できていないケース確認。Qが負だと起きやすい？

# Qが推定できていないケース確認
df_res2 %>% filter(is.na(Q_diff))

##       mu_true sigma_true    Q_true mu_est sigma_est Q_est mu_diff sigma_diff
## 1 10.75032179  0.9671487 -3.786987     NA        NA    NA      NA         NA
## 2 -0.05690626  2.2481857 -3.910603     NA        NA    NA      NA         NA
##   Q_diff
## 1     NA
## 2     NA

\(Q=4\)で同じことをやって、結果確認。青線はあるが、大きく外れていることが見て取れる。

サンプルサイズを増やせば解決するかもしれないので、\(Q=4\)のまま、サンプルサイズを100から10000にしてみる。一見、離れているように見えるが、縦軸をよく見ると、真値の4付近に来ている。Qがある程度以上大きいときは、サンプルサイズが大きくないと、真値付近で尤度関数が膨れないということなのかもしれない。尤度関数の形状によっては、そういうことが起きうるんだろうなと考えられる。これは推察に過ぎないが、Qが大きいときというのは、ガンマ分布の形状母数が1より小さい場合に近いかもしれない(ガンマ分布について言うと、逆に形状母数が大きい時は、つまり指数分布をたくさん足し合わせていることになるので、正規分布に近づいく。)。生存時間が高確率でかなり短く終わる分布だから上手く推定できないのかもしれない。

\(Q=-4\)で同じことをやって、結果確認。青線が途切れることが多いのが見て取れる。

同様に、サンプルサイズを増やせば解決するかもしれないので、\(Q=-4\)のまま、サンプルサイズを100から10000にしてみる。推定できているときは\(Q=4\)のときと同じで、近いところに来ている。しかし、推定がそもそもできていないのは多い。尤度関数が変な形をしているのかもしれない。どうすべきか分からない。

【学習予定】

生存時間解析のパラメトリックモデルのざっくりしたまとめを書こうとしてまだやっていないからそのうち。