Usando la bases de datos del paper “Zhang Z. Survival analysis in the presence of competing risks. Ann Transl Med 2017;5(3):47. doi: 10.21037/atm.2016.08.62”

install.packages("riskRegression")
library(riskRegression)
data(Melanoma)
str(Melanoma) 
data_melanoma <- data.frame(Melanoma) 
View(data_melanoma)

Tambien la podemos descargar en nuestro computador como una base de datos de excel

install.packages("writexl")
library(writexl)
write_xlsx(Melanoma, "Downloads/melanoma.xlsx")
install.packages("readxl")

Las variables de las base de datos son las siguientes:

str(data_melanoma)
## Classes 'tbl_df', 'tbl' and 'data.frame':    205 obs. of  12 variables:
##  $ time    : num  10 30 35 99 185 204 210 232 232 279 ...
##  $ status  : num  2 2 0 2 1 1 1 1 2 1 ...
##  $ event   : chr  "death.other.causes" "death.other.causes" "censored" "death.other.causes" ...
##  $ invasion: chr  "level.1" "level.0" "level.1" "level.0" ...
##  $ ici     : chr  "2" "0" "2" "2" ...
##  $ epicel  : chr  "present" "not present" "not present" "not present" ...
##  $ ulcer   : chr  "present" "not present" "not present" "not present" ...
##  $ thick   : num  6.76 0.65 1.34 2.9 12.08 ...
##  $ sex     : chr  "Male" "Male" "Male" "Female" ...
##  $ age     : num  76 56 41 71 52 28 77 49 60 68 ...
##  $ logthick: num  1.911 -0.431 0.293 1.065 2.492 ...
##  $ status2 : num  1 1 0 1 1 1 1 1 1 1 ...

Transformamos variables

## [1] "character"
## [1] "censored"                 "death.malignant.melanoma"
## [3] "death.other.causes"
## [1] "character"
## [1] "level.0" "level.1" "level.2"
## [1] "character"
## [1] "not present" "present"
## [1] "character"
## [1] "not present" "present"
## [1] "character"
## [1] "Female" "Male"

Instalamos los paquetes necesarios

install.packages("car")
install.packages("survival")
install.packages("KMsurv")
install.packages("survMisc")
install.packages("survminer")
install.packages("flexsurv")
install.packages("actuar")
install.packages("tidyverse")
install.packages("dplyr")
install.packages("ggplot2")
install.packages("ggfortify")
install.packages("htmlwidgets")

Posteriomente se cargan las librerias de los paquetes instalados previamente

library(car)
library(survival)
library(KMsurv)
library(survMisc)
library(survminer)
library(ggfortify)
library(flexsurv)
library(actuar)
library(tidyverse)
library(ggplot2)
library(dplyr)
library(htmlwidgets)

Analisis de las variables. Tablas de frecuencias

table(data_melanoma$event)
table(data_melanoma$sex)
table(data_melanoma$invasion)
table(data_melanoma$ici)
table(data_melanoma$epicel)
table(data_melanoma$ulcer)
table(data_melanoma$sex, data_melanoma$event)

Abreviada

attach(data_melanoma)
table(event)
table(event,sex)
library(rpivotTable)
rpivotTable(data = data_melanoma, rows="sex",cols = "event", vals = "time", aggregatorName = "Average",rendererName = "Table")

Supervivencia por todas las causas

table(data_melanoma$status2)

#Creamos el objeto supervivencia

library(survival)
data_melanoma.surv <- Surv(data_melanoma$time, data_melanoma$status2)

Determinar probabilidad de sobrevivir por cualquier causa: A 1 año por ejemplo

summary(survfit(Surv(time, status2) ~ 1, data = data_melanoma), times = 365.25)
## Call: survfit(formula = Surv(time, status2) ~ 1, data = data_melanoma)
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##   365    193      11    0.946  0.0158        0.916        0.978

En este caso la probabilidad de sobrevivir a un año es del 95%

Estimación Kaplan Meier muerte por todas las causas

data_melanoma.km <- survfit(data_melanoma.surv ~ 1, data = data_melanoma, type = "kaplan-meier")
data_melanoma.km
## Call: survfit(formula = data_melanoma.surv ~ 1, data = data_melanoma, 
##     type = "kaplan-meier")
## 
##       n  events  median 0.95LCL 0.95UCL 
##     205      71      NA    3338      NA
summary(data_melanoma.km)

