#Preliminares

#install.packages("survival")
#install.packages("KMsurv")
#install.packages("survMisc")
#install.packages("survminer")
#install.packages("flexsurv")
#install.packages("actuar")
#install.packages("dplyr")
#install.packages("moments")
#install.packages("rms")
#install.packages("SurvRegCensCov")

library(survival)  #Kaplan-Meier con la función survfit()
library(KMsurv)
library(survMisc)
library(survminer)  #graficar curva de supervivencia estimada con ggsurvplot()
library(ggfortify)
library(flexsurv) #estimar función de supervivencia con distribuciones paramétricas con la función flexsurvreg()
library(actuar)
library(dplyr)
library(readr)
library(moments)
library(car)
library("SurvRegCensCov")

Conceptos fundamentales

Ejemplos de modelos de duración

  • La duración de una huelga;

  • el tiempo que una persona permanece en el desempleo;

  • el periodo de tiempo que hay entre que una mujer tiene pareja hasta que nace su primer hijo o entre sus sucesivas maternidades;

  • el tiempo que transcurre entre las sucesivas compras de bienes duraderos,

¿Qué es un suceso?

Un suceso es el resultado de cualquier acontecimiento que ocurre en la vida de cualquier individuo o empresa.

  • Ej. Situación laboral de un individuo; situación financiera de una empresa.

  • Sinónimos: evento o acontecimiento

¿Qué es la duración de un suceso?

La medida del tiempo que transcurre entre su inicio y su final.

  • Datos de duración: aquellos que suministran información sobre la duración de un suceso

Estado:

aquella situación descrita por una o varias fundamentales que son permanentes a lo largo del tiempo.

Tránsito:

cuando esta o estas característicasfundamentales definidoras del estado se alteran. A medida que el periodo transcurre, siempre, en cada momento, existe un riesgo de que se produzca un tránsito. Hay una dualidad entre duración y tránsito

Fdp y de función de distribución(acumulada)

  • Consideremos una población de individuos. De uno de ellos observamos el tiempo hasta el tránsito (failure) o hasta la pérdida (censura).

  • \(f(t)\) es la función de densidad de probabilidad (fdp) de la variable \(T\) en el tiempo \(t\). Puede ser continua o discreta.

  • Su función de distribución acumulada \(F(t)\) es $ P(Tt)=_0^t f(s)ds$

Función de supervivencia

\[S(t)=P(T>t)=1-F(t)\] Es la probabilidad de que la longitud de un periodo sea al menos \(t\).

Tasa de riesgo (hazard rate)

Es la probabilidad de que el periodo finalice,i.e, que haya un tránsito, en el intervalo \(t +\Delta t\), condicional a que la longitud del periodo sea al menos \(t\). En otras palabras, la tasa de riesgo es la probabilidad de que algo ocurra después de la duración \(t\), condicionada a que al menos su duración sea \(t\).

\[h(t)=\lim_{\bigtriangleup t \to 0} \frac {P( t\leq T \leq t+\bigtriangleup t|T \geq t)}{\bigtriangleup t}=\lim_{\bigtriangleup t \to 0} \frac {F(t+\bigtriangleup t )-F(t)}{\bigtriangleup tS(t)}=\dfrac{f(t)}{S(t)}\]

Función de riesgo(Hazard function)

\[\lambda(t)=\frac{f(t)}{S(t)}\] Probabilidad instantánea de dejar un estado, condicional a sobrevivir al tiempo \(t\).

  • La función de riesgo recoge los valores de la tasa de riesgo para cualquier valor de \(t\)

  • Los modelos de duración se caracterizan por el modo en que especifican la tasa de riesgo.

Función de riesgo integrado

\[\Lambda(t)=\int_0^t\lambda(s)ds=-\ln S(t)\] ### Enfoques para estimar los modelos de duración

Ajuste del modelo Básicamente, existen tres enfoques para ajustar los modelos de supervivencia:

  • El primero y quizás el más sencillo es el enfoque paramétrico,en el que asumimos una forma funcional específica para la Baseline Hazard \(\lambda_0(t)\). Los ejemplos son modelos basados en las distribuciones exponencial, Weibull, gamma y F generalizada.

-Un segundo enfoque es lo que podría llamarse una estrategia flexible o semiparamétrica , en la que hacemos suposiciones leves sobre el peligro de referencia.\(\lambda_0(t)\). Específicamente, podemos subdividir el tiempo en intervalos razonablemente pequeños y asumir que el riesgo de línea de base es constante en cada intervalo, lo que conduce a un modelo exponencial por partes.

-El tercer enfoque es una estrategia no paramétrica que se centra en la estimación de los coeficientes de regresión \(\beta\) dejando la baseline Hazard \(\lambda_0(t)\) completamente sin especificar. Este enfoque se basa en una función de verosimilitud parcial propuesta por Cox (1972) en su artículo original.

Recomendación

Se sugiere estudiar el tema en el siguiente libro.

Ejemplo: Duración de las huelgas laborales

(Example 19.8 Survival Models for Strike Duration)

Ejemplo tomado de la sección 19.4, Green pág 862: MODELS FOR DURATION.

Datos (Kennan, J. “The Duration of Contract Strikes in U.S. Manufacturing.” Journal of Econometrics, 28, 1985, pp. 5–28.)

Se tienen registros de las duraciones, en días, de 62 huelgas que comenzaron en junio de los años 1968 a 1976. Cada una involucró al menos a 1,000 trabajadores y comenzó al vencimiento o reapertura de un contrato. Kennan informó la duración real. En su encuesta, Kiefer (1985), utilizando las mismas observaciones, censuró los datos a los 80 días para demostrar los efectos de la censura. Descarga los datos aquí.

