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.
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