Grafico KM

library(survminer)
ggsurvplot(fit = survfit(Surv(time, status2) ~ 1, data = data_melanoma), xlab = "Days", ylab = "Overall survival probability", legend.title = "Estimación", legend.labs = "Kaplan-Meier")

Mediana de sobrevida NO ha sido alcanzada en el tiempo de seguimiento

survfit(Surv(time, status2) ~ 1, data = data_melanoma)

#KM por ulcer

data_melanoma.km.by.ulcer <- survfit(data_melanoma.surv~ulcer, data=data_melanoma,  conf.type="log-log")
data_melanoma.km.by.ulcer
## Call: survfit(formula = data_melanoma.surv ~ ulcer, data = data_melanoma, 
##     conf.type = "log-log")
## 
##                     n events median 0.95LCL 0.95UCL
## ulcer=not present 115     23     NA      NA      NA
## ulcer=present      90     48   2256    1506      NA

#Graficando en plot

plot(data_melanoma.km.by.ulcer, xlab="Survival Time in Years", ylab="% Surviving", yscale=100, col=c("red","blue"),main="Survival Distributions by ulcer")
legend ("topright", title="Ulcer", c("No Present", "Present"),fill=c("red", "blue"))

#Graficando con ggplot

ggsurvplot(fit = data_melanoma.km.by.ulcer, pval = TRUE, conf.int = TRUE, risk.table = TRUE,risk.table.col = "strata", linetype = "strata", surv.median.line = "hv",ggtheme = theme_bw(), palette = c("#E7B800", "#2E9FDF"))

#KM por sex

data_melanoma.km.by.sex <- survfit(data_melanoma.surv~sex, data=data_melanoma,  conf.type="log-log")
data_melanoma.km.by.sex
## Call: survfit(formula = data_melanoma.surv ~ sex, data = data_melanoma, 
##     conf.type = "log-log")
## 
##              n events median 0.95LCL 0.95UCL
## sex=Female 126     35     NA    3458      NA
## sex=Male    79     36   3154    2061      NA
ggsurvplot(fit = data_melanoma.km.by.sex, pval = TRUE, conf.int = TRUE, risk.table = TRUE,risk.table.col = "strata", linetype = "strata", surv.median.line = "hv",ggtheme = theme_bw(), palette = c("#E7B800", "#2E9FDF"))
## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.

#Podemos realizar pruebas de significancia entre grupos usando una prueba de log-rank

#La prueba de log-rank pondera igualmente las observaciones durante todo el tiempo de seguimiento y es la forma más común de comparar los tiempos de supervivencia entre grupos

#Se usa el valor p del log-rank para estimar diferencias en el tiempo de supervivencia segun variables de interes

#Para lo anterior usamos la funcion survdiff

survdiff(Surv(time, status2) ~ sex, data = data_melanoma)
survdiff(Surv(time, status2) ~ ulcer, data = data_melanoma)
survdiff(Surv(time, status2) ~ invasion, data = data_melanoma)
survdiff(Surv(time, status2) ~ ici, data = data_melanoma)
survdiff(Surv(time, status2) ~ epicel, data = data_melanoma)

Para obtener solo el valor p del log rank que se puede hacer para cada variable

sd <- survdiff(Surv(time, status2) ~ sex, data = data_melanoma)
1 - pchisq(sd$chisq, length(sd$n) - 1)
## [1] 0.004953186
fit = survdiff(Surv(time, status2) ~ sex, data = data_melanoma)

#Es posible que queramos cuantificar el tamaño de un efecto para una sola variable, o incluir más de una variable en un modelo de regresión para tener en cuenta los efectos de múltiples variables.

#El modelo de regresion de COX es usado para modelos de regresion que tienen resultados de supervivencia

#COX  de causa especifica  por  sex
Coxsex <- coxph(Surv(time, status2) ~ sex, data = data_melanoma)
summary(Coxsex)
## Call:
## coxph(formula = Surv(time, status2) ~ sex, data = data_melanoma)
## 
##   n= 205, number of events= 71 
## 
##           coef exp(coef) se(coef)     z Pr(>|z|)   
## sexMale 0.6559    1.9269   0.2376 2.761  0.00577 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##         exp(coef) exp(-coef) lower .95 upper .95
## sexMale     1.927      0.519      1.21      3.07
## 
## Concordance= 0.585  (se = 0.031 )
## Likelihood ratio test= 7.51  on 1 df,   p=0.006
## Wald test            = 7.62  on 1 df,   p=0.006
## Score (logrank) test = 7.9  on 1 df,   p=0.005