strikes<-read_csv("./TableF19-2.csv")
attach(strikes)
summary(strikes)
##        T               PROD         
##  Min.   :  1.00   Min.   :-0.10443  
##  1st Qu.: 10.25   1st Qu.:-0.03143  
##  Median : 27.00   Median : 0.01138  
##  Mean   : 42.68   Mean   : 0.01102  
##  3rd Qu.: 51.25   3rd Qu.: 0.06450  
##  Max.   :216.00   Max.   : 0.07427
#View(strikes)
head(strikes)
## # A tibble: 6 x 2
##       T   PROD
##   <dbl>  <dbl>
## 1     7 0.0114
## 2     9 0.0114
## 3    13 0.0114
## 4    14 0.0114
## 5    26 0.0114
## 6    29 0.0114

La variable PROD es el residuo de una regresión lineal del logaritmo de la producción industrial en el tiempo, el tiempo al cuadrado y un conjunto de variables ficticias mensuales. Mide la producción industrial agregada menos los componentes tendenciales y estacionales.

Para estos datos necesitamos una variable de censura para indicar si un proceso aún está en curso. Abajo creamos la variable de censura, STATUS, para todas las duraciones por debajo de 80 días.

strikes$status<-strikes$T<80
detach(strikes)
attach(strikes)
#Creación de un objeto Surv: la combinación de información entre los tiempos y su censura.
Surv(strikes$T,strikes$status) 
##  [1]   7    9   13   14   26   29   52  130+   9   37   41   49   52  119+   3 
## [16]  17   19   28   72   99+ 104+ 114+ 152+ 153+ 216+  15   61   98+   2   25 
## [31]  85+   3   10    1    2    3    3    3    4    8   11   22   23   27   32 
## [46]  33   35   43   43   44  100+   5   49    2   12   12   21   21   27   38 
## [61]  42  117+
#Surv(Time,Event), toma como argumentos: Time como el tiempo hasta la observación y Event un indicador en donde es 1 si se observó el evento y 0 si es censura.

Los tiempos censurados se caracterizan por tener el símbolo + acompañado a su derecha.

Censura. Si disponemos de una muestra de observaciones de la duración de \(N\) sucesos, la censura ocurre cuando al medir la duración de algún suceso, éste todavía permanece en el estado inicial, sin haber transitado hasta otro estado. El resultado en este caso nos dice que el fenómeno ha durado al menos la medida tomada, sin ninguna información sobre la duración real del fenómeno hasta la transición. Tránsito que ocurrirá, con seguridad, un tiempo después de nuestra medida. Así, la censura es, en general, no informativa sobre la duración de un suceso.

Estimación no-paramétrica de la función de supervivencia

Función de supervivencia (survival function): es la probabilidad de que la longitud de un periodo sea al menos \(t\).

Estimador de Kaplan-Meier (sin covariables)

# Guardando el Objeto Surv
strikes.surv<- Surv(strikes$T,strikes$status)  #Creando objeto tipo Surv
strikes.km<- survfit(strikes.surv~1, data =strikes, type = "kaplan-meier")  #Estimación Kaplan Meier

# Sin Guardar el Objeto Surv
#strikes.km <- survfit(Surv(T,status)~1, data =strikes, type = "kaplan-meier")

La estimación de la función de supervivencia se lleva a cabo con la función summary().

summary(strikes.km)  #Resumen estadístico
## Call: survfit(formula = strikes.surv ~ 1, data = strikes, type = "kaplan-meier")
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     1     62       1    0.984  0.0160        0.953        1.000
##     2     61       3    0.935  0.0312        0.876        0.999
##     3     58       5    0.855  0.0447        0.772        0.947
##     4     53       1    0.839  0.0467        0.752        0.935
##     5     52       1    0.823  0.0485        0.733        0.923
##     7     51       1    0.806  0.0502        0.714        0.911
##     8     50       1    0.790  0.0517        0.695        0.898
##     9     49       2    0.758  0.0544        0.659        0.873
##    10     47       1    0.742  0.0556        0.641        0.859
##    11     46       1    0.726  0.0567        0.623        0.846
##    12     45       2    0.694  0.0585        0.588        0.818
##    13     43       1    0.677  0.0594        0.571        0.804
##    14     42       1    0.661  0.0601        0.553        0.790
##    15     41       1    0.645  0.0608        0.536        0.776
##    17     40       1    0.629  0.0613        0.520        0.762
##    19     39       1    0.613  0.0619        0.503        0.747
##    21     38       2    0.581  0.0627        0.470        0.717
##    22     36       1    0.565  0.0630        0.454        0.702
##    23     35       1    0.548  0.0632        0.438        0.687
##    25     34       1    0.532  0.0634        0.421        0.672
##    26     33       1    0.516  0.0635        0.406        0.657
##    27     32       2    0.484  0.0635        0.374        0.626
##    28     30       1    0.468  0.0634        0.359        0.610
##    29     29       1    0.452  0.0632        0.343        0.594
##    32     28       1    0.435  0.0630        0.328        0.578
##    33     27       1    0.419  0.0627        0.313        0.562
##    35     26       1    0.403  0.0623        0.298        0.546
##    37     25       1    0.387  0.0619        0.283        0.529
##    38     24       1    0.371  0.0613        0.268        0.513
##    41     23       1    0.355  0.0608        0.254        0.496
##    42     22       1    0.339  0.0601        0.239        0.480
##    43     21       2    0.306  0.0585        0.211        0.446
##    44     19       1    0.290  0.0576        0.197        0.428
##    49     18       2    0.258  0.0556        0.169        0.394
##    52     16       2    0.226  0.0531        0.142        0.358
##    61     14       1    0.210  0.0517        0.129        0.340
##    72     13       1    0.194  0.0502        0.116        0.322

Estimador de Kaplan-Meier (con covariable)

