ワイブル回帰
ワイブル分布を仮定して、形状パラメータ共通、尺度パラメータを説明変数で変わるモデル
\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)