El modelo se puede hacer por todas las variables de interes

#COX  de causa especifica  por ulcer
Coxulcer <- coxph(Surv(time, status2) ~ ulcer, data = data_melanoma)
summary(Coxulcer)
#COX  de causa especifica  por invasion
Coxinvasion <- coxph(Surv(time, status2) ~ invasion, data = data_melanoma)
summary(Coxinvasion)
#COX  de causa especifica  por thick
Coxthick <- coxph(Surv(time, status2) ~ thick, data = data_melanoma)
summary(Coxthick)

Si queremos ver con formato los resultados de la regresion de COX

install.packages("gtsummary")
remotes::install_github("rstudio/gt", ref = gtsummary::gt_sha)
library(gtsummary)
coxph(Surv(time, status2) ~ sex, data = data_melanoma) %>% gtsummary::tbl_regression(exp = TRUE)
Characteristic HR 95% CI p-value
sex
Female —
Male 1.93 1.21, 3.07 0.006

Si incluimos las variables significativas en el modelo de COX

library(survival)
Coxsexulcerinvasionthickage <- coxph(Surv(time, status2) ~ sex + ulcer + invasion + thick + age, data = data_melanoma)
summary(Coxsexulcerinvasionthickage)

En el modelo multivariado la presencia de ulcera, la edad y el grosor en mm se mantienen como siginficativos

Coxulceragethick <- coxph(Surv(time, status2) ~ ulcer + age + thick, data = data_melanoma)
summary(Coxulceragethick)
## Call:
## coxph(formula = Surv(time, status2) ~ ulcer + age + thick, data = data_melanoma)
## 
##   n= 205, number of events= 71 
## 
##                 coef exp(coef) se(coef)     z Pr(>|z|)    
## ulcerpresent 0.99892   2.71535  0.26587 3.757 0.000172 ***
## age          0.02306   1.02333  0.00782 2.949 0.003185 ** 
## thick        0.09887   1.10392  0.03323 2.975 0.002926 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##              exp(coef) exp(-coef) lower .95 upper .95
## ulcerpresent     2.715     0.3683     1.613     4.572
## age              1.023     0.9772     1.008     1.039
## thick            1.104     0.9059     1.034     1.178
## 
## Concordance= 0.734  (se = 0.031 )
## Likelihood ratio test= 44.94  on 3 df,   p=1e-09
## Wald test            = 44.75  on 3 df,   p=1e-09
## Score (logrank) test = 50.03  on 3 df,   p=8e-11

#El valor de interés de un modelo de regresión de Cox es un HR

#El HR representa la proporción de riesgos entre dos grupos en cualquier punto particular en el tiempo.

#El HR se interpreta como la tasa instantánea de ocurrencia del evento de interés en aquellos que aún están en riesgo de sufrir el evento.

library(survminer)
ggsurvplot(fit = survfit(Surv(time, status2) ~ ulcer, data = data_melanoma), conf.int = TRUE, risk.table.col = "strata", ggtheme = theme_bw(), palette = c("#E7B800", "#2E9FDF"),fun = "cumhaz")

Tambien podemos obtener el plot de eventos acumulados

data_melanoma.km.by.ulcer <- survfit(data_melanoma.surv~ulcer, data=data_melanoma,  conf.type="log-log")
data_melanoma.km.by.ulcer
ggsurvplot(data_melanoma.km.by.ulcer, fun="event", censor=F, xlab="Days")

ggsurvplot(fit = survfit(Surv(time, status2) ~ ulcer, data = data_melanoma), conf.int = TRUE, risk.table.col = "strata", ggtheme = theme_bw(), palette = c("#E7B800", "#2E9FDF"),fun = "event")

Riesgos Competitivos

#Cuando los sujetos tienen múltiples eventos posibles en un entorno de tiempo hasta el evento, la depedencia no conocida entre los eventos es el problema fundamental que requiere un abordaje especial

Primero instalamos los paquetes necesarias

install.packages("cmprsk")
install.packages("survival")
install.packages("gridExtra")
install.packages("mstate")

Cargamos las librerias