strikes.km2<- survfit(Surv(strikes$T,strikes$status) ~PROD, data=strikes, type = "kaplan-meier")  #Estimación Kaplan Meier por grupos
#summary(strikes.km2)  #Resumen estadístico
fortify(strikes.km2)  #Resumen la información
##    time n.risk n.event n.censor       surv    std.err     upper       lower
## 1     5      2       1        0 0.50000000 0.70710678 1.0000000 0.125048827
## 2    49      1       1        0 0.00000000        Inf        NA          NA
## 3    15      3       1        0 0.66666667 0.40824829 1.0000000 0.299507130
## 4    61      2       1        0 0.33333333 0.81649658 1.0000000 0.067278391
## 5    98      1       0        1 0.33333333 0.81649658 1.0000000 0.067278391
## 6     3     11       1        0 0.90909091 0.09534626 1.0000000 0.754133845
## 7    17     10       1        0 0.81818182 0.14213381 1.0000000 0.619248987
## 8    19      9       1        0 0.72727273 0.18463724 1.0000000 0.506446765
## 9    28      8       1        0 0.63636364 0.22792115 0.9947460 0.407097577
## 10   72      7       1        0 0.54545455 0.27524094 0.9355006 0.318033634
## 11   99      6       0        1 0.54545455 0.27524094 0.9355006 0.318033634
## 12  104      5       0        1 0.54545455 0.27524094 0.9355006 0.318033634
## 13  114      4       0        1 0.54545455 0.27524094 0.9355006 0.318033634
## 14  152      3       0        1 0.54545455 0.27524094 0.9355006 0.318033634
## 15  153      2       0        1 0.54545455 0.27524094 0.9355006 0.318033634
## 16  216      1       0        1 0.54545455 0.27524094 0.9355006 0.318033634
## 17    2      9       1        0 0.88888889 0.11785113 1.0000000 0.705557502
## 18   12      8       2        0 0.66666667 0.23570226 1.0000000 0.420028359
## 19   21      6       2        0 0.44444444 0.37267800 0.9226597 0.214088528
## 20   27      4       1        0 0.33333333 0.47140452 0.8397287 0.132317867
## 21   38      3       1        0 0.22222222 0.62360956 0.7544056 0.065459105
## 22   42      2       1        0 0.11111111 0.94280904 0.7051443 0.017508018
## 23  117      1       0        1 0.11111111 0.94280904 0.7051443 0.017508018
## 24    2      3       1        0 0.66666667 0.40824829 1.0000000 0.299507130
## 25   25      2       1        0 0.33333333 0.81649658 1.0000000 0.067278391
## 26   85      1       0        1 0.33333333 0.81649658 1.0000000 0.067278391
## 27    7      8       1        0 0.87500000 0.13363062 1.0000000 0.673381937
## 28    9      7       1        0 0.75000000 0.20412415 1.0000000 0.502701841
## 29   13      6       1        0 0.62500000 0.27386128 1.0000000 0.365400279
## 30   14      5       1        0 0.50000000 0.35355339 0.9998048 0.250048822
## 31   26      4       1        0 0.37500000 0.45643546 0.9173812 0.153289602
## 32   29      3       1        0 0.25000000 0.61237244 0.8302184 0.075281393
## 33   52      2       1        0 0.12500000 0.93541435 0.7818729 0.019984067
## 34  130      1       0        1 0.12500000 0.93541435 0.7818729 0.019984067
## 35    9      6       1        0 0.83333333 0.18257419 1.0000000 0.582654795
## 36   37      5       1        0 0.66666667 0.28867513 1.0000000 0.378606461
## 37   41      4       1        0 0.50000000 0.40824829 1.0000000 0.224630348
## 38   49      3       1        0 0.33333333 0.57735027 1.0000000 0.107507139
## 39   52      2       1        0 0.16666667 0.91287093 0.9974380 0.027849128
## 40  119      1       0        1 0.16666667 0.91287093 0.9974380 0.027849128
## 41    1     18       1        0 0.94444444 0.05716620 1.0000000 0.844338248
## 42    2     17       1        0 0.88888889 0.08333333 1.0000000 0.754942745
## 43    3     16       3        0 0.72222222 0.14617634 0.9618257 0.542307142
## 44    4     13       1        0 0.66666667 0.16666667 0.9242206 0.480885649
## 45    8     12       1        0 0.61111111 0.18802536 0.8834209 0.422739379
## 46   11     11       1        0 0.55555556 0.21081851 0.8398013 0.367517864
## 47   22     10       1        0 0.50000000 0.23570226 0.7935972 0.315021269
## 48   23      9       1        0 0.44444444 0.26352314 0.7449528 0.265158882
## 49   27      8       1        0 0.38888889 0.29546842 0.6939508 0.217932680
## 50   32      7       1        0 0.33333333 0.33333333 0.6406378 0.173438256
## 51   33      6       1        0 0.27777778 0.38005848 0.5850647 0.131883705
## 52   35      5       1        0 0.22222222 0.44095855 0.5273902 0.093636007
## 53   43      4       2        0 0.11111111 0.66666667 0.4104169 0.030080829
## 54   44      2       1        0 0.05555556 0.97182532 0.3732044 0.008270053
## 55  100      1       0        1 0.05555556 0.97182532 0.3732044 0.008270053
## 56    3      2       1        0 0.50000000 0.70710678 1.0000000 0.125048827
## 57   10      1       1        0 0.00000000        Inf        NA          NA
##      strata
## 1  -0.10443
## 2  -0.10443
## 3  -0.05467
## 4  -0.05467
## 5  -0.05467
## 6  -0.03957
## 7  -0.03957
## 8  -0.03957
## 9  -0.03957
## 10 -0.03957
## 11 -0.03957
## 12 -0.03957
## 13 -0.03957
## 14 -0.03957
## 15 -0.03957
## 16 -0.03957
## 17   -0.007
## 18   -0.007
## 19   -0.007
## 20   -0.007
## 21   -0.007
## 22   -0.007
## 23   -0.007
## 24  0.00535
## 25  0.00535
## 26  0.00535
## 27  0.01138
## 28  0.01138
## 29  0.01138
## 30  0.01138
## 31  0.01138
## 32  0.01138
## 33  0.01138
## 34  0.01138
## 35  0.02299
## 36  0.02299
## 37  0.02299
## 38  0.02299
## 39  0.02299
## 40  0.02299
## 41   0.0645
## 42   0.0645
## 43   0.0645
## 44   0.0645
## 45   0.0645
## 46   0.0645
## 47   0.0645
## 48   0.0645
## 49   0.0645
## 50   0.0645
## 51   0.0645
## 52   0.0645
## 53   0.0645
## 54   0.0645
## 55   0.0645
## 56  0.07427
## 57  0.07427

