gepuro.net
gepulog

データ分析エンジニアによる備忘録的ブログ

ワイブル回帰と比例ハザードモデルのサンプルデータ生成と推定

ワイブル回帰

ワイブル分布を仮定して、形状パラメータ共通、尺度パラメータを説明変数で変わるモデル

\begin{align*} F(t|x) = 1 - \exp\left( - \left( \frac{t}{\eta(x)} \right)^m \right) \\ \eta(x) = \exp(\beta * x) \end{align*}

ワイブル回帰のサンプルデータ

set.seed(0)
n <- 1000 # 台数
m <- 2 # 形状パラメータ
beta <- -1 # 係数
beta0 <- 10 # 切片
x <- rnorm(n)
data.df <- data.frame(tt = rweibull(n, shape = m, scale=exp(beta * x + beta0)), x = x, status=rep(1,n))

weibregで推定

> library(eha)
> weibreg(Surv(tt, status)~x, data=data.df)
Call:
weibreg(formula = Surv(tt, status) ~ x, data = data.df)

Covariate           Mean       Coef Exp(Coef)  se(Coef)    Wald p
x                  -0.927     2.017     7.513     0.059     0.000 

log(scale)                   10.016 22386.800     0.016     0.000 
log(shape)                    0.710     2.034     0.025     0.000 

Events                    1000 
Total time at risk        32057517 
Max. log. likelihood      -10619 
LR test statistic         1495 
Degrees of freedom        1 
Overall p-value           0

survregで推定

> library(survival)
> survreg(Surv(tt, status) ~ x, data=data.df, dist="weibull")
Call:
survreg(formula = Surv(tt, status) ~ x, data = data.df, dist = "weibull")

Coefficients:
(Intercept)           x 
 10.0162268  -0.9916226 

Scale= 0.491728 

Loglik(model)= -10618.8   Loglik(intercept only)= -11366.3
    Chisq= 1494.89 on 1 degrees of freedom, p= 0 
n= 1000 

比例ハザードモデル

ベースラインのハザード(微小時間における故障率)は共通で、ハザードが説明変数によって変動するモデル \begin{align*} \lambda(t|x) = \lambda_0(t) exp(\beta*x) \end{align*}

比例ハザードモデルのサンプルデータ

set.seed(0)
n <- 1000 # 台数
x <- rbinom(n, size=1, prob = 1/3)
tt <- runif(n , min=1, max=1000)
beta <- -1
# とりあえずベースラインはワイブル分布にしておいた
# ベースラインはノンパラなので、適当に
baseline <- pweibull(1:500, scale=100, shape=2)
hazard0 <- baseline * exp(beta * 0)
hazard1 <- baseline * exp(beta * 1)
data.df <- data.frame(tt = tt, status=ifelse(runif(tt, min=0, max=1) <= ifelse(x==0,hazard0[tt], hazard1[tt]), 1, 0), x = x)

coxphで推定

> library(survival)
> coxph(Surv(tt, status) ~ x, data=data.df)
Call:
coxph(formula = Surv(tt, status) ~ x, data = data.df)


    coef exp(coef) se(coef)     z       p
x -0.979     0.376     0.15 -6.51 7.7e-11

Likelihood ratio test=51.9  on 1 df, p=5.78e-13  n= 520, number of events= 333 
   (480 observations deleted due to missingness)

似てる記事

似てない記事