La regresión de Cox se utiliza en el análisis del tiempo de supervivencia para determinar la influencia de distintas variables en el tiempo de supervivencia. Las variables pueden ser cualquier mezcla de datos continuos, binarios o categóricos. A continuación, se utiliza el modelo de riesgos proporcionales de Cox para determinar el efecto sobre el tiempo de supervivencia. La regresión de Cox nos permite determinar los efectos de múltiples variables independientes sobre un resultado en el tiempo, ya sea para probar hipótesis sobre las variables independientes o para construir un modelo predictivo.

framin<- read.csv("framingham_dataset.csv", header = TRUE)
library("survival")
library("survminer")
colnames(framin)
##  [1] "X"        "randid"   "sex"      "totchol"  "age"      "sysbp"   
##  [7] "diabp"    "cursmoke" "cigpday"  "bmi"      "diabetes" "bpmeds"  
## [13] "heartrte" "glucose"  "educ"     "prevchd"  "prevap"   "prevmi"  
## [19] "prevstrk" "prevhyp"  "time"     "period"   "hdlc"     "ldlc"    
## [25] "death"    "angina"   "hospmi"   "mi_fchd"  "anychd"   "stroke"  
## [31] "cvd"      "hyperten" "timeap"   "timemi"   "timemifc" "timechd" 
## [37] "timestrk" "timecvd"  "timedth"  "timehyp"

Las variables mas importantes en este caso serian la variable ‘death’ que es el estatus del paciente (vivo o muerto) y el ‘timedth’ que nos indica cuantos dias pasaron desde el examen inicial hasta la muerte o de ultimo contacto con el paciente y como factores tomaremos el sexo, la edad, colesterol y la hipertension

res.cox <- coxph(Surv(timedth, death) ~ sex + age + totchol + hyperten, data = framin)
res.cox
## Call:
## coxph(formula = Surv(timedth, death) ~ sex + age + totchol + 
##     hyperten, data = framin)
## 
##                coef  exp(coef)   se(coef)       z      p
## sex      -0.6312757  0.5319128  0.0353980 -17.834 <2e-16
## age       0.0586254  1.0603779  0.0018443  31.787 <2e-16
## totchol   0.0010891  1.0010897  0.0003931   2.770 0.0056
## hyperten  0.0943203  1.0989117  0.0430538   2.191 0.0285
## 
## Likelihood ratio test=1358  on 4 df, p=< 2.2e-16
## n= 11218, number of events= 3410 
##    (409 observations deleted due to missingness)

En esta regresion, los datos nos indican que en la variable sex, ser mujer tiene un menor riesgo en comparacion con los hombres, en la edad el coeficiente positivo indica que con mayor edad mayor riesgo por cada año aumenta 1.06 el HR (Hazard ratio o taza de riesgo) dando asi la edad como predictor muy significativo estadisticamente, al igual que la edad el colesterol conforme mayor sea mayor sera el riesgon este caso por cada unidad de colesterol el riesgo aumenta 1% aunque pequeño es el efecto el p-value lo hace estadisticamente significativo y finalmente la hipertension es otro predictor significativo en comparacion de las personas que no tienen la enfermedad, aunque todas las variables son significativas las mas fuertes son la edad y el sexo. Ahora se realizara la misma tarea pero con datos imputados con el metodo MICE.

res.cox <- coxph(Surv(timedth, death) ~ sex + age + totchol + hyperten, data = framin_comp)
res.cox
## Call:
## coxph(formula = Surv(timedth, death) ~ sex + age + totchol + 
##     hyperten, data = framin_comp)
## 
##                coef  exp(coef)   se(coef)       z       p
## sex      -0.6268996  0.5342456  0.0347500 -18.040 < 2e-16
## age       0.0592295  1.0610187  0.0018149  32.636 < 2e-16
## totchol   0.0011454  1.0011460  0.0003838   2.984 0.00284
## hyperten  0.1098399  1.1160994  0.0424926   2.585 0.00974
## 
## Likelihood ratio test=1426  on 4 df, p=< 2.2e-16
## n= 11627, number of events= 3527

Para este caso todos los coeficientes aumentaron, pero el aumento mas notable se vio en la hipertension lo cual aumenta nuestro HR mas que en los datos no imputados, pero las lecturas sobre los efectos que tienen estas variables como predictor se mantienen iguales.

covariates <- c("age", "sex",  "totchol", "hyperten")
univ_formulas <- sapply(covariates,
                        function(x) as.formula(paste('Surv(timedth, death)~', x)))
                        
univ_models <- lapply( univ_formulas, function(x){coxph(x, data = framin_comp)})
# Extract data 
univ_results <- lapply(univ_models,
                       function(x){ 
                          x <- summary(x)
                          p.value<-signif(x$wald["pvalue"], digits=2)
                          wald.test<-signif(x$wald["test"], digits=2)
                          beta<-signif(x$coef[1], digits=2);#coeficient beta
                          HR <-signif(x$coef[2], digits=2);#exp(beta)
                          HR.confint.lower <- signif(x$conf.int[,"lower .95"], 2)
                          HR.confint.upper <- signif(x$conf.int[,"upper .95"],2)
                          HR <- paste0(HR, " (", 
                                       HR.confint.lower, "-", HR.confint.upper, ")")
                          res<-c(beta, HR, wald.test, p.value)
                          names(res)<-c("beta", "HR (95% CI for HR)", "wald.test", 
                                        "p.value")
                          return(res)
                          #return(exp(cbind(coef(x),confint(x))))
                         })
res <- t(as.data.frame(univ_results, check.names = FALSE))
as.data.frame(res)
##            beta HR (95% CI for HR) wald.test  p.value
## age       0.058      1.1 (1.1-1.1)      1100 2.7e-237
## sex       -0.54   0.58 (0.55-0.62)       250  9.6e-57
## totchol  0.0014            1 (1-1)        16  8.3e-05
## hyperten   0.33      1.4 (1.3-1.5)        61  5.2e-15