Graficación de la curva de supervivencia

ggsurvplot(fit =strikes.km, data =strikes, conf.int = TRUE, title = "Curva de Supervivencia", xlab = "Tiempo(días)", ylab = "Probabilidad de supervivencia", 
    legend.title = "Estimación",legend.labs = "Kaplan-Meier")

Estimación de la función de riesgo acumulado

La estimación de la función de riesgo acumulado se puede estimar de dos formas:

  • la primera es por Máxima Verosimilitud: tomando el menos logaritmo de la función de supervivencia

\[\lambda(t)=\frac{f(t)}{S(t)}=-\frac{d}{dt}\ln S(t)\] \[\begin{aligned} \Lambda (t)= \int_0^t \lambda(u)du&=-\ln S(u)\Big|_0^t \\ &=-(\ln S(t)-\ln S(0)) =-(\ln S(t)-\ln (1)) \\ &=-\ln S(t) \end{aligned}\]

  • y el segundo es utilizando el estimador Nelson-Aalen.
ggsurvplot(strikes.km, fun = "cumhaz", xlab = "Tiempo (días)", censor = T, 
    ylab = "Riesgo Acumulado", title = "Riesgo Acumulado", legend.title = "Duración de huelgas contractuales")
## Warning in if (censor & any(df$n.censor >= 1)) {: the condition has length > 1
## and only the first element will be used

Estimación de la media, mediana y percentiles de los tiempos de supervivencia

print(strikes.km, print.rmean = TRUE)
## Call: survfit(formula = strikes.surv ~ 1, data = strikes, type = "kaplan-meier")
## 
##          n     events     *rmean *se(rmean)     median    0.95LCL    0.95UCL 
##      62.00      50.00      60.50       9.89      27.00      21.00      41.00 
##     * restricted mean with upper limit =  216

Tanto como para la estimación de la mediana, tanto como los percentiles se utiliza la función , los argumentos de la función son los siguientes:

x : Objeto tipo Survfit. probs : Vector de cuantiles que se desea estimar. conft.int* : Indicador(TRUE/FALSE) para estimar el intervalo de confianza para los cuantiles.

#Mediana
quantile(strikes.km,0.5)$quantile  #La mediana es 27 días
## 50 
## 27

Estimación paramétrica de la función de supervivencia

Empecemos estimando algunos modelos de duración sin covariables. Esto sería útil para explorar la dependencia de la duración de un proceso.

¿Qué ocurre cuando la duración del suceso se prolonga?

  • En primer lugar, que a medida que dure más un suceso más riesgo habrá de que haya un tránsito: Dependencia de la duración positiva.

  • En segundo lugar, que a medida que dure más un suceso menor riesgo habrá de que haya un tránsito: Dependencia de la duración negativa.

Primero el modelo exponencial.

library(rms)
## Loading required package: Hmisc
## Loading required package: lattice
## Loading required package: Formula
## 
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:dplyr':
## 
##     src, summarize
## The following objects are masked from 'package:base':
## 
##     format.pval, units
## Loading required package: SparseM
## 
## Attaching package: 'SparseM'
## The following object is masked from 'package:base':
## 
##     backsolve
## 
## Attaching package: 'rms'
## The following objects are masked from 'package:car':
## 
##     Predict, vif
# Estimación del modelo exponencial

exp.model<-psm(Surv(T,status)~1,strikes,dist="exponential") #Parametric Survival Model
exp.model
## Parametric Survival Model: Exponential Distribution
##  
##  psm(formula = Surv(T, status) ~ 1, data = strikes, dist = "exponential")
##  
##              Model Likelihood    Discrimination    
##                    Ratio Test           Indexes    
##  Obs   62    LR chi2     0.00    R2       0.000    
##  Events50    d.f.           0    Dxy      0.000    
##  sigma  1                        g        0.000    
##                                  gr       1.000    
##  
##       Coef   S.E.   Wald Z Pr(>|Z|)
##  [1,] 3.9688 0.1414 28.06  <0.0001 
##

Calcule el AIC en el modelo exponencial para su uso posterior en la comparación de modelos.

ll.exp<- exp.model$loglik[2]
df.exp<- nrow(strikes)-exp.model$df.residual-1
aic.exp<- -2*ll.exp+2*df.exp+1
aic.exp
## [1] 497.8781

Ploteamos la función de sobrevivencia exponencial.

S<-Survival(exp.model)
times<-seq(0,200,by=1)
plot(times,S(times,exp.model$coef),xlab="Días",ylab="S(t)",main="Función de supervivencia exponencial",type='l') 

