rplotsurvival-analysisweibullparametric-equations

Plot survival and hazard function of survreg using curve()


I have the following survreg model:

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

I would like to plot the hazard function and the survival function based on the above estimates.
I don't want to use predict() or pweibull() (as presented here Parametric Survival or here SO question.

I would like to use the curve() function. Any ideas how I can accomplish this? It seems the Weibull function of the survreg uses other definitions of scale and shape than the usual (and different that for example rweibull).

UPDATE: I guess what I really require it to express hazard / survival as a function of the estimates Intercept, age (+ other potential covariates), Scale without using any ready made *weilbull function.


Solution

  • Your parameters are:

    scale=exp(Intercept+beta*x) in your example and lets say for age=40

    scale=283.7
    

    your shape parameter is the reciprocal of the scale that the model outputs

    shape=1/1.15
    

    Then the hazard is:

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

    The cumulative hazard function is:

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

    The Survival function is 1-the cumulative hazard function, so:

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

    Also check out the eha package, and the function hweibull and Hweibull

    Weibull Functions