Cox analysis

Objetivos

OBJETIVO GENERAL

Determinar los factores que intervienen en el pronóstico de sobrevida pulpar de dientes permanentes con desarrollo radicular incompleto que sufren algún tipo de fractura coronaria.

OBJETIVOS ESPECÍFICOS

• Determinar si la presencia de malos hábitos, traumatismos repetidos y control inadecuado de la higiene (por sÍ solos o asociados) influyen en el pronóstico de sobrevida pulpar de dientes permanentes con desarrollo radicular incompleto que sufren fracturas coronarias.
• Obtener un orden de relevancia de los factores que intervienen en el pronóstico de sobrevida de dientes permanentes con desarrollo radicular incompleto que sufren fracturas coronarias.
• Generar un protocolo de abordaje de dientes permanentes con desarrollo radicular incompleto que sufren fracturas coronarias.

Paquetes

Paquetes especÍficos

Dataset

df <- read_csv("https://docs.google.com/spreadsheets/d/e/2PACX-1vSHQ17f-felEDBHEx2uz0Jjst-_z6wHTfB7hSFHD1iX8XWwd6GYTOYBo5F8DXjLGdxpTk9BqUt_M0lG/pub?gid=1995037555&single=true&output=csv")
glimpse(df) 

Data manipulation

Creo una nueva variable: status

Pronóstico pulpar: Variable dependiente. Evaluado con radiografÍa periapical. - Buen pronóstico: (PS) sobrevivencia de la Pulpa - Sin cambio radiográfico, signos clÍnicos de vitalidad pulpar. - PCO: obliteración del conducto pulpar. - Detención del crecimiento radicular, conducto pulpar normal, longitud radicular disminuida en comparación con el homólogo.

• Mal pronóstico: (PN) necrosis pulpar, observando amplitud del conducto radicular, espesor disminuido de paredes radiculares con o sin lesión apical - Desarrollo radicular completo, con tratamiento de endodoncia. - Detención del crecimiento radicular, con tratamiento de endodoncia.

Entonces

Sobrevida si - PCO - Detención crecimiento radicular, pulpa vital - Sin cambio rx

Sobrevida no - Desarrollo Radicular completo, con endodoncia - Detención crecimiento radicular, con endodoncia

Variables importantes

• tiempo
• evento
• independientes

Tiempo

cambio las fechas a aaaa/mm/dd y calculo la diferencia en años

df <- df %>% 
  mutate(
    `FECHA TRAUMATISMO` = dmy(`FECHA TRAUMATISMO`), 
    `FECHA SOBREVIDA` = dmy(`FECHA SOBREVIDA`), 
    `FECHA DE NACIMIENTO` = dmy(`FECHA DE NACIMIENTO`)
  )

# creo time y edad

df <- df %>% 
  mutate(
    time = difftime(`FECHA SOBREVIDA`, `FECHA TRAUMATISMO`, units = "weeks") / 52.25, 
    edad =  difftime(`FECHA TRAUMATISMO`, `FECHA DE NACIMIENTO`, units = "weeks")/52.25 # edad en años
  )

Eventos

Creo una nueva variable que será el evento = 1 (si) en caso que haya algo negativo y no (0) en caso que no haya pasado nada malo

df <-  df %>% 
  mutate(
    Evento = case_when(
      `SOBREVIDA PULPAR` == "PCO" |
        `SOBREVIDA PULPAR` == "Detención crecimiento radicular, pulpa vital" |
        `SOBREVIDA PULPAR` == "Sin cambio rx" ~ 0, 
      TRUE ~ 1

    )
  )

Veo cuantos eventos hay

table(df$Evento)

Dataset clean

CORREGIDO Hay dos valores negativos, mientras los elimino

which(df\(time > 0) df[df\)time < 0, ] # ubico a la que es menor df\(time df[df\)time < 1, ] df <- df[-c(65), ] # elimino la 65 mientras tanto VERIFICAR!!!!