#Alternativamente puedes usar flexsurvreg() del pkg flexsurv
flex<- flexsurvreg(Surv(T,status) ~ 1, data =strikes, dist = "exp")  #Ajuste exponencial
flex
## Call:
## flexsurvreg(formula = Surv(T, status) ~ 1, data = strikes, dist = "exp")
## 
## Estimates: 
##       est      L95%     U95%     se     
## rate  0.01890  0.01432  0.02493  0.00267
## 
## N = 62,  Events: 50,  Censored: 12
## Total time at risk: 2646
## Log-likelihood = -248.4391, df = 1
## AIC = 498.8781
plot(flex, type = "survival", ci = F, main = "Probabilidad de Supervivencia", 
    conf.int = F, col = 1, col.obs = 2, xlab = "Tiempo", ylab = "Probabilidad de Supervivencia")

Gráfica de riesgio acumulado:

#Gráfica de riesgo acumulado
plot(flex, type = "cumhaz", ci = F, main = "Riesgo Acumulado", conf.int = F, 
    col = 3, col.obs = 4, xlab = "Tiempo", ylab = "Riesgo Acumulado")

Compación entre la estimación no-paramétrica y la estimación paramétrica

flexgg<- flexsurvreg(Surv(T,status)~1, data=strikes, dist = "exp") %>% 
    summary(type = "survival") %>% data.frame

kmgg<- survfit(Surv(T,status) ~ 1, data =strikes) %>% fortify

ggplot() + geom_line(aes(time, est, col = "Paramétrico"), data = flexgg) + geom_step(aes(time, 
    surv, col = "No Paramétrico"), data = kmgg) + labs(x="Tiempo(Días)", 
    y = "Probabilidad de Supervivencia", col = "Ajustes", title = "Comparación entre curvas de supervivencia")

flexch <- flexsurvreg(Surv(T,status)~1,data =strikes, dist = "exp") %>% 
    summary() %>% data.frame

kmch <-survfit(Surv(T,status)~1,data=strikes) %>% fortify

ggplot() + geom_line(aes(time,-log(est), col = "Paramétrico"), data = flexch) +geom_step(aes(time,-log(surv), col = "No Paramétrico"), data = kmch) + 
    labs(x = "Tiempo (Semanas)", y = "Riesgo Acumulado", col = "Ajustes", title = "Comparación entre Riesgo acumulado")

Estimación del modelo Weibull

weibull.model<-psm(Surv(T,status)~1,strikes, dist="weibull")
weibull.model
## Parametric Survival Model: Weibull Distribution
##  
##  psm(formula = Surv(T, status) ~ 1, data = strikes, dist = "weibull")
##  
##                   Model Likelihood    Discrimination    
##                         Ratio Test           Indexes    
##  Obs        62    LR chi2     0.00    R2       0.000    
##  Events     50    d.f.           0    Dxy      0.000    
##  sigma 1.33954                        g        0.000    
##                                       gr       1.000    
##  
##              Coef   S.E.   Wald Z Pr(>|Z|)
##  (Intercept) 3.9092 0.1900 20.57  <0.0001 
##  Log(scale)  0.2923 0.1157  2.53  0.0115  
## 
# Calculamos el AIC para su uso posterior en la comparación de los modelos

ll.weibull<- weibull.model$loglik[2]
df.weibull<- nrow(strikes)-weibull.model$df.residual-1
aic.weibull<- -2*ll.weibull+2*df.weibull+1
aic.weibull
## [1] 492.5162
# Ploteamos la función de supervivencia y la función hazard para el modelo Weibull

S<-Survival(weibull.model)

times<-seq(0,max(T),by=1)
par(mfrow=c(2,1))
plot(times, S(times,weibull.model$coef), xlab="Días",ylab="S(t)",main="Función de supervivencia Weibull",type='l')   # plots survival curve at X*Beta hat=?

H<-Hazard(weibull.model)

plot(times, H(times,weibull.model$coef), xlab="Días",ylab="H(t)",main="Función Hazard Weibull",type='l')   # plots survival curve at X*Beta hat=?

Estimación del moldeo logNormal

lognormal.model<-psm(Surv(T,status)~1,strikes,dist="lognormal")
lognormal.model
## Parametric Survival Model: Log Normal Distribution
##  
##  psm(formula = Surv(T, status) ~ 1, data = strikes, dist = "lognormal")
##  
##                    Model Likelihood    Discrimination    
##                          Ratio Test           Indexes    
##  Obs         62    LR chi2     0.00    R2       0.000    
##  Events      50    d.f.           0    Dxy      0.000    
##  sigma 1.539578                        g        0.000    
##                                        gr       1.000    
##  
##              Coef   S.E.   Wald Z Pr(>|Z|)
##  (Intercept) 3.2619 0.2008 16.24  <0.0001 
##  Log(scale)  0.4315 0.1053  4.10  <0.0001 
## 
#Calculamos el AIC para su uso posterior en la comparación de los modelos

ll.lognormal<-lognormal.model$loglik[2]
df.lognormal<-nrow(strikes)-lognormal.model$df.residual-1
aic.lognormal<- -2*ll.lognormal+2*df.lognormal+1
aic.lognormal
## [1] 483.8879
#Ploteamos la función de supervivencia y la función hazard para el modelo logNormal

S<-Survival(lognormal.model)
times <- seq(0,max(T),by=1)

par(mfrow=c(2,1))
plot(times,S(times,lognormal.model$coef), xlab="Días",ylab="S(t)",main="Función de supervivencia lognormal",type='l')   # plots survival curve at X*Beta hat=?

H<-Hazard(lognormal.model)

plot(times,H(times,lognormal.model$coef), xlab="Días",ylab="H(t)",main="Función Hazard lognormal",type='l')   # plots survival curve at X*Beta hat=?

Estimación del modelo logistic

