Primero se carga los paquetes necesarios.

library(tidyverse)
library(broom)
library(lubridate)
library(scales)
library(skimr)
library(janitor)
library(ggsci)
library(RColorBrewer)
library(arules)

library(survival)
library(survminer)
library(KMsurv)

Se establece el tema y los colores.

theme_set(theme_bw() + theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5)))
pastel1 <- brewer.pal(9, "Pastel1")
pastel2 <- brewer.pal(8, "Pastel2")

color_prin <- pastel1[3]
color_sec <- pastel1[2]

# display.brewer.pal(9, 'Pastel1') display.brewer.pal(8, 'Pastel2')

1 Análisis descriptivo

Se carga los datos:

data("larynx", package = "KMsurv")
help("larynx", package = "KMsurv")

Se trata de información acerca de 90 pacientes de con cáncer de laringe, las variables del conjunto de datos son:

  • stage: es la etapa de la enfermedad del paciente.

  • time: es el tiempo en meses transcurrido dentro del estudio, puede que el paciente falleciera o que se censurara.

  • age: es la edad del paciente al momento de recibir el diagnóstico de cáncer.

  • diagyr: es el año en que se diagnóstico la enfermedad.

  • delta: indica si el paciente murió (1) o fue censurado (0).

El conjunto de datos es el siguiente:

larynx_fct <- larynx %>%
    mutate(stage = as_factor(stage), stage = fct_relevel(stage, c("1", "2", "3",
        "4")), stage = fct_recode(stage, Etapa_1 = "1", Etapa_2 = "2", Etapa_3 = "3",
        Etapa_4 = "4"), stage = as_factor(str_replace(stage, pattern = "_", replacement = " ")),
        delta = as_factor(delta), delta = fct_recode(delta, Muerto = "1", Vivo = "0"))

larynx_fct

y en seguida se muestra un breve análisis descriptivo:

skim(larynx_fct)
Data summary
Name larynx_fct
Number of rows 90
Number of columns 5
_______________________
Column type frequency:
factor 2
numeric 3
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
stage 0 1 FALSE 4 Eta: 33, Eta: 27, Eta: 17, Eta: 13
delta 0 1 FALSE 2 Mue: 50, Viv: 40

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
time 0 1 4.20 2.63 0.1 2.00 4 6.2 10.7 ▇▇▆▅▂
age 0 1 64.61 10.80 41.0 57.00 65 72.0 86.0 ▃▅▇▆▃
diagyr 0 1 74.24 2.19 70.0 72.25 74 76.0 78.0 ▅▆▅▇▅

Se puede observar que:

  • La persona más joven tuvo 41 años y la más grande 86 años.

  • El primer diagnóstico se realizó en 1970 mientras que el último se llevó a cabo en 1978.

  • La primer persona en salir del estudio lo hizo en alrededor de 3 días, mientras que la última tardó 10 meses con 21 días.

Es interesante saber cuál es la etapa de la enfermedad más común:

larynx %>%
    group_by(stage) %>%
    summarise(total = n()) %>%
    ggplot(aes(stage, total, label = total)) + geom_col(fill = color_prin) + geom_label(position = position_dodge(1),
    hjust = -0.1) + scale_x_continuous(breaks = scales::breaks_width(width = 1)) +
    scale_y_continuous(breaks = scales::breaks_width(width = 2)) + coord_flip() +
    labs(title = "Cantidad de pacientes por etapa de la enfermedad", x = "Etapa de la enfermedad",
        y = "Número de observaciones")

La etapa 1 es la más frecuente y es más del doble de común que la etapa 4, la cual afortunadamente es la menos reportada. Además la etapa 3 es la segunda más frecuente, lo que indica que si bien muchos pacientes no avanzan en su enfermedad, hay bastantes que se acercan a la etapa más avanzada y esto puede tener un impacto en la cantidad de fallecidos.

Para revisar esto, se muestra la cantidad de censurados y muertos:

larynx_fct %>%
    group_by(delta) %>%
    summarise(total = n()) %>%
    ggplot(aes(delta, total, label = total)) + geom_col(fill = color_prin) + geom_label(position = position_dodge(1),
    hjust = -0.1) + scale_y_continuous(breaks = scales::breaks_width(width = 5)) +
    coord_flip() + labs(title = "Cantidad de pacientes por estado de salud", x = NULL,
    y = "Número de observaciones")

y se confirma que la mayoría de los pacientes fallecieron, sin embargo la diferencia no es considerable debido quizá a que la mayor parte de las personas solo tuvieron la primer etapa del cáncer.