which(df\(edad < 0) df[df\)edad < 1, ] # hay uno que tiene edad negativa, verificar, mientras lo saco del análisis df <- df[-c(77), ]

summary(df)

Survival analysis

Creo el objeto estándar de sobrevida

km <-  with(df, Surv(time, Evento == 1))
head(km, 50)

Análisis de Kaplan Meier

Ahora el análisis de KM

km_fit <- survfit(Surv(time, Evento) ~ 1, data = df)
summary(km_fit)

Estos datos se leen de la siguiente manera:

• en la columna time aparece el tiempo. Este es el eje X del gráfico;
• en n.risk el n de individuos en riesgo;
• en n.event la cantidad de eventos en ese perÍodo;
• en survival la función de sobrevida. Corresponde al eje Y del gráfico;
• luego el error std y los intervalos de confianza, que aparecen marcados en tono gris.

km_fit

autoplot(km_fit, 
         main = "Análisis de sobrevida Kaplan Meier para estado pulpar",
         xlab = "Años", ylab = "Probabilidad de sobrevida")

ggsurvplot(km_fit, main = "Survival curve",
   font.main = 18,
   font.x =  16,
   font.y = 16,
   font.tickslab = 14, 
   palette = "Dark Blue", 
   legend = "none", 
   ggtheme = theme_pubclean() ) +
        labs(y = "Probabilidad de sobrevida", x = "Años", 
             title = "Curva de Kaplan-Meier para sobrevida pulpar")

Ahora examino por algunas variables: sexo, dgo del trauma

autoplot(
  survfit(Surv(time, Evento) ~ GÉNERO, data = df), 
         main = "Análisis de sobrevida Kaplan Meier para estado pulpar por sexo",
         xlab = "Años", ylab = "Probabilidad de sobrevida"
) +
        theme(text = element_text(size=20))

autoplot(
  survfit(Surv(time, Evento) ~ DIAGNÓSTICO, data = df), 
          main = "Análisis de sobrevida Kaplan Meier para estado pulpar por tipo de trauma",
         xlab = "Años", ylab = "Probabilidad de sobrevida"
) +
        theme(text = element_text(size=20))

19 dic hacer solo para FCNC con y sin luxación

df_soloFCNC <- df %>% 
       filter(DIAGNÓSTICO %in% c("FCNC sin luxación", "FCNC con luxación")) 




autoplot(
  survfit(Surv(time, Evento) ~ DIAGNÓSTICO, data = df_soloFCNC), 
          main = "Análisis de sobrevida Kaplan Meier para estado pulpar por tipo de trauma",
         xlab = "Años", ylab = "Probabilidad de sobrevida"
) +
        theme(text = element_text(size=20))

Análisis de Cox

El propósito del modelo de Cox es evaluar simultáneamente el efecto de varios factores en la supervivencia. Nos permite examinar cómo diversos influyen en la tasa de un evento en particular que ocurre (por ejemplo, infección, muerte) en un momento particular. Esta tasa se conoce comúnmente como la tasa de riesgo (hazard rate). Las variables (o factores) del predictor generalmente se denominan covariables en la literatura de análisis de supervivencia.

h(t)=h0(t)×exp(b1x1+b2x2+…+bpxp) h(t)=h0(t)×exp(b1x1+b2x2+…+bpxp)

donde,

• t representa el tiempo de supervivencia
• h (t) h (t) es la función de riesgo determinada por un conjunto de covariables p (x1, x2, …, xpx1, x2, …, xp)
• los coeficientes (b1, b2, …, bpb1, b2, …, bp) miden el impacto (es decir, el tamaño del efecto) de las covariables. el término h0h0 se llama riesgo de referencia. Corresponde al valor del peligro si todos los xixi son iguales a cero (la cantidad exp (0) es igual a 1). La ‘t’ en h (t) nos recuerda que el peligro puede variar con el tiempo.