library(cmprsk)
library(survival)
library(gridExtra)
library(mstate)

Incidencia Acumulada

La funcion de incidencia acumulada CIF puede estimar diferentes causas de falla y permite comparaciones entre grupos

library(cmprsk)
cifdata_melanoma <- cuminc(ftime = data_melanoma$time, fstatus= data_melanoma$status, cencode = 0)

Diferencias por sexo

cifdata_melanomasex <- cuminc(ftime = data_melanoma$time, fstatus= data_melanoma$status, group = data_melanoma$sex)

Los pacientes masculinos tienen un mayor riesgo de muerte por melanoma y otras causas que las mujeres

La diferencia parece mayor para la falla por melanoma versus la falla por otras causas

plot(cifdata_melanomasex,col=1:4,xlab="Days")

Para realizar una prueba estadística formal para la diferencia entre los grupos se usa el estadística χ2

cifdata_melanomasex$Tests
##        stat        pv df
## 1 5.8140209 0.0158989  1
## 2 0.8543656 0.3553203  1

Los pacientes masculinos tienen más probabilidades de morir de melanoma que las mujeres (P = 0.016), pero no existe una diferencia significativa en el riesgo de mortalidad por otras causas para pacientes masculinos versus femeninos (P = 0.36)

Graficando las incidencias acumuladas

library(mstate)
CIFdata_melanoma <- Cuminc(time = data_melanoma$time, status = data_melanoma$status)
head(CIFdata_melanoma)

ggplot(CIFdata_melanoma, aes(time)) + geom_step(aes(y = `CI.1`, colour = "Muerte por melanoma"))+ geom_step(aes(y = `CI.1`+ `CI.2`, colour = "Muerte por otra causa")) + geom_step(aes(y = Surv, colour = "Overall survival")) + labs(x = "Tiempo (años)", y = "Proporcion", colour = "") + theme_classic()

Graficando las incidencias acumulada por sex con ggplot2

ggcompetingrisks(cifdata_melanomasex, palette = "lancet",legend = "top", ggtheme = theme_bw())

diferencias por ulcera

cifdata_melanomaulcer <- cuminc(ftime = data_melanoma$time, fstatus= data_melanoma$status, group = data_melanoma$ulcer)
plot(cifdata_melanomaulcer,col=1:4,xlab="Days")

cifdata_melanomaulcer$Tests
##        stat           pv df
## 1 26.120719 3.207240e-07  1
## 2  0.158662 6.903913e-01  1

Para graficar las incidencias acumuladas debemos crear el siguiente factor

data_melanoma$sex <- as.factor(data_melanoma$sex)
levels(data_melanoma$sex)
## [1] "Female" "Male"
cifdata_melanomasex <- Cuminc(time = "time", status = "status", data= data_melanoma, group = data_melanoma$sex)
head(cifdata_melanomasex)
##    group time      Surv        CI.1        CI.2      seSurv      seCI.1
## 1 Female   99 0.9920635 0.000000000 0.007936508 0.007904951 0.000000000
## 2 Female  232 0.9841270 0.000000000 0.015873016 0.011134482 0.000000000
## 3 Female  279 0.9761905 0.007936508 0.015873016 0.013581801 0.007904951
## 4 Female  295 0.9682540 0.015873016 0.015873016 0.015619031 0.011134482
## 5 Female  355 0.9603175 0.015873016 0.023809524 0.017390892 0.011134482
## 6 Female  386 0.9523810 0.023809524 0.023809524 0.018971883 0.013581801
##        seCI.2
## 1 0.007904951
## 2 0.011134482
## 3 0.011134482
## 4 0.011134482
## 5 0.013581801
## 6 0.013581801
library(gridExtra)
grid.arrange(
  ggplot(cifdata_melanomasex, aes(time)) +
    geom_step(aes(y = `CI.1`, colour = group)) +
    labs(x = "Time (dias)", y = "Proporcion", title = "Muerte por Melanoma") + 
    theme_test(),
  ggplot(cifdata_melanomasex, aes(time)) +
    geom_step(aes(y = `CI.2`, colour = group)) +
    labs(x = "Time (dias)", y = "Proporcion", title = "Muerte por otra causa") + 
    theme_test(),ncol = 2)

#Graficar incidencia acumulada por evento

#Instalar manualmente ‘ CumIncidence.R ’ from http://www.stat.unipg.it/luca/R