logistic.model <- psm(Surv(T, status) ~ 1, strikes, dist="logistic")
logistic.model
## Parametric Survival Model: Logistic Distribution
##  
##  psm(formula = Surv(T, status) ~ 1, data = strikes, dist = "logistic")
##  
##                    Model Likelihood    Discrimination    
##                          Ratio Test           Indexes    
##  Obs         62    LR chi2     0.00    R2       0.000    
##  Events      50    d.f.           0    Dxy      0.000    
##  sigma 28.09632                        g        0.000    
##                                        gr       1.000    
##  
##              Coef    S.E.   Wald Z Pr(>|Z|)
##  (Intercept) 36.9320 6.0694  6.08  <0.0001 
##  Log(scale)   3.3356 0.1250 26.68  <0.0001 
## 
#Calculamos el AIC para su uso posterior en la comparación de los modelos

ll.logistic<- logistic.model$loglik[2]
df.logistic<- nrow(strikes)-logistic.model$df.residual-1
aic.logistic<- -2*ll.logistic+2*df.logistic+1
aic.logistic
## [1] 566.4571
#Ploteamos la función de supervivencia y la función hazard para el modelo logistic

S<-Survival(logistic.model)
times<-seq(0,max(T),by=1)

par(mfrow=c(2,1))

plot(times,S(times,logistic.model$coef), xlab="Días",ylab="S(t)",main="Función de supervivencia logistic",type='l')   # plots survival curve at X*Beta hat=?

H<-Hazard(logistic.model)

plot(times,H(times,logistic.model$coef), xlab="Días",ylab="H(t)",main="Función Hazard logistico",type='l')   # plots survival curve at X*Beta hat=?

Estimación the loglogistic model

loglogistic.model<-psm(Surv(T,status)~1,strikes,dist="loglogistic")
loglogistic.model
## Parametric Survival Model: Log logistic Distribution
##  
##  psm(formula = Surv(T, status) ~ 1, data = strikes, dist = "loglogistic")
##  
##                     Model Likelihood    Discrimination    
##                           Ratio Test           Indexes    
##  Obs          62    LR chi2     0.00    R2       0.000    
##  Events       50    d.f.           0    Dxy      0.000    
##  sigma 0.8973768                        g        0.000    
##                                         gr       1.000    
##  
##              Coef    S.E.   Wald Z Pr(>|Z|)
##  (Intercept)  3.2431 0.1977 16.41  <0.0001 
##  Log(scale)  -0.1083 0.1188 -0.91  0.3619  
## 
#Calculamos el AIC para su uso posterior en la comparación de los modelos

ll.loglogistic<-loglogistic.model$loglik[2]
df.loglogistic<-nrow(strikes)-loglogistic.model$df.residual-1
aic.loglogistic<- -2*ll.loglogistic+2*df.loglogistic+1
aic.loglogistic
## [1] 484.2037
# Ploteamos la función de supervivencia y la función hazard para el modelo logistic

S<-Survival(loglogistic.model)
times<-seq(0,max(T),by=1)

par(mfrow=c(2,1))

plot(times,S(times,loglogistic.model$coef), xlab="Días",ylab="S(t)",main="Función de supervivencia loglogistic",type='l')   # plots survival curve at X*Beta hat=?

H<-Hazard(loglogistic.model)

plot(times, H(times,loglogistic.model$coef), xlab="Días",ylab="H(t)",main="Función Hazard loglogistic",type='l')   # plots survival curve at X*Beta hat=?

Comparación de modelos mediante AIC

Podemos usar el criterio AIC como herramienta para comparar los modelos. El modelo con el AIC más bajo es mejor.

Con base en esta comparación resulta que el modelo lognormal es el mejor modelo.

AIC.compare<- rbind(aic.exp,aic.weibull,aic.lognormal,aic.logistic,aic.loglogistic)

AIC.compare
##                     [,1]
## aic.exp         497.8781
## aic.weibull     492.5162
## aic.lognormal   483.8879
## aic.logistic    566.4571
## aic.loglogistic 484.2037

Comparación gráfica de los modelos

datosflex <-strikes
Dist<-c("exp", "weibull", "lognormal","llogis", "gompertz", "genf")
data.Surv<-Surv(datosflex$T,datosflex$status)

# Inicia SCRIPT#

model<-sapply(Dist, function(x) flexsurvreg(data.Surv ~1,dist = x), USE.NAMES = T,simplify = F)

model
## $exp
## Call:
## flexsurvreg(formula = data.Surv ~ 1, dist = x)
## 
## Estimates: 
##       est      L95%     U95%     se     
## rate  0.01890  0.01432  0.02493  0.00267
## 
## N = 62,  Events: 50,  Censored: 12
## Total time at risk: 2646
## Log-likelihood = -248.4391, df = 1
## AIC = 498.8781
## 
## 
## $weibull
## Call:
## flexsurvreg(formula = data.Surv ~ 1, dist = x)
## 
## Estimates: 
##        est      L95%     U95%     se     
## shape   0.7465   0.5951   0.9366   0.0864
## scale  49.8583  34.3535  72.3608   9.4753
## 
## N = 62,  Events: 50,  Censored: 12
## Total time at risk: 2646
## Log-likelihood = -244.7581, df = 2
## AIC = 493.5162
## 
## 
## $lognormal
## Call:
## flexsurvreg(formula = data.Surv ~ 1, dist = x)
## 
## Estimates: 
##          est    L95%   U95%   se   
## meanlog  3.262  2.868  3.656  0.201
## sdlog    1.540  1.252  1.893  0.162
## 
## N = 62,  Events: 50,  Censored: 12
## Total time at risk: 2646
## Log-likelihood = -240.4439, df = 2
## AIC = 484.8879
## 
## 
## $llogis
## Call:
## flexsurvreg(formula = data.Surv ~ 1, dist = x)
## 
## Estimates: 
##        est     L95%    U95%    se    
## shape   1.114   0.883   1.406   0.132
## scale  25.612  17.386  37.731   5.063
## 
## N = 62,  Events: 50,  Censored: 12
## Total time at risk: 2646
## Log-likelihood = -240.6019, df = 2
## AIC = 485.2037
## 
## 
## $gompertz
## Call:
## flexsurvreg(formula = data.Surv ~ 1, dist = x)
## 
## Estimates: 
##        est       L95%      U95%      se      
## shape  -0.02119  -0.03271  -0.00968   0.00588
## rate    0.03746   0.02549   0.05504   0.00736
## 
## N = 62,  Events: 50,  Censored: 12
## Total time at risk: 2646
## Log-likelihood = -238.7954, df = 2
## AIC = 481.5909
## 
## 
## $genf
## Call:
## flexsurvreg(formula = data.Surv ~ 1, dist = x)
## 
## Estimates: 
##        est        L95%       U95%       se       
## mu      2.89e+00   2.11e+00   3.68e+00   4.00e-01
## sigma   1.56e+00   1.24e+00   1.96e+00   1.81e-01
## Q      -5.34e-01  -1.49e+00   4.17e-01   4.85e-01
## P       5.20e-03   2.77e-46   9.74e+40   2.64e-01
## 
## N = 62,  Events: 50,  Censored: 12
## Total time at risk: 2646
## Log-likelihood = -239.8512, df = 4
## AIC = 487.7024

