sopravvivenza Plot e la funzione di pericolosità del survreg usando curva()

Question

ho il seguente modello 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 

Vorrei tracciare la pericolosità funzione e il sopravvivenza funzione in base alle stime di cui sopra.
Non voglio usare predict() o pweibull() (come presentato qui Parametric Survival o qui SO question.

Vorrei utilizzare la funzione curve(). Tutte le idee come posso fare questo? Sembra la funzione di Weibull della survreg utilizza altre definizioni di scala e forma rispetto al solito (e diversa che, per esempio rweibull)

UPDATE:. Credo che quello che realmente lo richiedono di esprimere pericolo/sopravvivenza in funzione delle stime Intercept, age (+ other potential covariates), Scale senza utilizzare alcungià pronto Funzione.

Answer 1

I suoi parametri sono:

scale=exp(Intercept+beta*x) nel tuo esempio e diciamo per età = 40

scale=283.7 

tuo parametro di forma è il reciproco della scala che gli output del modello

shape=1/1.15 

Quindi il pericolo è:

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

La funzione di rischio cumulativo è:

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

La funzione sopravvivenza è 1-la funzione di rischio cumulativo, quindi:

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

anche controllare il pacchetto eha, e la funzione hweibull e Hweibull

Answer 2

Il first link you provided ha in realtà una spiegazione chiara sulla teoria di come funziona, insieme a un esempio piacevole. (Grazie per questo, è una buona risorsa che userò nel mio lavoro.)

Per utilizzare la funzione curve, è necessario passare alcune funzioni come argomento. È vero che la famiglia di funzioni *weibull utilizza una parametrizzazione diversa per Weibull rispetto a survreg, ma può essere facilmente trasformata, come spiegato nel primo collegamento. Inoltre, dalla documentazione in survreg:

Esistono diversi modi per parametrizzare una distribuzione di Weibull. La funzione superstite la impregna in una famiglia di scala di posizione generale, che è una parametrizzazione diversa da rispetto alla funzione rweibull e spesso porta a confusione lo .

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

Ecco un'implementazione di quella semplice trasformazione:

# 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)) 

Capisco che lei ha citato nella sua risposta che non volevi usare la funzione pweibull, ma Immagino che non volessi usarlo perché utilizza una parametrizzazione diversa. In caso contrario, si potrebbe semplicemente scrivere il proprio versione di pweibull che utilizza tale parametrizzazione survreg s':

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') 

---

Come ho già detto nei commenti, non so il motivo per cui si dovrebbe fare questo (a meno che non si tratta di compiti a casa), ma si potrebbe handcode pweibull e dweibull se ti piace:

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) 

Tali definizioni vengono direttamente dalla ?dweibull. Ora basta avvolgere quelle funzioni più lente, non testate invece di pweibull e dweibull direttamente.