3

私の目標は、競合するリスク モデルのテストに使用できるデータ セットをシミュレートすることです。関数で簡単な例を試しているだけですsurvsim::crisk.simが、期待する結果にはなりません。

 require(survival)
 simulated_data <- survsim::crisk.sim(n = 100,
                                      foltime = 200,
                                      dist.ev = rep("weibull", 2),
                                      anc.ev = c(0.8, 0.9),
                                      beta0.ev = c(2, 4),
                                      anc.cens = 1,
                                      beta0.cens = 5,
                                      nsit = 2)

 model <- survreg(Surv(time, status) ~ 1 + strata(cause), data = simulated_data)

 exp(model$scale)

 ## cause=1  cause=2 
 ## 4.407839 2.576357

これらの数値は と同じであると予想されbeta0.evます。私が間違っているかもしれないことへのポインタ、または競合するリスクデータをシミュレートする方法に関するその他の提案。

補足: シミュレートされたデータのイベントが、リスクごとに異なるワイブル分布に従って発生することを望みます。データの階層とクラスターを指定できるようにしたいと考えています。打ち切りは、ワイブル分布またはベルヌーイ分布に従うことができます。

4

1 に答える 1

1

survreg指定された推定値を回復するには、原因固有の表記法を使用できます。

この例では、パラメーターを使用していますが、より正確な推定のためにより多くの患者を使用しています。

set.seed(101)
stack_data <- survsim::crisk.sim(n = 2000,
                                     foltime = 200,
                                     dist.ev = rep("weibull", 2),
                                     anc.ev = c(0.8, 0.9),
                                     beta0.ev = c(2, 4),
                                     anc.cens = 1,
                                     beta0.cens = 5,
                                     nsit = 2)

m1 <- survreg(Surv(time, cause==1) ~ 1, data =stack_data, dist = "weibull")
m2 <- survreg(Surv(time, cause==2) ~ 1, data = stack_data, dist = "weibull")

m1$coefficientsこれは、原因 1 の beta0.ev に近づきます

m2$coefficientsこれは、原因 2 の beta0.ev に近づきます

> m1$coefficients
(Intercept) 
   1.976449 
> m2$coefficients
(Intercept) 
   3.995716

m1$scaleこれは、原因 1 の anc.ev に近づきます

m2$scaleこれは、原因 2 の anc.ev に近づきます

> m1$scale
[1] 0.8088574
> m2$scale
[1] 0.8923334

残念ながら、これは均一な打ち切り、または不均一な打ち切りが少ない場合にのみ当てはまります (あなたの例のように)

打ち切りの危険性を高めると、切片は beta0.ev パラメータを表していません

set.seed(101)
stack_data <- survsim::crisk.sim(n = 2000,
                                     foltime = 200,
                                     dist.ev = rep("weibull", 2),
                                     anc.ev = c(0.8, 0.9),
                                     beta0.ev = c(2, 4),
                                     anc.cens = 1,
                                     beta0.cens = 2, #reduced from 5, increasing the hazard function for censoring rate
                                     nsit = 2)

m1 <- survreg(Surv(time, cause==1) ~ 1, data =stack_data, dist = "weibull")
m2 <- survreg(Surv(time, cause==2) ~ 1, data = stack_data, dist = "weibull")


> m1$coefficients
(Intercept) 
   1.531818 
> m2$coefficients
(Intercept) 
   3.553687 
> 
> m1$scale
[1] 0.8139497
> m2$scale
[1] 0.8910465
于 2020-03-04T11:27:00.087 に答える