Ahora si podemos hacer el analisis

fit<-CumIncidence (data_melanoma$time, data_melanoma$sex, data_melanoma$event, cencode = 1, xlab = "dias")

col<-c(2,3,4)
fit<-CumIncidence (data_melanoma$time, data_melanoma$sex, data_melanoma$event, cencode = 0, xlab = "dias", col= col)

fit<-CumIncidence (data_melanoma$time, data_melanoma$sex, data_melanoma$event, cencode = 0, xlab = "dias", t=c(0, 3, 6, 9, 12, 24, 36, 48), col=col)

par(mfrow=c(2,3))
fit<-CumIncidence (data_melanoma$time, data_melanoma$sex, data_melanoma$event, cencode = 0, xlab = "dias", col= col,level = 0.95)

par(mfrow=c(1,1))
fit<-CumIncidence (data_melanoma$time, data_melanoma$sex, data_melanoma$event, cencode = 0, xlab = "dias", col= col,level = 0.95)

#Modelo de regresion hazard de causa especifica

#Para ajustar el modelo de COX se introdujo el modelo de Fine Gray

#Modelo subdistribuciones o Fine gray en ek cual la incidencia acumulativa esta asociada con el hazard de causa especifica

install.packages("riskRegression")
install.packages("prodlim")
library (riskRegression)
library(prodlim)

#La funcion coxph() considera el status ==1 como evento y todo lo demas como censura por lo cual un modelo con las mismas variables difiere del obtenido por causa especifica

cshdata_melanoma <- coxph(Surv(time,status==1)~sex+age+invasion,data=data_melanoma)
summary(cshdata_melanoma)
## Call:
## coxph(formula = Surv(time, status == 1) ~ sex + age + invasion, 
##     data = data_melanoma)
## 
##   n= 205, number of events= 57 
## 
##                     coef exp(coef) se(coef)     z Pr(>|z|)    
## sexMale         0.663383  1.941349 0.266320 2.491 0.012741 *  
## age             0.009823  1.009871 0.008339 1.178 0.238840    
## invasionlevel.1 1.037168  2.821217 0.328241 3.160 0.001579 ** 
## invasionlevel.2 1.403225  4.068300 0.380744 3.685 0.000228 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                 exp(coef) exp(-coef) lower .95 upper .95
## sexMale             1.941     0.5151    1.1519     3.272
## age                 1.010     0.9902    0.9935     1.027
## invasionlevel.1     2.821     0.3545    1.4826     5.368
## invasionlevel.2     4.068     0.2458    1.9290     8.580
## 
## Concordance= 0.7  (se = 0.035 )
## Likelihood ratio test= 26.7  on 4 df,   p=2e-05
## Wald test            = 24.39  on 4 df,   p=7e-05
## Score (logrank) test = 26.85  on 4 df,   p=2e-05

Con la funcion CSC() se puede producir automaticamente modelos hazard de causa especifica

library (riskRegression)
CSH <-CSC(Hist(time,status)~age+invasion+sex, data=data_melanoma)
summary(CSH)

Con el modelo de regresion regresión ajustado, se puede predecir el riesgo individual con covariables dadas para lo cual instalamos el siguiente paquete

install.packages("pec")

Por ejemplo, quiero predecir el riesgo de un paciente masculino de 50 años de edad y que tiene un nivel de invasión 2.

library(prodlim)
library(cmprsk)
library(riskRegression)
library(pec)
CSH <-CSC(Hist(time,status)~age+invasion+sex, data=data_melanoma)
print(CSH)
pec:::predictEventProb(CSH,cause=1,newdata=data.frame(age=50,invasion=factor("level.1",levels=levels(data_melanoma$invasion)), sex=factor("Female",levels=levels(data_melanoma$sex))),time=c(1000,2000,3000))
#                    [,1]               [,2]          [,3]
#  [1,]     0.1276021    0.2390171    0.3360492

pec:::predictEventProb(CSH,cause=1,newdata=data.frame(age=50, invasion=factor("level.2", levels=levels(data_melanoma$invasion)),sex=factor("Male",levels=levels(data_melanoma$sex))),time=c(1000,2000,3000))
#                     [,1]             [,2]            [,3]
#  [1,]     0.3189477    0.5363142    0.6846546

Modelo Fine Grey