El modelo de Cox se puede escribir como una regresión lineal múltiple del logaritmo del peligro en las variables xixi, con el riesgo de la lÍnea de base como un término de “intercepción” que varÍa con el tiempo.

Las cantidades exp (bi) exp (bi) se denominan relaciones de riesgo (HR). Un valor de bibi mayor que cero, o equivalente a una razón de riesgo mayor que uno, indica que a medida que el valor de la novena covariable aumenta, el riesgo de evento aumenta y, por lo tanto, la duración de la supervivencia disminuye.

Dicho de otra manera, una razón de riesgo superior a 1 indica una covariable que está asociada positivamente con la probabilidad del evento y, por lo tanto, está asociada negativamente con la duración de la supervivencia.

En resumen,

HR = 1: sin efecto HR < 1: reducción en el riesgo HR > 1: aumento de riesgo

res.cox1 <- coxph(km ~ 
                    edad + 
                    `DIENTE AFECTADO` +
                    GÉNERO + 
                    DIAGNÓSTICO + 
                    `TRAUMATISMO REPETIDO` +
                    `TRAUMATISMO PREVIO`,
                  data =  df)

Chequeo por si se viola el riesgo proporcional (constante HR en el tiempo)

(res.zph1 <- cox.zph(res.cox1))

No se viola, asÍ que estamos OK

Ahora examino con detalle el modelo de Cox

summary(res.cox1)

Esto se lee de la siguiente manera

individuos que tienen más edad tienen más probabilidad (0.4841) de tener el evento, mientras que ser sexo masculino disminuye el riesgo (-0.0337) de tener el evento.

¿Cuanto?

UN aumento de Edad aumenta un 62% el riesgo de tener el evento, con un ic95% entre 05% y 149% (restale uno a los coeficientes)

cox_fit <- survfit(res.cox1)
autoplot(cox_fit)

aa_fit <-aareg(Surv(time, Evento) ~ 
                 `DIENTE AFECTADO` +
                    GÉNERO + 
                    DIAGNÓSTICO + 
                    `TRAUMATISMO REPETIDO` +
                    df$`TRAUMATISMO PREVIO`,
                  data =  df)

autoplot(aa_fit)

TODOS

modelo_0 <- coxph(km ~ 
                    df$GÉNERO + 
                    df$`DIENTE AFECTADO` +
                    df$DIAGNÓSTICO + 
                    df$TRATAMIENTO + 
                    # df$`TIEMPO HASTA LA ATENCIÓN (en horas)` + 
                    df$`ESTADO RADICULAR INICIAL` + 
                    df$`TRAUMATISMO REPETIDO` + 
                    df$`TRAUMATISMO PREVIO` + 
                    df$`INDICACIÓN DE HIGIENE` + 
                    df$`EVALUACIÓN DE HIGIENE` + 
                    df$`ONICOFAGIA - INTERPOSICIÓN DE OBJETOS` , 
                  data =  df)

summary(modelo_0)

TRAUMATISMO PREVIO

modelo 1

modelo 1.1 traumatismo repetido por si solo

modelo_1.1 <- coxph(km ~ 
                    # df$GÉNERO + 
                    # df$`DIENTE AFECTADO` +
                    # df$DIAGNÓSTICO + 
                    # df$TRATAMIENTO + 
                    # df$`TIEMPO HASTA LA ATENCIÓN (en horas)` + 
                    # df$`ESTADO RADICULAR INICIAL` + 
                     df$`TRAUMATISMO REPETIDO`,
                    # df$`TRAUMATISMO PREVIO` + 
                    # df$`INDICACIÓN DE HIGIENE` + 
                    # df$`EVALUACIÓN DE HIGIENE` + 
                    #df$`ONICOFAGIA - INTERPOSICIÓN DE OBJETOS` , 
                  data =  df)
summary(modelo_1.1)

modelo 1.2 traumatismo previo por si solo