También es importante conocer en qué edades hay más diagnósticos de la enfermedad:

intervalos <- c("41 a 49", "50 a 58", "59 a 67", "68 a 76", "77 a 87")

larynx %>%
    mutate(grupo_edad = discretize(age, method = "interval", breaks = 5, labels = intervalos)) %>%
    group_by(grupo_edad) %>%
    summarise(total = n()) %>%
    ggplot(aes(grupo_edad, total, label = total)) + geom_col(fill = color_prin) +
    geom_label(position = position_dodge(1), hjust = -0.1) + scale_y_continuous(breaks = scales::breaks_width(width = 2)) +
    coord_flip() + labs(title = "Edad de diagnóstico de cáncer", x = "Edad (años)",
    y = "Número de observaciones")

y es evidente que la mayoría de los pacientes tenían de 59 a 76 años cuando recibieron el diagnóstico; además se nota que son pocas las personas más jóvenes y más grandes (su suma es apenas igual al número de pacientes de 58 a 67 años).

Para conocer el año en que más diagnósticos hubo:

larynx_fct %>%
    group_by(diagyr) %>%
    summarise(total = n()) %>%
    ggplot(aes(diagyr, total, label = total)) + geom_col(fill = color_prin) + geom_label(position = position_dodge(1),
    hjust = -0.1) + scale_y_continuous(breaks = scales::breaks_width(width = 2)) +
    coord_flip() + labs(title = "Año de diagnóstico de cáncer", x = "Año", y = "Número de observaciones")

y se ve que en 1976 hubo más diagnósticos, además en 1970 solo hubo dos mientras que en 1978 únicamente 4.

Ahora se analiza la distribución del tiempo en el estudio. Primero se considera el estado de salud:

larynx_fct %>%
    ggplot(aes(time, fill = delta)) + geom_histogram(bins = 15) + scale_x_continuous(breaks = scales::breaks_width(width = 1)) +
    scale_y_continuous(breaks = scales::breaks_width(width = 1)) + scale_fill_manual(values = pastel1[c(3,
    2)]) + coord_flip() + labs(title = "Distribución del tiempo durante el estudio",
    fill = "Estado", x = "Tiempo (meses)", y = "Número de observaciones")

resalta que hay pacientes que no sobrevivieron ni un mes, que los que permanecieron vivos estuvieron al menos dos meses en el estudio y que los fallecidos no pasaron de los 8 meses. La mayoría de los pacientes estuvo alrededor de 4 meses pero fueron pocos los que fallecieron después de este tiempo.

Por otro lado, considerando la distribución del tiempo de observación respecto a las etapas del cáncer:

larynx_fct %>%
    ggplot(aes(time, fill = stage)) + geom_histogram(bins = 15) + scale_x_continuous(breaks = scales::breaks_width(width = 1)) +
    scale_y_continuous(breaks = scales::breaks_width(width = 2)) + scale_fill_manual(values = pastel1) +
    facet_wrap(~stage) + labs(title = "Distribución del tiempo durante el estudio por nivel de enfermdad",
    fill = "Nivel de la enfermedad", x = "Tiempo (meses)", y = "Número de observaciones")

Se puede decir que las personas que permanecieron menos tiempo tenían una etapa avanzada de la enfermedad y que las personas que se encontraban en una etapa temprana por lo regular permanecieron más de dos meses pero menos de 9.

Para terminar el análisis descriptivo se muestra la relación entre la edad de los pacientes y el tiempo que permanecieron en el estudio:

larynx_fct %>%
    ggplot(aes(age, time, colour = stage)) + geom_point() + scale_x_continuous(breaks = scales::breaks_width(width = 5)) +
    scale_y_continuous(breaks = scales::breaks_width(width = 2)) + scale_colour_manual(values = pastel1) +
    labs(title = "Relación entre la edad y el tiempo de estudio", colour = "Etapa de la enfermedad",
        x = "Edad (años)", y = "Tiempo (meses)")

No hay una relación clara entre ambas varales, por lo general las personas de edad cercana a la mediana padecieron la etapa 4 y los jóvenes solían solo llegar a la primer etapa. Además, quienes tenían una etapa avanzada de la enfermedad permanecieron menos tiempo en el estudio.

2 Estimación Kaplan-Meier

Todos los intervalos de confianza son al nivel 95% y se calcularon mediante el método “plain” que es el que se revisó en clase.

2.1 Población completa

