16

我有以下 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()（如此处介绍的参数生存或此处的SO 问题。

我想使用该curve()功能。有什么想法可以做到这一点吗？似乎 survreg 的 Weibull 函数使用了与通常不同的比例和形状定义（并且不同于例如 rweibull）。

更新：我想我真正需要它来表达作为估计函数的危险/生存，Intercept而不使用任何现成的函数。age (+ other potential covariates)Scale*weilbull

4

2 回答 2

11

您提供的第一个链接实际上对它的工作原理有一个清晰的解释，还有一个可爱的例子。（谢谢你，这是一个很好的资源，我将在自己的工作中使用。）

要使用该curve函数，您需要将某个函数作为参数传递。确实，*weibull函数族对 Weibull 使用不同的参数化而不是survreg，但它可以很容易地转换，如您的第一个链接所述。此外，从文档中survreg：

有多种方法可以参数化 Weibull 分布。survreg 函数将其嵌入到一般的位置尺度族中，这是与 rweibull 函数不同的参数化，并且经常导致混淆。

  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否则，您可以简单地使用 thatsurvreg的参数化编写您自己的版本：

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。

于 2013-04-26T15:07:54.837 回答
11

你的参数是：

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)

生存函数是 1-累积风险函数，因此：

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

威布尔函数

于 2013-04-26T16:51:04.197 回答