modelo_1.2 <- coxph(km ~ 
                    # df$GÉNERO + 
                    # df$`DIENTE AFECTADO` +
                    # df$DIAGNÓSTICO + 
                    # df$TRATAMIENTO + 
                    # df$`TIEMPO HASTA LA ATENCIÓN (en horas)` + 
                    # df$`ESTADO RADICULAR INICIAL` + 
                    #  df$`TRAUMATISMO REPETIDO`,
                     df$`TRAUMATISMO PREVIO`,
                    # df$`INDICACIÓN DE HIGIENE` + 
                    # df$`EVALUACIÓN DE HIGIENE` + 
                    #df$`ONICOFAGIA - INTERPOSICIÓN DE OBJETOS` , 
                  data =  df)
summary(modelo_1.2)

EVALUACIÓN DE HIGIENE

modelo 2: indicación y evaluación de higiene

modelo_2 <- coxph(km ~ 
                    # df$GÉNERO + 
                    # df$`DIENTE AFECTADO` +
                    # df$DIAGNÓSTICO + 
                    # df$TRATAMIENTO + 
                    # df$`TIEMPO HASTA LA ATENCIÓN (en horas)` + 
                    # df$`ESTADO RADICULAR INICIAL` + 
                    # df$`TRAUMATISMO REPETIDO`,
                    # df$`TRAUMATISMO PREVIO` + 
                     df$`INDICACIÓN DE HIGIENE` + 
                     df$`EVALUACIÓN DE HIGIENE`,
                    #df$`ONICOFAGIA - INTERPOSICIÓN DE OBJETOS` , 
                  data =  df)
summary(modelo_2)

modelo 2.1 indicación de higiene por si solo

modelo_2.1 <- coxph(km ~ 
                    # df$GÉNERO + 
                    # df$`DIENTE AFECTADO` +
                    # df$DIAGNÓSTICO + 
                    # df$TRATAMIENTO + 
                    # df$`TIEMPO HASTA LA ATENCIÓN (en horas)` + 
                    # df$`ESTADO RADICULAR INICIAL` + 
                    # df$`TRAUMATISMO REPETIDO`,
                    # df$`TRAUMATISMO PREVIO` + 
                     df$`INDICACIÓN DE HIGIENE`, 
                    # df$`EVALUACIÓN DE HIGIENE` + 
                    #df$`ONICOFAGIA - INTERPOSICIÓN DE OBJETOS` , 
                  data =  df)
summary(modelo_2.1)

modelo 2.2 evaluación de higiene por si solo

modelo_2.2 <- coxph(km ~ 
                    # df$GÉNERO + 
                    # df$`DIENTE AFECTADO` +
                    # df$DIAGNÓSTICO + 
                    # df$TRATAMIENTO + 
                    # df$`TIEMPO HASTA LA ATENCIÓN (en horas)` + 
                    # df$`ESTADO RADICULAR INICIAL` + 
                    # df$`TRAUMATISMO REPETIDO`,
                    # df$`TRAUMATISMO PREVIO` + 
                    # df$`INDICACIÓN DE HIGIENE` + 
                     df$`EVALUACIÓN DE HIGIENE`,
                    #df$`ONICOFAGIA - INTERPOSICIÓN DE OBJETOS` , 
                  data =  df)
summary(modelo_2.2)

ONICOFAGIA - INTERPOSICIÓN DE OBJETOS

modelo 3: malos hábitos

modelo_3 <- coxph(km ~ 
                    # df$GÉNERO + 
                    # df$`DIENTE AFECTADO` +
                    # df$DIAGNÓSTICO + 
                    # df$TRATAMIENTO + 
                    # df$`TIEMPO HASTA LA ATENCIÓN (en horas)` + 
                    # df$`ESTADO RADICULAR INICIAL` + 
                    # df$`TRAUMATISMO REPETIDO`,
                    # df$`TRAUMATISMO PREVIO` + 
                    # df$`INDICACIÓN DE HIGIENE` + 
                    # df$`EVALUACIÓN DE HIGIENE` + 
                    df$`ONICOFAGIA - INTERPOSICIÓN DE OBJETOS` , 
                  data =  df)
