Riesgos Competitivos
Omar Huang David Escobar Franco
6/9/2019
Para crear un R markdown (Rmd) lo primero que vamos a hacer es uprimir todas las alertas en el scrip con el siguiente codigo
options(warn=-1)Instalacion de paquetes
Vamos a instalar los siquientes paquetes, dentro del Rmd
install.packages("cmprsk")
install.packages("rpivotTable")
install.packages("mstate")
install.packages("gridExtra")abrimos las librerias
library(survival)
library(cmprsk)
library(mstate)y abrimos la base de datos de melanoma que estan basadas en el articulo melanoma la cual esta en stata (.dta)
library(haven)
melanoma <- read_dta("~/Documents/Maestria Epidemiologia/Bioestadistica IV/Competencia de Riesgos/malignantmelanoma.dta")con el siguiente comando vemos la estructura que compone la base de datos, sus detalles, por ejemplo
str(melanoma)## Classes 'tbl_df', 'tbl' and 'data.frame': 205 obs. of 5 variables:
## $ id : num 1 2 3 4 5 6 7 8 9 10 ...
## ..- attr(*, "format.stata")= chr "%8.0g"
## $ time : num 10 30 35 99 185 204 210 232 232 279 ...
## ..- attr(*, "format.stata")= chr "%8.0g"
## $ thick: num 3.84 -2.27 -1.58 -0.02 9.16 1.92 2.24 9.96 0.3 4.49 ...
## ..- attr(*, "label")= chr "Thickness of the tumor in mm minus the mean thickness"
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ sex : 'haven_labelled' num 1 1 1 0 1 1 1 1 0 0 ...
## ..- attr(*, "format.stata")= chr "%8.0g"
## ..- attr(*, "labels")= Named num 0 1
## .. ..- attr(*, "names")= chr "F" "M"
## $ cause: 'haven_labelled' num 2 2 0 2 1 1 1 1 2 1 ...
## ..- attr(*, "format.stata")= chr "%18.0g"
## ..- attr(*, "labels")= Named num 0 1 2
## .. ..- attr(*, "names")= chr "Alive" "Malignant Melanoma" "Other Causes"donde encontramos que la base de datos se compone de 5 variables, causa tiene 3 niveles y sexo 2, nos indica cada nivel como esta codificado
crear un factor
Los factores se componen de dos vectores osea 2D, y se crean de la siguiente manera.
Estos factores le vamos a dar el normbre que queramos y el nombre a cada nivel.
melanoma$causa <-factor(melanoma$cause, levels=c(0,1,2), labels=c("vivo", "Melanoma Maligno","otra causa" ))
melanoma$sexo <-factor(melanoma$sex, levels=c(0,1), labels=c("Famele", "Male"))Descripciones basicas y resumenes
La siguiente tabla mostrara el numero de observaciones por cada combinacion de evento y enfermedad, al igual que se puede subdividir en estrato como el sexo.
attach(melanoma)
table(causa)## causa
## vivo Melanoma Maligno otra causa
## 134 57 14table(causa, sexo)## sexo
## causa Famele Male
## vivo 91 43
## Melanoma Maligno 28 29
## otra causa 7 7La opcion average, determina el promedio. ### la media de tiempo
tapply(time, causa, mean)## vivo Melanoma Maligno otra causa
## 2620.672 1252.947 1338.286la media de tiempo por sexo
tapply(time, list(causa, sexo), mean)## Famele Male
## vivo 2641.011 2577.628
## Melanoma Maligno 1382.179 1128.172
## otra causa 1225.714 1450.857o se puede crear la tabla con la siguiente libreria
library(rpivotTable)
rpivotTable(data=melanoma, rows="causa",cols = "sexo", vals = "time", aggregatorName = "Average",rendererName = "Table" )round(tapply(time, list(causa, sexo), mean), digits = 2)## Famele Male
## vivo 2641.01 2577.63
## Melanoma Maligno 1382.18 1128.17
## otra causa 1225.71 1450.86en esta table podemos modificar cuales variables escoger,
Modelo de COX para riesgo competitivo
EL modelo cox para riesgo sproporcionales noc ayudara a calclcular el HR sin tener en cuenta las muertes por otras causas.
Crear el objeto
en los siguientes comando vamos a crear los objetos por las variables de interes, en este caso son: sexo y grosor
por sexo
Coxsex<- coxph(Surv(time, causa == "Melanoma Maligno") ~ sex, data = melanoma)
summary(Coxsex)## Call:
## coxph(formula = Surv(time, causa == "Melanoma Maligno") ~ sex,
## data = melanoma)
##
## n= 205, number of events= 57
##
## coef exp(coef) se(coef) z Pr(>|z|)
## sex 0.6622 1.9390 0.2651 2.498 0.0125 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## sex 1.939 0.5157 1.153 3.26
##
## Concordance= 0.59 (se = 0.034 )
## Likelihood ratio test= 6.15 on 1 df, p=0.01
## Wald test = 6.24 on 1 df, p=0.01
## Score (logrank) test = 6.47 on 1 df, p=0.01EL HR de sexo en el modelo de cox es de 1.94 con CI95%(1.153, 3.26), siendo significativo (pvalue= 0.01) para el modelo final (cox), y podemos intrerpretar que los hombres tienen 0.94 veces mas riesgo que las mujeres.
por grosor
Coxthick <- coxph(Surv(time, causa == "Melanoma Maligno") ~ thick, data = melanoma)
summary(Coxthick)## Call:
## coxph(formula = Surv(time, causa == "Melanoma Maligno") ~ thick,
## data = melanoma)
##
## n= 205, number of events= 57
##
## coef exp(coef) se(coef) z Pr(>|z|)
## thick 0.16024 1.17380 0.03126 5.126 2.96e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## thick 1.174 0.8519 1.104 1.248
##
## Concordance= 0.741 (se = 0.032 )
## Likelihood ratio test= 19.19 on 1 df, p=1e-05
## Wald test = 26.27 on 1 df, p=3e-07
## Score (logrank) test = 28.7 on 1 df, p=8e-08EL HR del grosor en el modelo de cox es de 1.174 con CI95%(1.10, 1,25), siendo significativo (pvalue= 2.96e-07) para el modelo final (cox), y podemos intrerpretar que el grosor que esta dado en milimetros (mm) tienen 0.17 veces mas riesgo por cada unidad.
modelando con las variables significativas
Si aguegamos al amodelo las dos variables se ve levemente alterado los valores, debemos ver si continuan siendo significativos.
Coxsexthick <- coxph(Surv(time, causa == "Melanoma Maligno") ~ sex + thick, data = melanoma)
summary(Coxsexthick)## Call:
## coxph(formula = Surv(time, causa == "Melanoma Maligno") ~ sex +
## thick, data = melanoma)
##
## n= 205, number of events= 57
##
## coef exp(coef) se(coef) z Pr(>|z|)
## sex 0.57411 1.77556 0.26526 2.164 0.0304 *
## thick 0.15910 1.17245 0.03268 4.869 1.12e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## sex 1.776 0.5632 1.056 2.986
## thick 1.172 0.8529 1.100 1.250
##
## Concordance= 0.719 (se = 0.035 )
## Likelihood ratio test= 23.82 on 2 df, p=7e-06
## Wald test = 28.77 on 2 df, p=6e-07
## Score (logrank) test = 32.2 on 2 df, p=1e-07ambos son significativos, podemos observar los Beta del modelo y sus respectivos HR con intervalos de confianza.
Graficar KM
creamos el objeto Surv y lo ajustamos antes de graficarlo por sexo
library(KMsurv)
melanoma.km <- survfit(Surv(melanoma$time, melanoma$cause) ~ 1, data = melanoma)
head(melanoma.km)## Call: survfit(formula = Surv(melanoma$time, melanoma$cause) ~ 1, data = melanoma)
##
## 134 observations deleted due to missingness
## n events median 0.95LCL 0.95UCL
## 71 14 3182 3154 NAes de notar la advertencia del status donde nos avisa que tiene mas de 2 eventos, luego creamos la grafica
library(ggplot2)
library(survminer)## Loading required package: ggpubr## Loading required package: magrittrggsurvplot(fit = melanoma.km, data = melanoma, conf.int = T, title = "Curva de Supervivencia",
xlab = "Tiempo", ylab = "Probabilidad de supervivencia", legend.title = "Estimación",
legend.labs = "Kaplan-Meier")
LogRank por sex
tambien podemos graficar el log rank de KM por sexo
melanoma.sex <- survfit(Surv(melanoma$time, melanoma$cause) ~ sexo, data = melanoma, conf.type = "log-log")
ggsurvplot(melanoma.sex, fun = "event", censor = F, xlab = "Time (dias)")
comparado con la creada en stata
Cumulative Incidence
las siguienes son las librerias que vamos a usar
library(cmprsk)
library(survival)
library(mstate)
library(ggplot2)Ajustando la funcion acumulada
pimero creamos la funcion de Incidencia acumulada por sexo, y luego con el Test corroboramos el valor P
cif <- cuminc(ftime = melanoma$time, fstatus= melanoma$causa, group = melanoma$sexo)
cifsex <- cuminc(ftime = melanoma$time, fstatus= melanoma$cause, group = melanoma$sex)
cifsex$Tests## stat pv df
## 1 5.8140209 0.0158989 1
## 2 0.8543656 0.3553203 1observamos el valor p de la diferencia de incidencias acumuladas por grupo de sexo,
CIF <- Cuminc(time = "time", status = "causa", data = melanoma)
head(CIF)## time Surv CI.Melanoma Maligno CI.otra causa seSurv
## 1 10 0.9951220 0.000000000 0.004878049 0.004866137
## 2 30 0.9902439 0.000000000 0.009756098 0.006864869
## 3 99 0.9853417 0.000000000 0.014658295 0.008400806
## 4 185 0.9804395 0.004902198 0.014658295 0.009684268
## 5 204 0.9755373 0.009804395 0.014658295 0.010805597
## 6 210 0.9706351 0.014706593 0.014658295 0.011811062
## seCI.Melanoma Maligno seCI.otra causa
## 1 0.000000000 0.004866137
## 2 0.000000000 0.006864869
## 3 0.000000000 0.008400806
## 4 0.004890166 0.008400806
## 5 0.006898683 0.008400806
## 6 0.008428185 0.008400806ggplot(CIF, aes(time)) +
geom_step(aes(y = `CI.Melanoma Maligno`, colour = "Muerte por melanoma"))+ geom_step(aes(y = `CI.Melanoma Maligno`+ `CI.otra causa`, colour = "Total")) + geom_step(aes(y = Surv, colour = "Overall survival")) + labs(x = "Tiempo (años)", y = "Proporcion", colour = "") + theme_classic()
para graficar las incidencias acumuladas debemos crear el siguiente factor
cif_sex <- Cuminc(time = "time", status = "causa", group = melanoma$sexo, data= melanoma )
head(cifsex)## Tests:
## stat pv df
## 1 5.8140209 0.0158989 1
## 2 0.8543656 0.3553203 1
## Estimates and Variances:
## $est
## 1000 2000 3000 4000 5000
## 0 1 0.08730159 0.18077594 0.23565169 0.28424490 0.28424490
## 1 1 0.19237175 0.31009828 0.42453587 0.42453587 NA
## 0 2 0.03174603 0.03983516 0.05220642 0.08538385 0.08538385
## 1 2 0.03814124 0.06693942 0.06693942 0.13474271 NA
##
## $var
## 1000 2000 3000 4000 5000
## 0 1 0.0006378135 0.0012450462 0.0018102025 0.002755577 0.002755577
## 1 1 0.0020223293 0.0028196248 0.0042695603 0.004269560 NA
## 0 2 0.0002459647 0.0003073878 0.0004529114 0.001528480 0.001528480
## 1 2 0.0004727379 0.0008614343 0.0008614343 0.002950698 NAlibrary(gridExtra)
grid.arrange(
ggplot(cif_sex, aes(time)) +
geom_step(aes(y = `CI.Melanoma Maligno`, colour = group)) +
labs(x = "Time (dias)", y = "Proporcion", title = "Muerte por Melanoma") +
theme_classic(),
ggplot(cif_sex, aes(time)) +
geom_step(aes(y = `CI.otra causa`, colour = group)) +
labs(x = "Time (dias)", y = "Proporcion", title = "Muerte por otra causa") +
theme_classic(),
ncol = 2
)
cif$Tests## stat pv df
## vivo 3.9808203 0.04602114 1
## Melanoma Maligno 5.6905837 0.01705618 1
## otra causa 0.8264053 0.36331405 1Graficar incidencia acumulada por evento
Para crear unos graficos con mas detalles descargamos el siguiente paquete CumIncidence
source("~/Documents/Maestria Epidemiologia/Bioestadistica IV/Competencia de Riesgos/CumIncidence.R")
fit<-CumIncidence (melanoma$time, melanoma$sexo, melanoma$causa, cencode = 0, xlab = "dias")##
## +-------------------------------------------------------------------+
## | Cumulative incidence function estimates from competing risks data |
## +-------------------------------------------------------------------+
## Test equality across groups:
## Statistic p-value df
## Famele 0.3539 0.8378091 2
## Male 15.5665 0.0004166 2
##
## Estimates at time points:
## 0 1000 2000 3000 4000 5000
## vivo Famele 0 0.000000 0.2313 0.4328 0.6119 0.6716
## Melanoma Maligno Famele 0 0.192982 0.3860 0.4561 NA NA
## otra causa Famele 0 0.285714 0.3571 0.4286 NA NA
## vivo Male 0 0.007463 0.1119 0.2015 0.2910 0.3209
## Melanoma Maligno Male 0 0.263158 0.4211 0.5088 NA NA
## otra causa Male 0 0.214286 0.3571 0.3571 NA NA
##
## Standard errors:
## 0 1000 2000 3000 4000 5000
## vivo Famele 0 0.000000 0.03659 0.04306 0.04255 0.04152
## Melanoma Maligno Famele 0 0.052928 0.06585 0.06799 NA NA
## otra causa Famele 0 0.127385 0.13640 0.14338 NA NA
## vivo Male 0 0.007463 0.02736 0.03485 0.03961 0.04094
## Melanoma Maligno Male 0 0.059005 0.06656 0.06877 NA NA
## otra causa Male 0 0.114884 0.13732 0.13732 NA NA
col<-c(2,3,4)
fit<-CumIncidence (melanoma$time, melanoma$sexo, melanoma$causa, cencode = 0, xlab = "dias", col= col)##
## +-------------------------------------------------------------------+
## | Cumulative incidence function estimates from competing risks data |
## +-------------------------------------------------------------------+
## Test equality across groups:
## Statistic p-value df
## Famele 0.3539 0.8378091 2
## Male 15.5665 0.0004166 2
##
## Estimates at time points:
## 0 1000 2000 3000 4000 5000
## vivo Famele 0 0.000000 0.2313 0.4328 0.6119 0.6716
## Melanoma Maligno Famele 0 0.192982 0.3860 0.4561 NA NA
## otra causa Famele 0 0.285714 0.3571 0.4286 NA NA
## vivo Male 0 0.007463 0.1119 0.2015 0.2910 0.3209
## Melanoma Maligno Male 0 0.263158 0.4211 0.5088 NA NA
## otra causa Male 0 0.214286 0.3571 0.3571 NA NA
##
## Standard errors:
## 0 1000 2000 3000 4000 5000
## vivo Famele 0 0.000000 0.03659 0.04306 0.04255 0.04152
## Melanoma Maligno Famele 0 0.052928 0.06585 0.06799 NA NA
## otra causa Famele 0 0.127385 0.13640 0.14338 NA NA
## vivo Male 0 0.007463 0.02736 0.03485 0.03961 0.04094
## Melanoma Maligno Male 0 0.059005 0.06656 0.06877 NA NA
## otra causa Male 0 0.114884 0.13732 0.13732 NA NA
fit<-CumIncidence (melanoma$time, melanoma$sexo, melanoma$causa, cencode = 0, xlab = "dias", t=c(0, 3, 6, 9, 12, 24, 36, 48), col=col)##
## +-------------------------------------------------------------------+
## | Cumulative incidence function estimates from competing risks data |
## +-------------------------------------------------------------------+
## Test equality across groups:
## Statistic p-value df
## Famele 0.3539 0.8378091 2
## Male 15.5665 0.0004166 2
##
## Estimates at time points:
## 0 3 6 9 12 24 36 48
## vivo Famele 0 0 0 0 0.00000 0.00000 0.000000 0.000000
## Melanoma Maligno Famele 0 0 0 0 0.00000 0.00000 0.000000 0.000000
## otra causa Famele 0 0 0 0 0.00000 0.00000 0.000000 0.000000
## vivo Male 0 0 0 0 0.00000 0.00000 0.007463 0.007463
## Melanoma Maligno Male 0 0 0 0 0.00000 0.00000 0.000000 0.000000
## otra causa Male 0 0 0 0 0.07143 0.07143 0.142857 0.142857
##
## Standard errors:
## 0 3 6 9 12 24 36 48
## vivo Famele 0 0 0 0 0.00000 0.00000 0.000000 0.000000
## Melanoma Maligno Famele 0 0 0 0 0.00000 0.00000 0.000000 0.000000
## otra causa Famele 0 0 0 0 0.00000 0.00000 0.000000 0.000000
## vivo Male 0 0 0 0 0.00000 0.00000 0.007463 0.007463
## Melanoma Maligno Male 0 0 0 0 0.00000 0.00000 0.000000 0.000000
## otra causa Male 0 0 0 0 0.07143 0.07143 0.097208 0.097208
# # The resulting output shows both estimates of cumulative incidence and s.e. # evaluated at the given time points.
# # The CumIncidence() function allows for the pointwise confidence intervals, # by simply adding a further argument, level, where we specify the desired # confidence level. For example, we may compute 95% pointwise confidence # interval at our selected time points as follows: #
Para agregar los intervalos de confianza
par(mfrow=c(2,3))
fit<-CumIncidence (melanoma$time, melanoma$sexo, melanoma$causa, cencode = 0, xlab = "dias", col= col,level = 0.95)##
## +-------------------------------------------------------------------+
## | Cumulative incidence function estimates from competing risks data |
## +-------------------------------------------------------------------+
## Test equality across groups:
## Statistic p-value df
## Famele 0.3539 0.8378091 2
## Male 15.5665 0.0004166 2
##
## Estimates at time points:
## 0 1000 2000 3000 4000 5000
## vivo Famele 0 0.000000 0.2313 0.4328 0.6119 0.6716
## Melanoma Maligno Famele 0 0.192982 0.3860 0.4561 NA NA
## otra causa Famele 0 0.285714 0.3571 0.4286 NA NA
## vivo Male 0 0.007463 0.1119 0.2015 0.2910 0.3209
## Melanoma Maligno Male 0 0.263158 0.4211 0.5088 NA NA
## otra causa Male 0 0.214286 0.3571 0.3571 NA NA
##
## Standard errors:
## 0 1000 2000 3000 4000 5000
## vivo Famele 0 0.000000 0.03659 0.04306 0.04255 0.04152
## Melanoma Maligno Famele 0 0.052928 0.06585 0.06799 NA NA
## otra causa Famele 0 0.127385 0.13640 0.14338 NA NA
## vivo Male 0 0.007463 0.02736 0.03485 0.03961 0.04094
## Melanoma Maligno Male 0 0.059005 0.06656 0.06877 NA NA
## otra causa Male 0 0.114884 0.13732 0.13732 NA NA
##
## 95% pointwise confidence intervals:
##
## , , vivo Famele
##
## 0 1000 2000 3000 4000 5000
## lower NaN NaN 0.1638 0.3475 0.5230 0.5829
## upper NaN NaN 0.3059 0.5151 0.6893 0.7456
##
## , , Melanoma Maligno Famele
##
## 0 1000 2000 3000 4000 5000
## lower NaN 0.1022 0.2586 0.3202 NA NA
## upper NaN 0.3053 0.5117 0.5822 NA NA
##
## , , otra causa Famele
##
## 0 1000 2000 3000 4000 5000
## lower NaN 0.08074 0.1188 0.1593 NA NA
## upper NaN 0.53599 0.6079 0.6765 NA NA
##
## , , vivo Male
##
## 0 1000 2000 3000 4000 5000
## lower NaN 0.0006702 0.06553 0.1381 0.2161 0.2426
## upper NaN 0.0375319 0.17212 0.2735 0.3700 0.4016
##
## , , Melanoma Maligno Male
##
## 0 1000 2000 3000 4000 5000
## lower NaN 0.1564 0.2901 0.3678 NA NA
## upper NaN 0.3827 0.5463 0.6334 NA NA
##
## , , otra causa Male
##
## 0 1000 2000 3000 4000 5000
## lower NaN 0.04749 0.1176 0.1176 NA NA
## upper NaN 0.45898 0.6094 0.6094 NA NApar(mfrow=c(1,1))
fit<-CumIncidence (melanoma$time, melanoma$sexo, melanoma$causa, cencode = 0, xlab = "dias", col= col,level = 0.95)##
## +-------------------------------------------------------------------+
## | Cumulative incidence function estimates from competing risks data |
## +-------------------------------------------------------------------+
## Test equality across groups:
## Statistic p-value df
## Famele 0.3539 0.8378091 2
## Male 15.5665 0.0004166 2
##
## Estimates at time points:
## 0 1000 2000 3000 4000 5000
## vivo Famele 0 0.000000 0.2313 0.4328 0.6119 0.6716
## Melanoma Maligno Famele 0 0.192982 0.3860 0.4561 NA NA
## otra causa Famele 0 0.285714 0.3571 0.4286 NA NA
## vivo Male 0 0.007463 0.1119 0.2015 0.2910 0.3209
## Melanoma Maligno Male 0 0.263158 0.4211 0.5088 NA NA
## otra causa Male 0 0.214286 0.3571 0.3571 NA NA
##
## Standard errors:
## 0 1000 2000 3000 4000 5000
## vivo Famele 0 0.000000 0.03659 0.04306 0.04255 0.04152
## Melanoma Maligno Famele 0 0.052928 0.06585 0.06799 NA NA
## otra causa Famele 0 0.127385 0.13640 0.14338 NA NA
## vivo Male 0 0.007463 0.02736 0.03485 0.03961 0.04094
## Melanoma Maligno Male 0 0.059005 0.06656 0.06877 NA NA
## otra causa Male 0 0.114884 0.13732 0.13732 NA NA
##
## 95% pointwise confidence intervals:
##
## , , vivo Famele
##
## 0 1000 2000 3000 4000 5000
## lower NaN NaN 0.1638 0.3475 0.5230 0.5829
## upper NaN NaN 0.3059 0.5151 0.6893 0.7456
##
## , , Melanoma Maligno Famele
##
## 0 1000 2000 3000 4000 5000
## lower NaN 0.1022 0.2586 0.3202 NA NA
## upper NaN 0.3053 0.5117 0.5822 NA NA
##
## , , otra causa Famele
##
## 0 1000 2000 3000 4000 5000
## lower NaN 0.08074 0.1188 0.1593 NA NA
## upper NaN 0.53599 0.6079 0.6765 NA NA
##
## , , vivo Male
##
## 0 1000 2000 3000 4000 5000
## lower NaN 0.0006702 0.06553 0.1381 0.2161 0.2426
## upper NaN 0.0375319 0.17212 0.2735 0.3700 0.4016
##
## , , Melanoma Maligno Male
##
## 0 1000 2000 3000 4000 5000
## lower NaN 0.1564 0.2901 0.3678 NA NA
## upper NaN 0.3827 0.5463 0.6334 NA NA
##
## , , otra causa Male
##
## 0 1000 2000 3000 4000 5000
## lower NaN 0.04749 0.1176 0.1176 NA NA
## upper NaN 0.45898 0.6094 0.6094 NA NAModelo de Fine-Grey
Evaluar efecto de las covariables curvas de incidencia acumulada y probar la diferencia entre ambos por la prueba de log-rank modificada (Gray)
es necesario instalar estos dos paquetes
install.packages("riskRegression")
install.packages("prodlim")las siguientes librerias vamos a usar
library (riskRegression)## Loading required package: data.table## Loading required package: prodlim## riskRegression version 2019.01.29library(prodlim)En modelo de riesgo proporcionales de FYG
SHR (Hazard Ratio Subdistribucion)
ahora para riesgos de subdistribucion por sexo
covsex<-model.matrix(~sex,data=melanoma)[,-1]
crr.sex<-crr(melanoma$time,melanoma$cause,cov1=covsex)
summary(crr.sex)## Competing Risks Regression
##
## Call:
## crr(ftime = melanoma$time, fstatus = melanoma$cause, 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 Hazar ratio de subdistribucion Sex SHR: 1.89 con IC95% (1.13, 3.17) con un p-value significativo (0.016).
ahora para riesgos de subdistribucion por grosor
covthick<-model.matrix(~thick,data=melanoma)[,-1]
crr.thick<-crr(melanoma$time,melanoma$cause,cov1=covthick)
summary(crr.thick)## Competing Risks Regression
##
## Call:
## crr(ftime = melanoma$time, fstatus = melanoma$cause, 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 Hazar ratio de subdistribucion Thick SHR: 1.16 con IC95% (1.1, 1.23) con un p-value significativo (1.1e-07)
Modelo final con las dos variables
cov<-model.matrix(~sex+thick,data=melanoma)[,-1]
crr.model<-crr(melanoma$time,melanoma$cause,cov1=cov)
summary(crr.model)## Competing Risks Regression
##
## Call:
## crr(ftime = melanoma$time, fstatus = melanoma$cause, cov1 = cov)
##
## coef exp(coef) se(coef) z p-value
## sex 0.515 1.67 0.2702 1.90 5.7e-02
## thick 0.144 1.15 0.0288 4.99 6.2e-07
##
## exp(coef) exp(-coef) 2.5% 97.5%
## sex 1.67 0.598 0.985 2.84
## thick 1.15 0.866 1.091 1.22
##
## Num. cases = 205
## Pseudo Log-likelihood = -276
## Pseudo likelihood ratio test = 20.5 on 2 df,Y el final de el modelo por riesgo proporcionales con las dos veriables podemos encontrar los beta respectivos para sexo y grosor, ambos con un valor de p significativos pero observamos que en el modelo juntos el intervalo de confianza de sexo pasa por el 1. por lo cual podemos considerar que para un modelo predictivo pierde sustento esta variable.
Es de anotar que de aqui para adelante se puede crear un modelo predictivo como en el articulo mensionado, sin emgargo no tenemos todas las variables que usaron.