SHdata_melanoma <- FGR(Hist(time,status)~sex+age+thick+invasion,data=data_melanoma)
SHdata_melanoma
## 
## Right-censored response of a competing.risks model
## 
## No.Observations: 205 
## 
## Pattern:
##          
## Cause     event right.censored
##   1          57              0
##   2          14              0
##   unknown     0            134
## 
## 
## Fine-Gray model: analysis of cause 1 
## 
## Competing Risks Regression
## 
## Call:
## FGR(formula = Hist(time, status) ~ sex + age + thick + invasion, 
##     data = data_melanoma, cause = "1")
## 
##                    coef exp(coef) se(coef)    z p-value
## sexMale         0.53506      1.71  0.28077 1.91   0.057
## age             0.00322      1.00  0.00948 0.34   0.730
## thick           0.09548      1.10  0.05735 1.66   0.096
## invasionlevel.1 0.85756      2.36  0.36053 2.38   0.017
## invasionlevel.2 0.74593      2.11  0.63757 1.17   0.240
## 
##                 exp(coef) exp(-coef)  2.5% 97.5%
## sexMale              1.71      0.586 0.985  2.96
## age                  1.00      0.997 0.985  1.02
## thick                1.10      0.909 0.983  1.23
## invasionlevel.1      2.36      0.424 1.163  4.78
## invasionlevel.2      2.11      0.474 0.604  7.36
## 
## Num. cases = 205
## Pseudo Log-likelihood = -272 
## Pseudo likelihood ratio test = 27.6  on 5 df,
## 
## Convergence: TRUE

Cuando comparamos los HR en los dos modelos obtenemos HR 1.94 vs 1.87 para analisis de riesgos competitivos refeljando supuestos diferentes para riesgos que compiten

Usando la funcion crr en el paquete cmprsk podemos tambien obtener la subdistribucion Hazard Ratio

Subdistribucion por la variable ulcer

covulcer<-model.matrix(~ulcer,data=data_melanoma)[,-1]
crr.ulcer<-crr(data_melanoma$time,data_melanoma$status,cov1=covulcer)
summary(crr.ulcer)
## Competing Risks Regression
## 
## Call:
## crr(ftime = data_melanoma$time, fstatus = data_melanoma$status, 
##     cov1 = covulcer)
## 
##           coef exp(coef) se(coef)   z p-value
## covulcer1 1.41      4.11    0.288 4.9 9.5e-07
## 
##           exp(coef) exp(-coef) 2.5% 97.5%
## covulcer1      4.11      0.243 2.34  7.24
## 
## Num. cases = 205
## Pseudo Log-likelihood = -273 
## Pseudo likelihood ratio test = 26.2  on 1 df,

Los Hazard ratio de subdistribucion ulcer SHR: 4.11 con IC95% (2.34, 7.24) con un p-value significativo 9.5e-07

covthick<-model.matrix(~thick,data=data_melanoma)[,-1]
crr.thick<-crr(data_melanoma$time,data_melanoma$status,cov1=covthick)
summary(crr.thick)
## Competing Risks Regression
## 
## Call:
## crr(ftime = data_melanoma$time, fstatus = data_melanoma$status, 
##     cov1 = covthick)
## 
##            coef exp(coef) se(coef)    z p-value
## covthick1 0.149      1.16    0.028 5.32 1.1e-07
## 
##           exp(coef) exp(-coef) 2.5% 97.5%
## covthick1      1.16      0.862  1.1  1.23
## 
## Num. cases = 205
## Pseudo Log-likelihood = -277 
## Pseudo likelihood ratio test = 16.8  on 1 df,

Los Hazard ratio de subdistribucion thick SHR: 1.16 con IC95% (1.1, 1.23) con un p-value significativo 1.1e-07

covage<-model.matrix(~age,data=data_melanoma)[,-1]
crr.age<-crr(data_melanoma$time,data_melanoma$status,cov1=covage)
summary(crr.age)
## Competing Risks Regression
## 
## Call:
## crr(ftime = data_melanoma$time, fstatus = data_melanoma$status, 
##     cov1 = covage)
## 
##          coef exp(coef) se(coef)    z p-value
## covage1 0.015      1.02  0.00944 1.59    0.11
## 
##         exp(coef) exp(-coef)  2.5% 97.5%
## covage1      1.02      0.985 0.996  1.03
## 
## Num. cases = 205
## Pseudo Log-likelihood = -284 
## Pseudo likelihood ratio test = 3.18  on 1 df,