Se ajusta el modelo Kaplan-Meier:

km <- survfit(Surv(time, delta) ~ 1, type = "kaplan-meier", data = larynx, conf.type = "plain")

La función de sobrevivencia \(\hat{S}\) junto con la mediana del tiempo de vida es:

km_g <- ggsurvplot(fit = km, data = larynx, conf.int = TRUE, surv.median.line = "hv",
    break.time.by = 1, legend = "none", risk.table = TRUE, cumevents = TRUE, cumcensor = TRUE,
    cumcensor.title = "Censuras acumuladas", tables.height = 0.21, ggtheme = theme_bw(),
    palette = color_prin)

km_g$plot <- km_g$plot + labs(title = "Estimación Kaplan-Meier", x = "Tiempo (meses)") +
    theme(axis.title.y = element_blank())

km_g$table <- km_g$table + labs(title = "Pacientes en riesgo", x = "Tiempo (meses)",
    y = "Probabilidad de superviencia")

km_g$cumevents <- km_g$cumevents + labs(title = "Eventos acumulados", x = "Tiempo (meses)") +
    theme(axis.title.y = element_blank())

km_g$ncensor.plot <- km_g$ncensor.plot + labs(title = "Censuras acumuladas", x = "Tiempo (meses)") +
    theme(axis.title.y = element_blank())

km_g$plot

También resulta ilustrativo mostrar la gráfica junto con los expuestos al riesgo de muerte, el número de decesos y las censuras acumuladas a lo largo del tiempo:

print(km_g, surv.plot.height = 2)

2.2 Etapa de la enfermedad

larynx_fct_etapa <- larynx %>%
    mutate(stage = as_factor(stage), stage = fct_relevel(stage, c("1", "2", "3",
        "4")), stage = fct_recode(stage, Etapa_1 = "1", Etapa_2 = "2", Etapa_3 = "3",
        Etapa_4 = "4"), stage = as_factor(str_replace(stage, pattern = "_", replacement = " ")))

km_stage <- survfit(Surv(time, delta) ~ stage, type = "kaplan-meier", data = larynx_fct_etapa,
    conf.type = "plain")

km_stage_g <- ggsurvplot(fit = km_stage, data = larynx_fct_etapa, conf.int = TRUE,
    surv.median.line = "hv", break.time.by = 1, risk.table = TRUE, cumevents = TRUE,
    cumcensor = TRUE, cumcensor.title = "Censuras acumuladas", tables.height = 0.21,
    ggtheme = theme_bw(), palette = pastel1, title = "Estimación Kaplan-Meir por etapa de enfermedad",
    legend.title = "Etapa cáncer")

km_stage_g$plot <- km_stage_g$plot + labs(title = "Estimación Kaplan-Meier", x = "Tiempo (meses)",
    y = "Probabilidad de supervivencia")
km_stage_g$plot

Es claro que las etapas 1 y 2 tienen un gran parecido y que poseen una mayor probabilidad de supervivencia que las otras etapas, así que se combina los niveles 1 y 2 de la enfermedad en la categoría “temprana”, mientras que los niveles 3 y 4 se unen en la categoría “avanzada”:

larynx_fct_etapa_colapso <- larynx_fct_etapa %>%
    mutate(stage = fct_collapse(stage, etapa_temprana = c("Etapa 1", "Etapa 2"),
        etapa_avanzada = c("Etapa 3", "Etapa 4")))

km_stage_colapso <- survfit(Surv(time, delta) ~ stage, type = "kaplan-meier", data = larynx_fct_etapa_colapso,
    conf.type = "plain")

km_stage_colapso_g <- ggsurvplot(fit = km_stage_colapso, data = larynx_fct_etapa_colapso,
    conf.int = TRUE, surv.median.line = "hv", break.time.by = 1, risk.table = TRUE,
    cumevents = TRUE, cumcensor = TRUE, cumcensor.title = "Censuras acumuladas",
    tables.height = 0.21, ggtheme = theme_bw(), palette = pastel1, legend.title = "Etapa",
    legend.labs = c("Temprana", "Avanzada"))


km_stage_colapso_g$plot <- km_stage_colapso_g$plot + labs(title = "Estimación Kaplan-Meier",
    x = "Tiempo (meses)", y = "Probabilidad de supervivencia")

km_stage_colapso_g$plot

y se obtiene una diferencia sustancial, las dos curvas no se cortan y sus intervalos de confianza apenas y se sobreponen, por lo que se concluye que la etapa del cáncer es una variable que afecta al tiempo de supervivencia, sin embargo habrá que realizar las pruebas de hipótesis apropiadas para comprobarlo con certeza.