summary(modelo_3)

ESTADO RADICULAR INICIAL

modelo 4: estado radicular moorres

modelo_4 <- coxph(km ~ 
                    # df$GÉNERO + 
                    # df$`DIENTE AFECTADO` +
                    # df$DIAGNÓSTICO + 
                    # df$TRATAMIENTO + 
                    # df$`TIEMPO HASTA LA ATENCIÓN (en horas)` + 
                     df$`ESTADO RADICULAR INICIAL`,
                    # df$`TRAUMATISMO REPETIDO`,
                    # df$`TRAUMATISMO PREVIO` + 
                    # df$`INDICACIÓN DE HIGIENE` + 
                    # df$`EVALUACIÓN DE HIGIENE` + 
                    #df$`ONICOFAGIA - INTERPOSICIÓN DE OBJETOS` , 
                  data =  df)
summary(modelo_4)

Modelo 99 que sugiero

modelo_99 <- coxph(km ~ 
                   GÉNERO + 
                    # df$`DIENTE AFECTADO` +
                   DIAGNÓSTICO + 
                    # df$TRATAMIENTO + 
                    # df$`TIEMPO HASTA LA ATENCIÓN (en horas)` + 
                   `ESTADO RADICULAR INICIAL` + 
                    `TRAUMATISMO REPETIDO` +
                   `TRAUMATISMO PREVIO` + 
                    # df$`INDICACIÓN DE HIGIENE` + 
                    `EVALUACIÓN DE HIGIENE` + 
                    `ONICOFAGIA - INTERPOSICIÓN DE OBJETOS` , 
                  data =  df)
summary(modelo_99)

Chequeo por si se viola el riesgo proporcional (constante HR en el tiempo)

(res.zph99 <- cox.zph(modelo_99))

cox_fit_99 <- survfit(modelo_99)
autoplot(cox_fit_99)

Ojo, la mediana no toca el 50%, por lo que da median = NA

cox_fit_99

summary(cox_fit_99)

aa_fit_99 <-aareg(Surv(time, Evento) ~ 
                  GÉNERO + 
                    # df$`DIENTE AFECTADO` +
                     DIAGNÓSTICO + 
                    # df$TRATAMIENTO + 
                    # df$`TIEMPO HASTA LA ATENCIÓN (en horas)` + 
                     `ESTADO RADICULAR INICIAL` + 
                     `TRAUMATISMO REPETIDO` +
                     `TRAUMATISMO PREVIO` + 
                    # df$`INDICACIÓN DE HIGIENE` + 
                     `EVALUACIÓN DE HIGIENE` + 
                     `ONICOFAGIA - INTERPOSICIÓN DE OBJETOS` , 
                  data =  df)

autoplot(aa_fit_99)

Tamaño muestral y potencia

hr = 2 # hazard ratio
hr0 = 1
pE = 0.8 # probabilidad del evento
pA = 0.5 # proporción en el grupo A
alpha = 0.05
beta = 0.20
(n=((qnorm(1-alpha/2)+qnorm(1-beta))/(log(hr)-log(hr0)))^2/(pA*(1-pA)*pE))
ceiling(n) # 82
(Power=pnorm((log(hr)-log(hr0))*sqrt(n*pA*(1-pA)*pE)-qnorm(1-alpha/2)))

Referencias

https://analyticsbuddhu.wordpress.com/2017/04/17/how-to-do-coxproportional-hazard-regression-modelling-using-r/

https://www.openintro.org/download.php?file=survival_analysis_in_R&referrer=/stat/surv.php

https://courses.nus.edu.sg/course/stacar/internet/st3242/handouts/notes3.pdf

Chow S, Shao J, Wang H. 2008. Sample Size Calculations in Clinical Research. 2nd Ed. Chapman & Hall/CRC Biostatistics Series. page 177.