Los Hazard ratio de subdistribucion age SHR: 1.02 con IC95% (0,996, 1.03) con un p-value no significativo 0.11

covinvasion<-model.matrix(~invasion,data=data_melanoma)[,-1]
crr.invasion<-crr(data_melanoma$time,data_melanoma$status,cov1=covinvasion)
summary(crr.invasion)
## Competing Risks Regression
## 
## Call:
## crr(ftime = data_melanoma$time, fstatus = data_melanoma$status, 
##     cov1 = covinvasion)
## 
##                 coef exp(coef) se(coef)    z p-value
## invasionlevel.1 1.09      2.99    0.312 3.50 0.00046
## invasionlevel.2 1.39      4.03    0.399 3.49 0.00048
## 
##                 exp(coef) exp(-coef) 2.5% 97.5%
## invasionlevel.1      2.99      0.335 1.62  5.50
## invasionlevel.2      4.03      0.248 1.84  8.82
## 
## Num. cases = 205
## Pseudo Log-likelihood = -277 
## Pseudo likelihood ratio test = 18.3  on 2 df,

Los Hazard ratio de subdistribucion invasionlevel1 SHR:2.99 con IC95% (1,62 5,50) con un p-value significativo 0.00046 Los Hazard ratio de subdistribucion invasionlevel1 SHR:4.03 con IC95% (1,84 8,82) con un p-value significativo 0.00048

covsex<-model.matrix(~sex,data=data_melanoma)[,-1]
crr.sex<-crr(data_melanoma$time,data_melanoma$status,cov1=covsex)
summary(crr.sex)
## Competing Risks Regression
## 
## Call:
## crr(ftime = data_melanoma$time, fstatus = data_melanoma$status, 
##     cov1 = covsex)
## 
##          coef exp(coef) se(coef)    z p-value
## covsex1 0.637      1.89    0.263 2.42   0.016
## 
##         exp(coef) exp(-coef) 2.5% 97.5%
## covsex1      1.89      0.529 1.13  3.17
## 
## Num. cases = 205
## Pseudo Log-likelihood = -283 
## Pseudo likelihood ratio test = 5.69  on 1 df,

Los Hazard ratio de subdistribucion age SHR: 1.89 con IC95% (1.13 3.17) con un p-value significativo 0.016

Podemos escoger las variables significativas para obtener un modelo

cov<-model.matrix(~sex+ulcer+thick+invasion,data=data_melanoma)[,-1]
crr.modeldata_melanoma <- crr(data_melanoma$time,data_melanoma$status,cov1=cov)

Prediccion del Modelo

newdata_melanoma<-data.frame(sex=factor(c("Male","Male","Female"),levels=levels(data_melanoma$sex)),ulcer=factor(c("not present","present","not present"),levels=levels(data_melanoma$ulcer)),thick=c(0.5,6.1,3.4),invasion =factor(c("level.0","level.1","level.2"), levels=levels(data_melanoma$invasion)))
newdata_melanoma
##      sex       ulcer thick invasion
## 1   Male not present   0.5  level.0
## 2   Male     present   6.1  level.1
## 3 Female not present   3.4  level.2

Porque sex, ulcer e invasion son variables factor necesitan ser transformadas en varibales dummy con model.matrix()

dummydata_melanoma.new <- model.matrix(~sex+ulcer+thick+invasion,data=newdata_melanoma)[,-1]
dummydata_melanoma.new
##   sexMale ulcerpresent thick invasionlevel.1 invasionlevel.2
## 1       1            0   0.5               0               0
## 2       1            1   6.1               1               0
## 3       0            0   3.4               0               1

Prediccion con las nuevos datos

pred<-predict(crr.modeldata_melanoma,dummydata_melanoma.new)

Grafico prediccion

plot(pred,lty=1:3,col=1:3,xlab="Failure time(days)",ylab="Cumulative incidence function")
legend("topleft",c("Male,ulcer_not present,thick=0.5,invasion2","Male,ulcer_present,thick=6.1,invasion1","Female,ulcer_not present,thick=3.4,invasion2"),lty=1:3,col=1:3)

Para la predicción y el grafico se pueden usar tambien el paquete riskRegression. La función riskRegression () en presencia de riesgos competitivos