2.3 Grupos de edad

mediana_edad <- median(larynx$age)

La edad mediana es 65, así que se separa a los pacientes en base a si son menores o mayores a esta edad:

larynx_edad_mediana <- larynx %>%
    mutate(grupo_edad = as_factor(ifelse(age <= 65, "Menor o igual a 65", "Mayor a 65")))

km_edad <- survfit(Surv(time, delta) ~ grupo_edad, type = "kaplan-meier", data = larynx_edad_mediana,
    conf.type = "plain")

km_edad_g <- ggsurvplot(fit = km_edad, data = larynx_edad_mediana, conf.int = TRUE,
    surv.median.line = "hv", break.time.by = 1, risk.table = TRUE, cumevents = TRUE,
    cumcensor = TRUE, cumcensor.title = "Censuras acumuladas", tables.height = 0.21,
    ggtheme = theme_bw(), palette = pastel1, legend.title = "Grupo de edad", legend.labs = c("Mayor a 65",
        "Menor a 65"))


km_edad_g$plot <- km_edad_g$plot + labs(title = "Estimación Kaplan-Meier", x = "Tiempo (meses)",
    y = "Probabilidad de supervivencia")

km_edad_g$plot

pero las dos curvas no parecen tener una diferencia sustancial. Se probará separando a la población en 3 grupos

larynx_edad2 <- larynx %>%
    mutate(grupo_edad = as_factor(ifelse(age < 60, "Menor a 60", ifelse(age < 70,
        "De 60 a 70", "Mayor a 70"))))

km_edad2 <- survfit(Surv(time, delta) ~ grupo_edad, type = "kaplan-meier", data = larynx_edad2,
    conf.type = "plain")

km_edad2_g <- ggsurvplot(fit = km_edad2, data = larynx_edad2, conf.int = TRUE, surv.median.line = "hv",
    break.time.by = 1, risk.table = TRUE, cumevents = TRUE, cumcensor = TRUE, cumcensor.title = "Censuras acumuladas",
    tables.height = 0.21, ggtheme = theme_bw(), palette = pastel1, legend.title = "Grupo de edad")


km_edad2_g <- km_edad2_g + labs(title = "Estimación Kaplan-Meier", x = "Tiempo (meses)",
    y = "Probabilidad de supervivencia")

km_edad2_g$plot

Pero en este caso tampoco resulta sustancial la separación, por lo que se concluye que la edad no es una variable que afecte el tiempo de supervivencia, sin embargo habrá que realizar las pruebas de hipótesis apropiadas para comprobarlo con certeza.

3 Pruebas de hipótesis

3.1 Etapa de la enfermedad

Realizando las pruebas en base a la etapa del cáncer:

survdiff(Surv(time, delta) ~ stage, rho = 0, data = larynx_fct_etapa)
## Call:
## survdiff(formula = Surv(time, delta) ~ stage, data = larynx_fct_etapa, 
##     rho = 0)
## 
##                N Observed Expected (O-E)^2/E (O-E)^2/V
## stage=Etapa 1 33       15    22.57     2.537     4.741
## stage=Etapa 2 17        7    10.01     0.906     1.152
## stage=Etapa 3 27       17    14.08     0.603     0.856
## stage=Etapa 4 13       11     3.34    17.590    19.827
## 
##  Chisq= 22.8  on 3 degrees of freedom, p= 5e-05
survdiff(Surv(time, delta) ~ stage, rho = 1, data = larynx_fct_etapa)
## Call:
## survdiff(formula = Surv(time, delta) ~ stage, data = larynx_fct_etapa, 
##     rho = 1)
## 
##                N Observed Expected (O-E)^2/E (O-E)^2/V
## stage=Etapa 1 33     9.26    15.78     2.694      6.48
## stage=Etapa 2 17     4.66     7.18     0.884      1.47
## stage=Etapa 3 27    12.77    10.06     0.730      1.35
## stage=Etapa 4 13     9.15     2.82    14.219     18.78
## 
##  Chisq= 23.1  on 3 degrees of freedom, p= 4e-05

Tanto la prueba de Mantel-Haenszel (\(\rho = 0\)) como la de Peto-Peto (\(\rho = 1\)) rechazan la hipótesis nula de igualdad en las funciones de supervivencia (con p-values considerablemente menores al nivel de significancia de 5%), por lo que se reafirma la aseveración hecha al analizar las gráficas de que la variable de etapa de la enfermedad tiene una influencia considerable en la supervivencia de la población.