Curvas de supervivencia

plot(model[[1]], ci = F, conf.int = F, lty = 2, main = "Ajuste Paramétrico", 
    xlab = "Tiempo", ylab = "Probabilidad de Supervivencia")
for (i in 2:length(Dist)) plot(model[[i]], ci = F, conf.int = F, add = T, col = i + 
    1, lty = i)
## Warning in if (!add) {: the condition has length > 1 and only the first element
## will be used

## Warning in if (!add) {: the condition has length > 1 and only the first element
## will be used

## Warning in if (!add) {: the condition has length > 1 and only the first element
## will be used

## Warning in if (!add) {: the condition has length > 1 and only the first element
## will be used

## Warning in if (!add) {: the condition has length > 1 and only the first element
## will be used
legend("topright", c("KM", Dist), lty = 1:(length(Dist) + 1), col = 1:(length(Dist) + 
    1))

Curvas de riesgo acumulado

plot(model[[1]], ci = F, conf.int = F, lty = 2, main = "Ajuste Paramétrico", 
    xlab = "Tiempo", ylab = "Riesgo Acumulado", type = "cumhaz")
for (i in 2:length(Dist)) plot(model[[i]], ci = F, conf.int = F, add = T, col = i + 
    1, lty = i, type = "cumhaz")
## Warning in if (!add) {: the condition has length > 1 and only the first element
## will be used

## Warning in if (!add) {: the condition has length > 1 and only the first element
## will be used

## Warning in if (!add) {: the condition has length > 1 and only the first element
## will be used

## Warning in if (!add) {: the condition has length > 1 and only the first element
## will be used

## Warning in if (!add) {: the condition has length > 1 and only the first element
## will be used
legend("topright", c("KM", Dist), lty = 1:(length(Dist) + 1), col = 1:(length(Dist) + 
    1))

Resumen tabla de criterios de información

datosflex<-strikes
data.Surv <- Surv(datosflex$T, datosflex$status)
Dist <- c("exp", "weibull","lognormal", "llogis", "gompertz", "genf")


# Inicia Script
model<-sapply(Dist, function(x) flexsurvreg(data.Surv ~ 1+PROD, data = datosflex,dist = x), USE.NAMES = T, simplify = F)

AIC.model<-sapply(model, function(x) c(AIC = AIC(x), BIC = BIC(x), LogLik = logLik(x)), 
    simplify = T)
#AIC.model

# Ordenados por menor AIC
AIC.model[, order(AIC.model["AIC", ])]
##         gompertz lognormal    llogis      genf   weibull       exp
## AIC     476.2356  479.3866  480.3795  483.2334  484.5031  485.5954
## BIC     482.6170  485.7680  486.7609  493.8691  490.8845  489.8496
## LogLik -235.1178 -236.6933 -237.1897 -236.6167 -239.2516 -240.7977

Podemos incluir covariables en cada unos de estos modelos paramétricos

En ocasiones nos interesa modelar el efecto de la presencia de un vector de covariables o variables explicativas sobre el tiempo de supervivencia.

Reestimamos el modelo Weibull incorporando la covariable PROD

weibull.model<-psm(Surv(T,status)~1+PROD,strikes,dist="weibull")
weibull.model
## Parametric Survival Model: Weibull Distribution
##  
##  psm(formula = Surv(T, status) ~ 1 + PROD, data = strikes, dist = "weibull")
##  
##                     Model Likelihood    Discrimination    
##                           Ratio Test           Indexes    
##  Obs         62    LR chi2     11.01    R2       0.163    
##  Events      50    d.f.            1    Dxy      0.224    
##  sigma 1.215848    Pr(> chi2) 0.0009    g        0.712    
##                                         gr       2.038    
##  
##              Coef     S.E.   Wald Z Pr(>|Z|)
##  (Intercept)   4.0221 0.1862 21.60  <0.0001 
##  PROD        -13.6972 4.0794 -3.36  0.0008  
##  Log(scale)    0.1954 0.1164  1.68  0.0930  
## 
# Ploteamos la función de supervivencia y la función Hazard del modelo Weibull en la media xbeta 

xbar<-c(1,mean(PROD))

xbarbeta<-xbar%*%weibull.model$coefficients

S<-Survival(weibull.model)

#times<-seq(0,max(T),by=1)
times<-seq(0,200,by=1)
par(mfrow=c(2,1))

plot(times, S(times,xbarbeta), xlab="Días",ylab="S(t)",main="Función de supervivencia Weibull",type='l')   # plots survival curve at X*Beta hat
## Warning in t.trans - lp: Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.
H<-Hazard(weibull.model)

plot(times,H(times,xbarbeta),xlab="Días",ylab="H(t)",main="Función Hazard Weibull",type='l')   # plots survival curve at X*Beta hat
## Warning in t.trans - lp: Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

