使用curve()绘制survreg的生存和危险函数

Question

我拥有以下的survreg模型:

Call:
survreg(formula = Surv(time = (ev.time), event = ev) ~ age, 
    data = my.data, dist = "weib")
             Value Std. Error    z        p
(Intercept) 4.0961     0.5566 7.36 1.86e-13
age         0.0388     0.0133 2.91 3.60e-03
Log(scale)  0.1421     0.1208 1.18 2.39e-01
Scale= 1.15 

Weibull distribution

我希望基于上述估计结果绘制危险函数和生存函数，但不想使用predict()或pweibull()（如此处Parametric Survival或SO问题）。

我想使用curve()函数。你有什么想法吗？似乎survreg的Weibull函数使用了与通常定义（例如rweibull）不同的比例和形状。

更新：我想要的实际上是根据估计值Intercept、age（+其他潜在协变量）和Scale来表达危险/生存函数，而不使用任何现成的*weilbull函数。

- user1252482

Answer 1

你提供的第一个链接实际上有关于这个理论的清晰解释和美妙的例子。（谢谢你提供这个好资源，我会在自己的工作中使用它。）

为了使用curve函数，你需要将一些函数作为参数传递。虽然*weibull函数族使用不同的参数化来描述Weibull函数，但正如你在第一个链接中所解释的那样，可以很容易地转换。此外，在survreg的文档中也有相关说明：

There are multiple ways to parameterize a Weibull distribution. The survreg function imbeds it in a general location-scale familiy, which is a different parameterization than the rweibull function, and often leads to confusion.

  survreg's scale  =    1/(rweibull shape)
  survreg's intercept = log(rweibull scale)

这是一个简单转换的实现例子：

# The parameters
intercept<-4.0961
scale<-1.15

par(mfrow=c(1,2),mar=c(5.1,5.1,4.1,2.1)) # Make room for the hat.
# S(t), the survival function
curve(pweibull(x, scale=exp(intercept), shape=1/scale, lower.tail=FALSE), 
      from=0, to=100, col='red', lwd=2, ylab=expression(hat(S)(t)), xlab='t',bty='n',ylim=c(0,1))
# h(t), the hazard function
curve(dweibull(x, scale=exp(intercept), shape=1/scale)
      /pweibull(x, scale=exp(intercept), shape=1/scale, lower.tail=FALSE), 
      from=0, to=100, col='blue', lwd=2, ylab=expression(hat(h)(t)), xlab='t',bty='n')
par(mfrow=c(1,1),mar=c(5.1,4.1,4.1,2.1))

我理解您在答案中提到不想使用pweibull函数，但我猜测您的原因是它使用了不同的参数化方法。否则，您可以编写自己的版本的pweibull函数，该函数使用survreg的参数化方法：

my.weibull.surv<-function(x,intercept,scale) pweibull(x,scale=exp(intercept),shape=1/scale,lower.tail=FALSE)
my.weibull.haz<-function(x,intercept,scale) dweibull(x, scale=exp(intercept), shape=1/scale) / pweibull(x,scale=exp(intercept),shape=1/scale,lower.tail=FALSE)

curve(my.weibull.surv(x,intercept,scale),1,100,lwd=2,col='red',ylim=c(0,1),bty='n')
curve(my.weibull.haz(x,intercept,scale),1,100,lwd=2,col='blue',bty='n')

---

正如我在评论中提到的，我不知道你为什么要这样做（除非这是作业），但如果你愿意的话，你可以手动编写pweibull和dweibull：

my.dweibull <- function(x,shape,scale) (shape/scale) * (x/scale)^(shape-1) * exp(- (x/scale)^shape)
my.pweibull <- function(x,shape,scale) exp(- (x/scale)^shape)

这些定义直接来自于?dweibull。现在只需要将那些速度较慢、未经测试的函数包装起来，而不是直接使用pweibull和dweibull。

- nograpes

---

感谢您详细的电子邮件，但我不想使用任何*weibull函数。是否可能将风险表达为Intercept、age（+其他协变量）和scale的函数？ - user1252482

嗯，也许你错过了我的最后一段话，在那里我向你展示如何编写一个仅仅是pweibull包装器的函数。我不知道为什么你想重写pweibull，因为它是用C编码的，速度非常快，并且经过了很好的测试。除非这只是作业。无论如何，我会向你展示如何手动编写pweibull和dweibull。 - nograpes

我认为OP没有理解R输出和风险/生存函数之间的数学联系。 - ECII

Answer 2

您的参数是：

在您的示例中，scale=exp(Intercept+beta*x)，假设年龄为40岁。

scale=283.7

您的形状参数是模型输出比例的倒数。

shape=1/1.15

然后存在的危险是：

curve((shape/scale)*(x/scale)^(shape-1), from=0,to=12,ylab=expression(hat(h)(t)), col="darkblue",xlab="t", lwd=5)

累积危险函数是：

curve((x/scale)^(shape), from=0,to=12,ylab=expression(hat(F)(t)), col="darkgreen",xlab="t", lwd=5)

生存函数是累积危险函数的补函数，因此：

curve(1-((x/scale)^(shape)), from=0,to=12,ylab=expression(hat(S)(t)), col="darkred",xlab="t", lwd=5, ylim=c(0,1))

还要查看 eha 包，以及函数 hweibull 和 Hweibull

- ECII