Cuando se colapsa los niveles 1 y 2 en una categoría y los niveles 3 y 4 en otra, se obtiene un resultado similar:

survdiff(Surv(time, delta) ~ stage, rho = 0, data = larynx_fct_etapa_colapso)
## Call:
## survdiff(formula = Surv(time, delta) ~ stage, data = larynx_fct_etapa_colapso, 
##     rho = 0)
## 
##                       N Observed Expected (O-E)^2/E (O-E)^2/V
## stage=etapa_temprana 50       22     32.6      3.43      10.1
## stage=etapa_avanzada 40       28     17.4      6.42      10.1
## 
##  Chisq= 10.1  on 1 degrees of freedom, p= 0.001
survdiff(Surv(time, delta) ~ stage, rho = 1, data = larynx_fct_etapa_colapso)
## Call:
## survdiff(formula = Surv(time, delta) ~ stage, data = larynx_fct_etapa_colapso, 
##     rho = 1)
## 
##                       N Observed Expected (O-E)^2/E (O-E)^2/V
## stage=etapa_temprana 50     13.9     23.0      3.56      13.1
## stage=etapa_avanzada 40     21.9     12.9      6.35      13.1
## 
##  Chisq= 13.1  on 1 degrees of freedom, p= 3e-04

3.2 Grupos de edad

Primero se considera la separación en dos grupos de acuerdo a la mediana de la edad:

survdiff(Surv(time, delta) ~ grupo_edad, rho = 0, data = larynx_edad_mediana)
## Call:
## survdiff(formula = Surv(time, delta) ~ grupo_edad, data = larynx_edad_mediana, 
##     rho = 0)
## 
##                                N Observed Expected (O-E)^2/E (O-E)^2/V
## grupo_edad=Mayor a 65         44       27     25.1     0.146     0.298
## grupo_edad=Menor o igual a 65 46       23     24.9     0.147     0.298
## 
##  Chisq= 0.3  on 1 degrees of freedom, p= 0.6
survdiff(Surv(time, delta) ~ grupo_edad, rho = 1, data = larynx_edad_mediana)
## Call:
## survdiff(formula = Surv(time, delta) ~ grupo_edad, data = larynx_edad_mediana, 
##     rho = 1)
## 
##                                N Observed Expected (O-E)^2/E (O-E)^2/V
## grupo_edad=Mayor a 65         44     18.7     17.8    0.0436     0.115
## grupo_edad=Menor o igual a 65 46     17.2     18.1    0.0429     0.115
## 
##  Chisq= 0.1  on 1 degrees of freedom, p= 0.7

Tanto la prueba de Mantel-Haenszel (\(\rho = 0\)) como la de Peto-Peto (\(\rho = 1\)) no rechazan la hipótesis nula de igualdad en las funciones de supervivencia (con p-values 6% y 7%, respectivamente), por lo que se reafirma la aseveración hecha al analizar las gráficas de que la variable de grupo de edad no tiene influencia en la supervivencia de la población.

Cuando se considera la otra clasificación de grupo de dad en tres grupos se obtiene resultados similares.

survdiff(Surv(time, delta) ~ grupo_edad, rho = 0, data = larynx_edad2)
## Call:
## survdiff(formula = Surv(time, delta) ~ grupo_edad, data = larynx_edad2, 
##     rho = 0)
## 
##                        N Observed Expected (O-E)^2/E (O-E)^2/V
## grupo_edad=Mayor a 70 31       22     15.3     2.978      4.37
## grupo_edad=Menor a 60 28       12     15.2     0.687      1.01
## grupo_edad=De 60 a 70 31       16     19.5     0.630      1.06
## 
##  Chisq= 4.4  on 2 degrees of freedom, p= 0.1
survdiff(Surv(time, delta) ~ grupo_edad, rho = 1, data = larynx_edad2)
## Call:
## survdiff(formula = Surv(time, delta) ~ grupo_edad, data = larynx_edad2, 
##     rho = 1)
## 
##                        N Observed Expected (O-E)^2/E (O-E)^2/V
## grupo_edad=Mayor a 70 31    15.97     11.2     2.077      3.99
## grupo_edad=Menor a 60 28     9.23     11.1     0.303      0.58
## grupo_edad=De 60 a 70 31    10.64     13.6     0.654      1.42
## 
##  Chisq= 4  on 2 degrees of freedom, p= 0.1