¿Qué nos dicen los coeficientes

#z = survreg(Surv(tiempo, censuras) ~ 1,dist="weibull") #pkg survival
z<-survreg(Surv(T,status)~PROD,strikes,dist="weibull")
summary(z)
## 
## Call:
## survreg(formula = Surv(T, status) ~ PROD, data = strikes, dist = "weibull")
##               Value Std. Error     z       p
## (Intercept)   4.022      0.186 21.60 < 2e-16
## PROD        -13.697      4.079 -3.36 0.00079
## Log(scale)    0.195      0.116  1.68 0.09302
## 
## Scale= 1.22 
## 
## Weibull distribution
## Loglik(model)= -239.3   Loglik(intercept only)= -244.8
##  Chisq= 11.01 on 1 degrees of freedom, p= 9e-04 
## Number of Newton-Raphson Iterations: 5 
## n= 62
#WeibullReg(Surv(T,status)~1+PROD,strikes)

# Parámetros a y b de la distribución de Weibull
alfa<-z$icoef[[1]]
sigma<- z$scale
c(1/sigma, exp(alfa))
## [1]  0.8224712 49.8582984
# Tiempo medio y coeficiente de información de Akaike
#media <- exp(alfa)*gamma(1 + sigma)
#aic <- -2*z$loglik[1]+2*2
#c(media, aic)

Interpretación 1/Scale = Weibull shape parameter. (1/1.22=0.8196)

Interpretación AFT(Accelerated failure time):

Los coeficientes son logaritmos de razones de tiempos de supervivencia, por lo que un coeficiente positivo significa una supervivencia más larga.

interpretación PH(proportional hazards model)

Para traducir el coeficiente en un modelo AFT \((\alpha_j)\) al de un modelo PH \((\beta_j)\),

\((\beta_j =-\alpha_jp)\)

donde \((p)\) es el parámetro de forma.

Para este ejemplo, el coeficiente se multiplica por -1, luego se multiplica por el parámetro de forma (1/parámetro de escala = 1/1.22 para el modelo exponencial).

\(exp(-0,516*-1*1/1.22) = 75135.56\)

Esta es la razón de riesgo que compara x y x+1.

Otro diagnóstico para la especificación del modelo utiliza la función de riesgo acumulativo.

El modelo bien especificado tendrá una función de riesgo acumulativo que comienza en cero y aumenta de forma lineal. Por ejemplo, aquí comparamos las funciones de riesgo acumulativo para los modelos weibull versus lognormal

# Reestimamos el modelo Weibull

weibull.model<-psm(Surv(T,status)~1,strikes, dist="weibull")
weibull.model
## Parametric Survival Model: Weibull Distribution
##  
##  psm(formula = Surv(T, status) ~ 1, data = strikes, dist = "weibull")
##  
##                   Model Likelihood    Discrimination    
##                         Ratio Test           Indexes    
##  Obs        62    LR chi2     0.00    R2       0.000    
##  Events     50    d.f.           0    Dxy      0.000    
##  sigma 1.33954                        g        0.000    
##                                       gr       1.000    
##  
##              Coef   S.E.   Wald Z Pr(>|Z|)
##  (Intercept) 3.9092 0.1900 20.57  <0.0001 
##  Log(scale)  0.2923 0.1157  2.53  0.0115  
## 
# Reestimamos el modelo Lognormal 

lognormal.model <-psm(Surv(T, status) ~ 1, strikes, dist="lognormal")
lognormal.model
## Parametric Survival Model: Log Normal Distribution
##  
##  psm(formula = Surv(T, status) ~ 1, data = strikes, dist = "lognormal")
##  
##                    Model Likelihood    Discrimination    
##                          Ratio Test           Indexes    
##  Obs         62    LR chi2     0.00    R2       0.000    
##  Events      50    d.f.           0    Dxy      0.000    
##  sigma 1.539578                        g        0.000    
##                                        gr       1.000    
##  
##              Coef   S.E.   Wald Z Pr(>|Z|)
##  (Intercept) 3.2619 0.2008 16.24  <0.0001 
##  Log(scale)  0.4315 0.1053  4.10  <0.0001 
## 
# Ploteamos la función de riesgo integrada(cumulative hazard function) para el modelo Weibull

S<-Survival(weibull.model)

times<-seq(0,max(T),by=1)

par(mfrow=c(2,1))

plot(times,-log(S(times,weibull.model$coef)), xlab="Días",ylab=" Función de riesgo integrado",main="Función de riesgo integrado Weibull",type='l')   
# plots survival curve at X*Beta hat

 

#Ploteamos la función de riesgo integrado(the integrated hazard function) para el modelo Lognormal

S<-Survival(lognormal.model)

times<-seq(0,max(T),by=1)

plot(times,-log(S(times,lognormal.model$coef)), xlab="Days",ylab="Función de riesgo integrado",main="Función de riesgo integrado del modelo lognormal",type='l')   # plots survival curve at X*Beta hat

# Basados en esta comparación, resulta que el modelo Weibull model is the best since it comes closest to ya que se acerca más a una formar una línea recta

TABLE 19.5 Green: Survival Distributions

Modelo de Riesgo proporcionales de Cox

coxph(Surv(T,status)~PROD,data=strikes)
## Call:
## coxph(formula = Surv(T, status) ~ PROD, data = strikes)
## 
##           coef exp(coef)  se(coef)     z       p
## PROD     9.217 10069.139     3.410 2.703 0.00686
## 
## Likelihood ratio test=7.63  on 1 df, p=0.005735
## n= 62, number of events= 50

Interpretación

Un incremento en una unidad en x=PROD aumenta el riesgo de evento en un factor de 10069.139(10069.139 veces más riesgo en comparación con la línea de base).