reg<-riskRegression(Hist(time, status) ~ sex+ulcer+thick+invasion, data = data_melanoma, cause = 1,link="prop") 
plot(reg,newdata=newdata_melanoma)
legend("topleft",c("Male,ulcer_not present,thick=0.5,invasion2","Male,ulcer_present,thick=6.1,invasion1","Female,ulcer_not present,thick=3.4,invasion2"),lty=1:3,col=1:3)

Diagnostico del modelo

Se debe precisar el cumplimiento del suspuesto de proporcionalidad. En el cual la funcion de riesgo es proporcional a dos factores pronosticos distintos y se mantienen a lo largo del tiempo.

La incidencia acumulada funciona a diferentes niveles de invasión, estableciendo el sexo en masculino, ulcera presente y thick en 0.5 mm. No hay evidencia de violación del supuesto de proporcionalidad para la variable invasión .

checkdata_melanoma <- data.frame(sex=factor(c("Male","Male","Male"), levels=levels(data_melanoma$sex)), ulcer=factor(c("present","present","present"),levels=levels(data_melanoma$ulcer)), thick=c(0.5,0.5,0.5), invasion=factor(c("level.0","level.1","level.2"), levels=levels(data_melanoma$invasion))) 
plot(reg,newdata=checkdata_melanoma,lty=1:3,col=1:3)
text(2200,1,"Covariates sex='Male'; ulcer='present'; thick=0.5")
legend("topleft",c("invasion.level0","invasion. level1","invasion.level2"),lty=1:3,col=1:3)

En el modelo de regresion un metodo diferente para comprobar el supuesto de proporcionalidad es incluir una covariable dependiente del tiempo

crr.time<-crr(data_melanoma$time,data_melanoma$status,cov1=cov,cov2=cov[,1],tf=function(t) t)
summary(crr.time)
## Competing Risks Regression
## 
## Call:
## crr(ftime = data_melanoma$time, fstatus = data_melanoma$status, 
##     cov1 = cov, cov2 = cov[, 1], tf = function(t) t)
## 
##                      coef exp(coef) se(coef)      z p-value
## sexMale          0.949401      2.58 0.504025  1.884  0.0600
## ulcerpresent     1.006569      2.74 0.321485  3.131  0.0017
## thick            0.055870      1.06 0.065175  0.857  0.3900
## invasionlevel.1  0.635717      1.89 0.372436  1.707  0.0880
## invasionlevel.2  0.678247      1.97 0.619732  1.094  0.2700
## cov[, 1]1*tf1   -0.000387      1.00 0.000349 -1.107  0.2700
## 
##                 exp(coef) exp(-coef)  2.5% 97.5%
## sexMale              2.58      0.387 0.962  6.94
## ulcerpresent         2.74      0.365 1.457  5.14
## thick                1.06      0.946 0.931  1.20
## invasionlevel.1      1.89      0.530 0.910  3.92
## invasionlevel.2      1.97      0.508 0.585  6.64
## cov[, 1]1*tf1        1.00      1.000 0.999  1.00
## 
## Num. cases = 205
## Pseudo Log-likelihood = -266 
## Pseudo likelihood ratio test = 39.5  on 6 df,

En el modelo Fine-Gray la covariable sexo variable en el tiempo no muestra significación estadística (P = 0.27), lo que indica que el efecto del sexo es constante en el tiempo

La librería RiskRegression proporciona una solución para modelar una covariable que varia en el tiempo.

reg.time<-riskRegression(Hist(time, status) ~ sex + ulcer + thick + strata(invasion), data = data_melanoma, cause = 1,link="prop")
plotEffects(reg.time,formula=~invasion)

Se grafican los efectos dependientes del tiempo en el modelo de regresión de Fine-Gray para la mortalidad. La curva y el intervalo de confianza del 95% correspondiente se dibujan con un método no paramétrico

Los residuales de Schoenfeld son otro metodo para verificar el modelo

par(mfrow=c(2,2))
for(j in 1:ncol(crr.modeldata_melanoma$res)) {scatter.smooth(crr.modeldata_melanoma$uft, crr.modeldata_melanoma$res[,j], main=names(crr.modeldata_melanoma$coef)[j], xlab = "Failure time", ylab ="Schoenfeld residuals") }

El nivel de invasion 1 tiene una pendiente diferente de 0, lo cual indica residuales no son constantes a traves del tiempo ,indicando que no se cumple el supuesto de proporcionalidad