Działka przetrwanie i funkcja niebezpieczeństwa survreg za pomocą krzywej()

Question

Mam następujący 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 

chciałbym wykreślić zagrożenia funkcję i przeżycia funkcję w oparciu o powyższe szacunki .
Nie chcę używać predict() lub pweibull() (przedstawione tutaj Parametric Survival lub tutaj SO question.

Chciałbym użyć funkcji curve(). Jakieś pomysły jak można to osiągnąć? Wydaje funkcji Weibulla z survreg wykorzystuje inne definicje skali i kształtu niż zwykle (i różne, że na przykład rweibull)

UPDATE. Chyba co naprawdę wymaga to wyrazić zagrożenia/przetrwanie jako funkcja szacunków Intercept, age (+ other potential covariates), Scale bez użycia gotowego *weilbull Funkcja.

Answer 1

Twoje parametry:

scale=exp(Intercept+beta*x) w swoim przykładzie i powiedzmy, że do wieku = 40

scale=283.7 

Twój parametr kształtu jest odwrotnością skali, że wyjścia modelu

shape=1/1.15 

Następnie zagrożenie jest następujące:

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

Funkcja zagrożenia kumulują się:

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

funkcji przeżycia, jest 1-funkcja zagrożenia zbiorcze, tak:

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

również sprawdzić pakiet eha i funkcję hweibull i Hweibull

Answer 2

Rzeczywiście ma jasne wyjaśnienie teorii, jak to działa, wraz z pięknym przykładem. (Dziękuję za to, jest to niezły zasób, którego użyję w mojej pracy).

Aby użyć funkcji curve, musisz przekazać jakąś funkcję jako argument. To prawda, że ​​rodzina funkcji używa innej parametryzacji dla Weibulla niż survreg, ale może być łatwo przekształcona, jak wyjaśniono w pierwszym linku. Ponadto z dokumentacji pod numerem survreg:

Istnieje wiele sposobów parametryzowania rozkładu Weibulla. Funkcja survreg umieszcza ją w rodzinie o ogólnym przybliżeniu, która jest inną parametryzacją niż funkcja rweibull i często prowadzi do pomyłki w przypadku .

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

Oto implementacja tej prostej transformacji:

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

Rozumiem, że pan wspomniał w swojej odpowiedzi, że nie chcesz korzystać z funkcji pweibull, ale Zgaduję, że nie chcesz go używać, ponieważ używa innej parametryzacji. W przeciwnym razie, można po prostu napisać własną wersję pweibull korzystające z tego survreg jest parametryzacja:

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

---

Jak wspomniano w komentarzach, nie wiem dlaczego chcesz to zrobić (chyba że jest to praca domowa), ale można handcode pweibull i dweibull jeśli chcesz:

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) 

definicje te pochodzą prosto z ?dweibull. Teraz po prostu zapisz te, wolniejsze, nieprzetestowane funkcje zamiast pweibull i dweibull bezpośrednio.