Supervivencia y mortalidad COVID-19
MDatos
2022-10-12
#Analisis de supervivencia en DATOS ABIERTOS COVID-19 MÉXICO ##Esta base ##de datos se descargó el día 27/09/2022
Tenemos la opción del paquete readr, pero es mejor que
lo abramos con el paquete arrow. Sin no lo tienes instalalo
con la función istall.packages() Par atrabajar con
arrow, necesitamos la función read_csv_arrow y
ubicamos la dirección de nuestro archivo csv (no olvidemos que deben
usarse los “/”, para navegar entre carpetas)
library(arrow)## Warning: package 'arrow' was built under R version 4.2.3##
## Attaching package: 'arrow'## The following object is masked from 'package:utils':
##
## timestampX220927COVID19MEXICO <- read_csv_arrow("C:/Users/fidel/OneDrive - CINVESTAV/datos abiertos COVID19 México Mayo 2021/analisis R/220927COVID19MEXICO.csv", skip_empty_rows = FALSE, na = c("NA"))
#library(readr)
#X220927COVID19MEXICO <- read_csv("C:/Users/fidel/OneDrive - CINVESTAV/datos abiertos COVID19 México Mayo 2021/analisis R/220927COVID19MEXICO.csv")
covid19mx<-X220927COVID19MEXICO #R es un lenguaje basado en objetos, entonces lo que estamos creando es un objeto nuevo con el nombre del archivo (sin extención) que hemos descargado de la pagina, usando el operador de asignación: "<-" menorque menos.
#str(covid19mx) #str() function in R Language is used for compactly displaying the internal structure of a R object.
head(covid19mx)#con esta función :head nos muestra los primeros 6 datos## FECHA_ACTUALIZACION ID_REGISTRO ORIGEN SECTOR ENTIDAD_UM SEXO ENTIDAD_NAC
## 1 2022-09-27 180725 2 9 9 2 9
## 2 2022-09-27 06fce8 1 12 7 1 7
## 3 2022-09-27 1a4a8d 1 12 23 2 27
## 4 2022-09-27 1e6dad 1 12 7 1 6
## 5 2022-09-27 06e951 1 12 24 2 24
## 6 2022-09-27 0f43ef 1 12 22 1 9
## ENTIDAD_RES MUNICIPIO_RES TIPO_PACIENTE FECHA_INGRESO FECHA_SINTOMAS
## 1 9 12 2 2022-01-19 2022-01-17
## 2 7 59 1 2022-03-05 2022-03-05
## 3 23 8 1 2022-07-22 2022-07-20
## 4 7 12 1 2022-08-23 2022-08-21
## 5 24 28 1 2022-07-07 2022-07-01
## 6 22 14 1 2022-01-20 2022-01-16
## FECHA_DEF INTUBADO NEUMONIA EDAD NACIONALIDAD EMBARAZO HABLA_LENGUA_INDIG
## 1 9999-99-99 2 2 33 1 97 2
## 2 9999-99-99 97 2 35 1 98 2
## 3 9999-99-99 97 2 49 1 97 2
## 4 9999-99-99 97 2 34 1 2 2
## 5 9999-99-99 97 2 40 1 97 2
## 6 9999-99-99 97 2 56 1 2 2
## INDIGENA DIABETES EPOC ASMA INMUSUPR HIPERTENSION OTRA_COM CARDIOVASCULAR
## 1 2 2 2 2 2 2 2 2
## 2 2 2 2 2 2 2 2 2
## 3 2 2 2 2 2 2 2 2
## 4 2 2 2 2 2 2 2 2
## 5 2 2 2 2 2 1 2 2
## 6 2 2 2 2 2 2 2 2
## OBESIDAD RENAL_CRONICA TABAQUISMO OTRO_CASO TOMA_MUESTRA_LAB RESULTADO_LAB
## 1 2 2 2 2 2 97
## 2 2 2 2 2 2 97
## 3 2 2 2 2 2 97
## 4 2 2 2 2 2 97
## 5 2 2 2 2 2 97
## 6 2 2 2 1 1 2
## TOMA_MUESTRA_ANTIGENO RESULTADO_ANTIGENO CLASIFICACION_FINAL MIGRANTE
## 1 1 2 7 99
## 2 1 2 7 99
## 3 1 2 7 99
## 4 2 97 6 99
## 5 1 2 7 99
## 6 1 2 7 99
## PAIS_NACIONALIDAD PAIS_ORIGEN UCI
## 1 México 97 2
## 2 México 97 97
## 3 México 97 97
## 4 México 97 97
## 5 México 97 97
## 6 México 97 97Para fines demostrativos de este tutorial, vamos a hacer un muestreo aleatorio. Pues el procesamiento de esta cantidad de datos podría tardar algo de tiempo, si tu quieres analizar el total de estos datos, te sugiero omitir este paso.
##Muestreo aleatorio simple
df<- sample(1:nrow(covid19mx), 500000)
dbcovid19mx<-covid19mx[df,]#cargaremos los paquetes necesarios con los que vamos atrabajar
#Paquetes:
library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.2.3## Warning: package 'ggplot2' was built under R version 4.2.3## Warning: package 'tibble' was built under R version 4.2.3## Warning: package 'tidyr' was built under R version 4.2.3## Warning: package 'readr' was built under R version 4.2.3## Warning: package 'purrr' was built under R version 4.2.3## Warning: package 'dplyr' was built under R version 4.2.3## Warning: package 'forcats' was built under R version 4.2.3## Warning: package 'lubridate' was built under R version 4.2.3## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.1 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ lubridate::duration() masks arrow::duration()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errorslibrary(lubridate)
library(gtsummary)## Warning: package 'gtsummary' was built under R version 4.2.3library(dlookr)## Warning in !is.null(rmarkdown::metadata$output) && rmarkdown::metadata$output
## %in% : 'length(x) = 3 > 1' in coercion to 'logical(1)'##
## Attaching package: 'dlookr'
##
## The following object is masked from 'package:tidyr':
##
## extract
##
## The following object is masked from 'package:base':
##
## transformlibrary(gtable)## Warning: package 'gtable' was built under R version 4.2.3library(gt)
library(ggpubr)## Warning: package 'ggpubr' was built under R version 4.2.3library(janitor)## Warning: package 'janitor' was built under R version 4.2.3##
## Attaching package: 'janitor'
##
## The following objects are masked from 'package:stats':
##
## chisq.test, fisher.test##ahora haremos una recodificación de nuestras variables que nos interesan, para esto necesitamos el diccionario que descargamos y que acompaña a nuestra base de datos
##de nuevo, exploraremos nuestras columnas y base de datos
#vamos a seleccionar unas variables de aqui: como SEXO,EDAD,TIPO_PACIENTE, NEUMONIA, DIABETES,OBESIDAD,
#CLASIFICACION_FINAL, con lo que vamos a trabajar como parte del analisis exploratorio
#str(dbcovid19mx)##vamos a renombrar las columnas en nuestra base de datos, para esto usamos la función “rename”. ###creamos una base de datos,con ayuda de un operador, denominado pipe o pipa que se escribe con la tecla control+shift+M “%>%”(estamos trabajando aqui con el paquete tydiverse) ##haremos una sobre escritura sobre el objeto con el mismo nombre que hemos puesto: “dbcovid19mx”, usando el comando
dbcovid19mx <- dbcovid19mx %>% rename("Genero" = "SEXO", "Edad" = "EDAD", "Paciente" = "TIPO_PACIENTE",
"Neumonia" = "NEUMONIA", "Diabetes" = "DIABETES", "Obesidad" ="OBESIDAD")
dbcovid19mx %>% select(Genero, Edad, Paciente, Neumonia, Diabetes, Obesidad)## # A tibble: 500,000 × 6
## Genero Edad Paciente Neumonia Diabetes Obesidad
## <int> <int> <int> <int> <int> <int>
## 1 2 47 1 2 2 2
## 2 2 28 1 2 2 2
## 3 1 39 1 2 2 2
## 4 2 8 1 2 2 2
## 5 1 60 1 2 2 2
## 6 1 42 1 2 2 2
## 7 2 44 1 2 2 2
## 8 2 84 2 1 2 2
## 9 1 27 1 2 2 2
## 10 2 22 1 2 2 2
## # ℹ 499,990 more rows##haremos una recodificación de las variables seleccionadas
#recodificar variables de interes: sobreescribimos sobre el mismo objeto usando operador pipe y con la función mutate() y recode(), presentes en el paquete dplyr de tydiverse
#Sobre escribimos el objeto ya antes formado.
dbcovid19mx <- dbcovid19mx %>% mutate(Genero=recode(Genero, `1` = "Femenino",
`2` = "Masculino"),
Paciente=recode(Paciente, `1` = "Ambulatorio",
`2` = "Hospitalizado"),
Neumonia=recode(Neumonia, `1`= "Si",
`2`="No",
`97`="No",
`98`="No",
`99`="No"),
Diabetes=recode(Diabetes, `1`= "Si",
`2`="No",
`97`="No",
`98`="No",
`99`="No"),
Obesidad=recode(Obesidad,`1`= "Si",
`2`="No",
`97`="No",
`98`="No",
`99`="No"),
CLASIFICACION_FINAL=recode(CLASIFICACION_FINAL, `1`="No",
`2`="No",
`3`="si",
`4`="No",
`5`="No",
`6`="No",
`7`="No"))
dbcovid19mx %>% select(Genero, Edad, Paciente, Neumonia, Diabetes, Obesidad)## # A tibble: 500,000 × 6
## Genero Edad Paciente Neumonia Diabetes Obesidad
## <chr> <int> <chr> <chr> <chr> <chr>
## 1 Masculino 47 Ambulatorio No No No
## 2 Masculino 28 Ambulatorio No No No
## 3 Femenino 39 Ambulatorio No No No
## 4 Masculino 8 Ambulatorio No No No
## 5 Femenino 60 Ambulatorio No No No
## 6 Femenino 42 Ambulatorio No No No
## 7 Masculino 44 Ambulatorio No No No
## 8 Masculino 84 Hospitalizado Si No No
## 9 Femenino 27 Ambulatorio No No No
## 10 Masculino 22 Ambulatorio No No No
## # ℹ 499,990 more rows#Siquieres redocificar otras variables, sobre este codigo lo puedes hacer, pero tienes que fijarte que el nombre de la variable sea el adecuado##Pretendemos hacer un analisis de sobrevida y analizar la mortalidad del COVID-19, en este ejercicio en R, vamos a considerar una variable que es la que nos puede dar la pista de que sucedio con el paciente que es representado por cada fila, en este caso la variable “FECHA_DEF”.
tenemos que limpiar la base de datos y eliminar las fechas
9999-99-99, aspi como convertir esta variable en fecha, vamos a indicar
que son espacios vacios. Puede usar el lubridatepaquete para convertir
los valores 9999-99-99 a <NA> y
conservar las fechas ya presentes se conservan.
library(dplyr)
library(lubridate)
# condicionalmente convierte las fechas que son "9999-99-99" a NA
dbcovid19mx <- dbcovid19mx %>%
mutate(FECHA_DEF = if_else(FECHA_DEF == "9999-99-99", NA_character_, FECHA_DEF))
# convierte las fechas restantes al formato deseado
dbcovid19mx$FECHA_DEF <- as.Date(dbcovid19mx$FECHA_DEF, format = "%Y-%m-%d")##como podemos apreciar en la base de datos, tiene el simbolo
##con esta pista, vamos a crear una nueva columna en nuestra base de datos que se va a llamar “estatus”: en la que vamos a definir: ## paciente con fecha de defunción = paciente que No sobrevivio. ##Paciente sin fecha de defunción = NA = Paciente que sobrevivio.
#Survival days
dbcovid19mx <- dbcovid19mx %>% mutate(estatus = case_when(is.na(FECHA_DEF) ~ "Sobrevivió")) #Para la creación de la nueva columna requerimos de la función mutate(), dfiniremos inmediatamente el nombre de nuestra nueva columna además de la función "case_when" es un espacio que esta en NA, con la función "is.na" selecionando en el parentesis cual es la columna que nos interesa y de donde haremos la transformación.
#Primero definiremos la variable estatus: no sobrevivio
dbcovid19mx <- dbcovid19mx %>% replace_na(list(estatus = "No sobrevivió")) #con la función "remplace_na"
dbcovid19mx %>% select(FECHA_DEF, estatus)## # A tibble: 500,000 × 2
## FECHA_DEF estatus
## <date> <chr>
## 1 NA Sobrevivió
## 2 NA Sobrevivió
## 3 NA Sobrevivió
## 4 NA Sobrevivió
## 5 NA Sobrevivió
## 6 NA Sobrevivió
## 7 NA Sobrevivió
## 8 NA Sobrevivió
## 9 NA Sobrevivió
## 10 NA Sobrevivió
## # ℹ 499,990 more rows#Ahora que ya tenemos definida esta variable tan importante, que es la que nos va a determinar la supervivencia, podemos hacer el analisis exploratorio.
##En primer lugar, indagaremos, cuantos pacientes sobrevivieron y cuantos no sobrevivierón. ##Siempre, digo, una gráfica dice más que mil palabras, entonces, haremos un grafico de barras. En el que incluiremos proporciones entre los pacientes que sobrevivieron y no sobrevivieron en esta muestra de cienmil pacientes de la base de datos de COVID-19
#Análisis exploratorio
#En primer lugar, calcularemos proporciones, agrupando, usaremos la función "group_by" y seleccionaremos la variable a agrupar, en este caso "estatus"
estatus <- dbcovid19mx %>%
group_by(estatus) %>%
summarise(counts = n()) %>% #haremos un resumen del conteno usando la función n() y calcularemos las proporciones (porcentaje) con la siguiente formula
mutate(prop = round(counts*100/sum(counts), 1)) # crearemos una nueva columna, con mutate, solicitando un redondeo, multiplicaremos el conteo por 100 entre el total
estatus## # A tibble: 2 × 3
## estatus counts prop
## <chr> <int> <dbl>
## 1 No sobrevivió 3120 0.6
## 2 Sobrevivió 496880 99.4##Pero esto se mira mejor si lo graficamos #Para lo cual usando el paquete que previamente cargamos (en tidyverse) ggplot2, crearemos unagrafica de barras usando el siguiente codigo, iniciamos con la función : ##ggplot,
ggplot(estatus, aes(x = estatus, y = prop, fill= estatus)) + #la función ggplot, indicamos la estetita que queremos que tenga nuestro grafico, en este caso queremos que en el eje de las x, se coloquen las barras que muestran los porcentajes y en y las proporciones divivida por estatus (fill), con el signo más vamos agregando capas a nuestro ggplot:
geom_bar(stat = "identity") + #agregaremos los datos con la función geom_bar que nos ayuda a definir que tipo de grafico en este caso es un grafico de barras
geom_text(aes(label = prop), vjust = -0.3, fontface=2, size= 6) +
theme_pubclean() + ylab("%") + xlab ("Mortalidad")+
theme(text = element_text(size = 16, face="bold"), axis.text = element_text(size = 16, face="bold"),
legend.text = element_text(size = 12),)+ scale_fill_brewer(palette="Dark2")+
theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(legend.position = "none")
¿Qué te parece sí esto lo vemos en una curva?
Primero, identificamos dos cosas importantes. La variable que continene la fecha de defuncion, esta la colocaremos en el eje de las x, con el paquete: lubridate. podemos usar funciones como scale_x_date, generando puntos de corte ( “2 week”) fue el que se uso y las etiquetas %d = día, %b= mes y %y = año. Dejaremos solo el mes y el año.
library(lubridate)
ggplot(data = dbcovid19mx, aes(x = as.Date(FECHA_DEF))) +
theme_bw() +
geom_histogram(binwidth = 2, colour = "black", fill = "#3e78d6", size = 0.2)+
geom_density(aes(y = ..density.. * (nrow(dbcovid19mx) * 0.01)), colour = "red", size= .6)+
scale_x_date(date_breaks = "2 week",
name = "Date") +
scale_y_continuous(breaks = seq(from = 0, to =100, by = 20), name = "Number of cases") + theme_classic()+
theme(axis.text.x = element_text(angle = 90, hjust = 1))## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Removed 496880 rows containing non-finite values (`stat_bin()`).## Warning: Removed 496880 rows containing non-finite values (`stat_density()`).
glimpse(covid19mx)## Rows: 5,705,577
## Columns: 40
## $ FECHA_ACTUALIZACION <date> 2022-09-27, 2022-09-27, 2022-09-27, 2022-09-27,…
## $ ID_REGISTRO <chr> "180725", "06fce8", "1a4a8d", "1e6dad", "06e951"…
## $ ORIGEN <int> 2, 1, 1, 1, 1, 1, 1, 2, 1, 2, 2, 2, 1, 2, 2, 2, …
## $ SECTOR <int> 9, 12, 12, 12, 12, 12, 12, 6, 12, 12, 12, 6, 12,…
## $ ENTIDAD_UM <int> 9, 7, 23, 7, 24, 22, 16, 9, 18, 9, 9, 9, 20, 14,…
## $ SEXO <int> 2, 1, 2, 1, 2, 1, 1, 1, 1, 1, 1, 2, 2, 1, 2, 2, …
## $ ENTIDAD_NAC <int> 9, 7, 27, 6, 24, 9, 16, 9, 18, 9, 9, 9, 20, 32, …
## $ ENTIDAD_RES <int> 9, 7, 23, 7, 24, 22, 16, 9, 18, 9, 9, 9, 20, 14,…
## $ MUNICIPIO_RES <int> 12, 59, 8, 12, 28, 14, 67, 15, 7, 2, 7, 16, 67, …
## $ TIPO_PACIENTE <int> 2, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, …
## $ FECHA_INGRESO <date> 2022-01-19, 2022-03-05, 2022-07-22, 2022-08-23,…
## $ FECHA_SINTOMAS <date> 2022-01-17, 2022-03-05, 2022-07-20, 2022-08-21,…
## $ FECHA_DEF <chr> "9999-99-99", "9999-99-99", "9999-99-99", "9999-…
## $ INTUBADO <int> 2, 97, 97, 97, 97, 97, 2, 97, 97, 97, 97, 2, 97,…
## $ NEUMONIA <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, …
## $ EDAD <int> 33, 35, 49, 34, 40, 56, 69, 52, 24, 41, 28, 54, …
## $ NACIONALIDAD <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ EMBARAZO <int> 97, 98, 97, 2, 97, 2, 2, 2, 2, 2, 2, 97, 97, 2, …
## $ HABLA_LENGUA_INDIG <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
## $ INDIGENA <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
## $ DIABETES <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
## $ EPOC <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
## $ ASMA <int> 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, …
## $ INMUSUPR <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
## $ HIPERTENSION <int> 2, 2, 2, 2, 1, 2, 1, 2, 2, 1, 2, 1, 2, 2, 2, 2, …
## $ OTRA_COM <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
## $ CARDIOVASCULAR <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
## $ OBESIDAD <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, …
## $ RENAL_CRONICA <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, …
## $ TABAQUISMO <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, …
## $ OTRO_CASO <int> 2, 2, 2, 2, 2, 1, 2, 2, 2, 1, 1, 2, 1, 2, 2, 2, …
## $ TOMA_MUESTRA_LAB <int> 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 1, 2, 2, 1, 2, …
## $ RESULTADO_LAB <int> 97, 97, 97, 97, 97, 2, 97, 97, 97, 97, 97, 4, 97…
## $ TOMA_MUESTRA_ANTIGENO <int> 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 1, …
## $ RESULTADO_ANTIGENO <int> 2, 2, 2, 97, 2, 2, 2, 97, 1, 2, 2, 97, 2, 1, 2, …
## $ CLASIFICACION_FINAL <int> 7, 7, 7, 6, 7, 7, 7, 6, 3, 7, 7, 2, 7, 3, 7, 7, …
## $ MIGRANTE <int> 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, …
## $ PAIS_NACIONALIDAD <chr> "México", "México", "México", "México", "México"…
## $ PAIS_ORIGEN <chr> "97", "97", "97", "97", "97", "97", "97", "97", …
## $ UCI <int> 2, 97, 97, 97, 97, 97, 2, 97, 97, 97, 97, 2, 97,…Ahora veremos una curva clasificados entre sobrevivientes y no sobrevivientes, basandonos en la fecha de inicio de los sintomas:
library(lubridate)
library(epiR); library(ggplot2); library(scales); library(zoo)## Warning: package 'epiR' was built under R version 4.2.3## Loading required package: survival## Package epiR 2.0.59 is loaded## Type help(epi.about) for summary information## Type browseVignettes(package = 'epiR') to learn how to use epiR for applied epidemiological analyses## ##
## Attaching package: 'scales'## The following object is masked from 'package:purrr':
##
## discard## The following object is masked from 'package:readr':
##
## col_factor##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(plotly)##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:arrow':
##
## schema## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutggplot(data = dbcovid19mx, aes(x = as.Date(FECHA_SINTOMAS), group=estatus, fill=estatus)) +
theme_bw() +
geom_histogram(binwidth = 7, colour = "black", size = 0.1)+
scale_x_date(date_breaks = "2 week",labels = date_format("%d %b"),
name = "Date") +
scale_y_continuous(breaks = seq(from = 0, to =100000, by = 5000), name = "Number of cases") + theme_classic()+
theme(axis.text.x = element_text(angle = 90, hjust = 1)) +scale_fill_manual(values = c("#d46a6a", "#738ca6"), name = "") +
geom_vline(aes(xintercept = as.numeric(as.Date("01/05/2021", format = "%d/%m/%Y"))),
linetype = "dashed") +
theme(legend.text = element_text(size = 12, color = "black", face = "bold"))+theme(text = element_text(size = 18)) 
##Parte del analisis exploratorio, es la creación de tablas en las cuales podemos incluir valores estadísticos como medias, medianas desviaciones estandar, intervalos de confianza y valores p ##utilizaremos el paquete gtsummary, que ya abrimos previamente cuando cargamos los paquetes para crear una tabla en la cual describamos las variables que hemos seleccionado y son de nnuestro interes con la función “tbl_summary()”, seleccionamos previamente que variables nos interesa explorar
dbcovid19mx %>% select(Genero, Edad, Paciente, Neumonia, Diabetes, Obesidad, estatus) %>% tbl_summary()| Characteristic | N = 500,0001 |
|---|---|
| Genero | |
| Femenino | 283,450 (57%) |
| Masculino | 216,550 (43%) |
| Edad | 36 (26, 49) |
| Paciente | |
| Ambulatorio | 481,352 (96%) |
| Hospitalizado | 18,648 (3.7%) |
| Neumonia | |
| No | 492,109 (98%) |
| Si | 7,891 (1.6%) |
| Diabetes | |
| No | 470,285 (94%) |
| Si | 29,715 (5.9%) |
| Obesidad | |
| No | 470,504 (94%) |
| Si | 29,496 (5.9%) |
| estatus | |
| No sobrevivió | 3,120 (0.6%) |
| Sobrevivió | 496,880 (99%) |
| 1 n (%); Median (IQR) | |
##Podemos hacer una división en columnas, considerando la variable estatus, que nos define que pacientes sobervivieron y no, usando la función tbl_summary y dividiendo con la función “by”= colocamos la columna que sería nuestra variable predictora o dependiente
dbcovid19mx %>% select(Genero, Edad, Paciente, Neumonia, Diabetes, Obesidad, CLASIFICACION_FINAL, estatus) %>% tbl_summary(by=estatus) %>% add_p()| Characteristic | No sobrevivió, N = 3,1201 | Sobrevivió, N = 496,8801 | p-value2 |
|---|---|---|---|
| Genero | <0.001 | ||
| Femenino | 1,307 (42%) | 282,143 (57%) | |
| Masculino | 1,813 (58%) | 214,737 (43%) | |
| Edad | 69 (57, 80) | 36 (26, 49) | <0.001 |
| Paciente | <0.001 | ||
| Ambulatorio | 88 (2.8%) | 481,264 (97%) | |
| Hospitalizado | 3,032 (97%) | 15,616 (3.1%) | |
| Neumonia | <0.001 | ||
| No | 1,427 (46%) | 490,682 (99%) | |
| Si | 1,693 (54%) | 6,198 (1.2%) | |
| Diabetes | <0.001 | ||
| No | 1,982 (64%) | 468,303 (94%) | |
| Si | 1,138 (36%) | 28,577 (5.8%) | |
| Obesidad | <0.001 | ||
| No | 2,742 (88%) | 467,762 (94%) | |
| Si | 378 (12%) | 29,118 (5.9%) | |
| CLASIFICACION_FINAL | <0.001 | ||
| No | 1,093 (35%) | 242,694 (49%) | |
| si | 2,027 (65%) | 254,186 (51%) | |
| 1 n (%); Median (IQR) | |||
| 2 Pearson's Chi-squared test; Wilcoxon rank sum test | |||
#Con esto ya organizado podemos hacer analisis estadísticos basandonos en comparaciones utilizando pruebas adecuadas.
#El paquete gtsummary, nos permite además agregar valores de p, basandose en los supuestos de normalidad, homocedasticidad para la selección d ela prueba adecuada##LIMPIEZA DE DATOS ##Como pueden darse cuenta, hay algunas variables que tienen espacios vacios, esos debemos de limpiarlos, antes de disponernos a hacer analisis de tipo estadísticos ##Vamos a crear una nueva variable que se llame filtradad covid19mxfilt, ##Además cabe señalar que existe una variable en la base de datos que va a ser crucial para saber que solo estamos incluyendo a pacientes que tuvieron COVID-19, pues esta se define como aquellos que fueron positivos a SARS-COV2 confirmado por PCR, excluiremos entonces a los pacientes que fueron negativos, usando la función “filter ()”
#filtrar espacios vacios
#excluimos pacientes negativos
covid19mxfilt <- dbcovid19mx %>% filter(Neumonia != "" & Diabetes != "" & Obesidad !="" & CLASIFICACION_FINAL== "si")
covid19mxfilt## # A tibble: 256,213 × 41
## FECHA_ACTUALIZACION ID_REGISTRO ORIGEN SECTOR ENTIDAD_UM Genero ENTIDAD_NAC
## <date> <chr> <int> <int> <int> <chr> <int>
## 1 2022-09-27 g126507 2 4 16 Masculi… 16
## 2 2022-09-27 g0d8fe5 2 4 15 Femenino 15
## 3 2022-09-27 gz56518 2 12 9 Masculi… 9
## 4 2022-09-27 g14393e 1 4 11 Femenino 11
## 5 2022-09-27 g11716a 2 4 13 Masculi… 13
## 6 2022-09-27 9a3816 2 4 26 Masculi… 26
## 7 2022-09-27 g30a96c 2 4 19 Femenino 19
## 8 2022-09-27 g106b9e 2 4 14 Masculi… 14
## 9 2022-09-27 g08a261 2 12 29 Femenino 29
## 10 2022-09-27 804e94 2 12 28 Femenino 30
## # ℹ 256,203 more rows
## # ℹ 34 more variables: ENTIDAD_RES <int>, MUNICIPIO_RES <int>, Paciente <chr>,
## # FECHA_INGRESO <date>, FECHA_SINTOMAS <date>, FECHA_DEF <date>,
## # INTUBADO <int>, Neumonia <chr>, Edad <int>, NACIONALIDAD <int>,
## # EMBARAZO <int>, HABLA_LENGUA_INDIG <int>, INDIGENA <int>, Diabetes <chr>,
## # EPOC <int>, ASMA <int>, INMUSUPR <int>, HIPERTENSION <int>, OTRA_COM <int>,
## # CARDIOVASCULAR <int>, Obesidad <chr>, RENAL_CRONICA <int>, …#con esta acción se reducira el tamaño de muestra a = 256,962 pacientes
#volvemos a hacer la tabla pero con el nuevo objeto que acabos de crear
covid19mxfilt %>% select(Genero, Edad, Paciente, Neumonia, Diabetes, Obesidad,CLASIFICACION_FINAL, estatus) %>% tbl_summary(by=estatus) %>% add_overall() %>% add_p()## There was an error in 'add_p()/add_difference()' for variable 'CLASIFICACION_FINAL', p-value omitted:
## Error in stats::chisq.test(x = c("si", "si", "si", "si", "si", "si", "si", : 'x' and 'y' must have at least 2 levels| Characteristic | Overall, N = 256,2131 | No sobrevivió, N = 2,0271 | Sobrevivió, N = 254,1861 | p-value2 |
|---|---|---|---|---|
| Genero | <0.001 | |||
| Femenino | 146,280 (57%) | 849 (42%) | 145,431 (57%) | |
| Masculino | 109,933 (43%) | 1,178 (58%) | 108,755 (43%) | |
| Edad | 37 (27, 49) | 71 (58, 81) | 37 (27, 49) | <0.001 |
| Paciente | <0.001 | |||
| Ambulatorio | 249,045 (97%) | 0 (0%) | 249,045 (98%) | |
| Hospitalizado | 7,168 (2.8%) | 2,027 (100%) | 5,141 (2.0%) | |
| Neumonia | <0.001 | |||
| No | 251,550 (98%) | 775 (38%) | 250,775 (99%) | |
| Si | 4,663 (1.8%) | 1,252 (62%) | 3,411 (1.3%) | |
| Diabetes | <0.001 | |||
| No | 240,989 (94%) | 1,271 (63%) | 239,718 (94%) | |
| Si | 15,224 (5.9%) | 756 (37%) | 14,468 (5.7%) | |
| Obesidad | <0.001 | |||
| No | 239,353 (93%) | 1,787 (88%) | 237,566 (93%) | |
| Si | 16,860 (6.6%) | 240 (12%) | 16,620 (6.5%) | |
| CLASIFICACION_FINAL | ||||
| si | 256,213 (100%) | 2,027 (100%) | 254,186 (100%) | |
| 1 n (%); Median (IQR) | ||||
| 2 Pearson's Chi-squared test; Wilcoxon rank sum test | ||||
#Ahora así, podemos hacer analisis estadístico agregando la función “add_p ()”, que nos va a agregar el valor de p a la tabla para las diferentes comparaciones
covid19mxfilt %>% select(Genero, Edad, Paciente, Neumonia, Diabetes, Obesidad, estatus) %>% tbl_summary(by=estatus) %>% add_overall() %>% add_p()| Characteristic | Overall, N = 256,2131 | No sobrevivió, N = 2,0271 | Sobrevivió, N = 254,1861 | p-value2 |
|---|---|---|---|---|
| Genero | <0.001 | |||
| Femenino | 146,280 (57%) | 849 (42%) | 145,431 (57%) | |
| Masculino | 109,933 (43%) | 1,178 (58%) | 108,755 (43%) | |
| Edad | 37 (27, 49) | 71 (58, 81) | 37 (27, 49) | <0.001 |
| Paciente | <0.001 | |||
| Ambulatorio | 249,045 (97%) | 0 (0%) | 249,045 (98%) | |
| Hospitalizado | 7,168 (2.8%) | 2,027 (100%) | 5,141 (2.0%) | |
| Neumonia | <0.001 | |||
| No | 251,550 (98%) | 775 (38%) | 250,775 (99%) | |
| Si | 4,663 (1.8%) | 1,252 (62%) | 3,411 (1.3%) | |
| Diabetes | <0.001 | |||
| No | 240,989 (94%) | 1,271 (63%) | 239,718 (94%) | |
| Si | 15,224 (5.9%) | 756 (37%) | 14,468 (5.7%) | |
| Obesidad | <0.001 | |||
| No | 239,353 (93%) | 1,787 (88%) | 237,566 (93%) | |
| Si | 16,860 (6.6%) | 240 (12%) | 16,620 (6.5%) | |
| 1 n (%); Median (IQR) | ||||
| 2 Pearson's Chi-squared test; Wilcoxon rank sum test | ||||
##Una vez teniendo encuenta estos datos como parte del analisis exploratorio haremos el analisis de sobrevida basado en curvas de kaplan Meier (escribir) ##Para el siguiente analisis requeriremos otros paquetes:
library("survival")
library("ggpubr")
library("survminer")##
## Attaching package: 'survminer'## The following object is masked from 'package:survival':
##
## myeloma##Lo importante en el analisis de supervivencia es tener una variable que defina los días, el tiempo de supervivencia, para lo cual tomaremos de referencia dos columnas ##la fecga de defunción (FECHA_DEF) y fecha de inicio de sintomas, que marca el día 0 en el que se inició la enfermedad (FECHA_SINTOMAS) ##Creamos una nueva columna, entonces, usando la función mutate y a la cual nombraremos “futime”, es una variable númerica por lo tanto a defininoms como as.double y esta creada por la diferencia entre la fecha de defunción y la fecha de inicio de sintomas
covid19mxfilt <- covid19mxfilt %>% mutate(futime = as.double(FECHA_DEF-FECHA_SINTOMAS))
summary(covid19mxfilt$futime)#con la función summary podemos hacer un analisis exploratorio rápido para encontrar los estadisticos que se muestran## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.00 5.00 10.00 11.67 16.00 83.00 254186##Ademas de definir el tiempo de supervivencia, es importante también definir la variable descenlace y codificarla como 0 y 1, en donde el 0 será igual a Sobrevivió y 1 igual a no sobrevivió, esto se hace creando una nueva columna a la cual llamaremos “evento”
# create the event var which is 1 if the patient died and 0 if he was right censored
covid19mxfilt <- covid19mxfilt %>% mutate(evento = ifelse(is.na(estatus) | estatus == "Sobrevivió", 0, 1))
#cross tabulate the new event var and the outcome var from which it was created
# to make sure the code did what it was intended to
covid19mxfilt %>% select(estatus, evento, futime)## # A tibble: 256,213 × 3
## estatus evento futime
## <chr> <dbl> <dbl>
## 1 Sobrevivió 0 NA
## 2 Sobrevivió 0 NA
## 3 Sobrevivió 0 NA
## 4 Sobrevivió 0 NA
## 5 Sobrevivió 0 NA
## 6 Sobrevivió 0 NA
## 7 Sobrevivió 0 NA
## 8 Sobrevivió 0 NA
## 9 Sobrevivió 0 NA
## 10 Sobrevivió 0 NA
## # ℹ 256,203 more rowscovid19mxfilt %>% tabyl(evento, estatus)## evento No sobrevivió Sobrevivió
## 0 0 254186
## 1 2027 0ANALISIS DE SUPERVIVENCIA
graficas de kaplan-Meier
covidmxsuv <- covid19mxfilt
total<-survfit(Surv(futime, evento) ~ 1, data = covidmxsuv)
total## Call: survfit(formula = Surv(futime, evento) ~ 1, data = covidmxsuv)
##
## 254186 observations deleted due to missingness
## n events median 0.95LCL 0.95UCL
## [1,] 2027 2027 10 9 10##Curvas de kaplan meier
ggsurvplot(
fit = survfit(Surv(futime, evento) ~ 1, data = covidmxsuv),
xlab = "Días",
ylab = "Probabilidad de supervivencia total",
surv.median.line = "hv")
##Podemos definir que probabilidad existion a no sobrevivir si los pacientes se encontraban con diferentes condiciones por ejemplo la presencia de neumonía, diabetes u obesidad, así como en genero
gen <- survfit(Surv(futime, evento) ~ Genero, data=covidmxsuv)
neu <- survfit(Surv(futime, evento) ~ Neumonia, data = covidmxsuv)
dm <- survfit( Surv(futime, evento) ~ Diabetes, data = covidmxsuv)
obe <- survfit(Surv(futime, evento) ~ Obesidad , data = covidmxsuv)##Genero
gen## Call: survfit(formula = Surv(futime, evento) ~ Genero, data = covidmxsuv)
##
## 254186 observations deleted due to missingness
## n events median 0.95LCL 0.95UCL
## Genero=Femenino 849 849 9 9 10
## Genero=Masculino 1178 1178 10 9 11##Curva KM y valor de p para la diferencia entre curvas que depende del genero, agregamos valor de p con la función “pval”
library(ggthemes)
ggsurvplot(gen,
xlab = "Días",
ylab = "Probabilidad de supervivencia total",
surv.median.line = "hv",
pval= TRUE, pval.coord = c(100, 0.25),risk.table = TRUE,fontsize = 5, break.time.by = 20, legend.title = "Genero", legend.labs = c("Femenino", "Másculino"))
##Neumonia
neu## Call: survfit(formula = Surv(futime, evento) ~ Neumonia, data = covidmxsuv)
##
## 254186 observations deleted due to missingness
## n events median 0.95LCL 0.95UCL
## Neumonia=No 775 775 9 8 9
## Neumonia=Si 1252 1252 10 10 11##Curva KM y valor de p para la diferencia entre curvas que depende del genero, agregamos valor de p con la función “pval”
library(ggthemes)
ggsurvplot(neu,
xlab = "Días",
ylab = "Probabilidad de supervivencia total",
surv.median.line = "hv",
pval= TRUE, pval.coord = c(100, 0.25),risk.table = TRUE,fontsize = 5, break.time.by = 20, legend.title = "Neumonia", legend.labs = c("No", "Si"))
##diabetes
dm## Call: survfit(formula = Surv(futime, evento) ~ Diabetes, data = covidmxsuv)
##
## 254186 observations deleted due to missingness
## n events median 0.95LCL 0.95UCL
## Diabetes=No 1271 1271 10 9 10
## Diabetes=Si 756 756 9 9 10##Curva KM y valor de p para la diferencia entre curvas que depende del genero, agregamos valor de p con la función “pval”
library(ggthemes)
ggsurvplot(dm,
xlab = "Días",
ylab = "Probabilidad de supervivencia total",
surv.median.line = "hv",
pval= TRUE, pval.coord = c(100, 0.25),risk.table = TRUE,fontsize = 5, break.time.by = 20, legend.title = "Diabetes", legend.labs = c("No", "Si"))
##Obesidad
obe## Call: survfit(formula = Surv(futime, evento) ~ Obesidad, data = covidmxsuv)
##
## 254186 observations deleted due to missingness
## n events median 0.95LCL 0.95UCL
## Obesidad=No 1787 1787 10 9 10
## Obesidad=Si 240 240 10 10 12##Curva KM y valor de p para la diferencia entre curvas que depende del genero, agregamos valor de p con la función “pval”
ggsurvplot(obe,
xlab = "Días",
ylab = "Probabilidad de supervivencia total",
surv.median.line = "hv",
pval= TRUE, pval.coord = c(100, 0.25),risk.table = TRUE,fontsize = 5, break.time.by = 20, legend.title = "Obesidad", legend.labs = c("No", "Si"))
##Cox proportional hazard model
#fit the model
linelistsurv_cox <- coxph(Surv(futime, evento) ~ Genero + Edad + Neumonia + Diabetes + Obesidad, data = covidmxsuv)
linelistsurv_cox## Call:
## coxph(formula = Surv(futime, evento) ~ Genero + Edad + Neumonia +
## Diabetes + Obesidad, data = covidmxsuv)
##
## coef exp(coef) se(coef) z p
## GeneroMasculino -5.568e-02 9.458e-01 4.518e-02 -1.232 0.21780
## Edad 8.788e-05 1.000e+00 1.313e-03 0.067 0.94662
## NeumoniaSi -1.502e-01 8.605e-01 4.599e-02 -3.266 0.00109
## DiabetesSi 8.164e-02 1.085e+00 4.631e-02 1.763 0.07790
## ObesidadSi -1.432e-01 8.665e-01 6.970e-02 -2.055 0.03986
##
## Likelihood ratio test=19.85 on 5 df, p=0.001333
## n= 2027, number of events= 2027
## (254186 observations deleted due to missingness)##HR model
summary(linelistsurv_cox)## Call:
## coxph(formula = Surv(futime, evento) ~ Genero + Edad + Neumonia +
## Diabetes + Obesidad, data = covidmxsuv)
##
## n= 2027, number of events= 2027
## (254186 observations deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## GeneroMasculino -5.568e-02 9.458e-01 4.518e-02 -1.232 0.21780
## Edad 8.788e-05 1.000e+00 1.313e-03 0.067 0.94662
## NeumoniaSi -1.502e-01 8.605e-01 4.599e-02 -3.266 0.00109 **
## DiabetesSi 8.164e-02 1.085e+00 4.631e-02 1.763 0.07790 .
## ObesidadSi -1.432e-01 8.665e-01 6.970e-02 -2.055 0.03986 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## GeneroMasculino 0.9458 1.0573 0.8657 1.0334
## Edad 1.0001 0.9999 0.9975 1.0027
## NeumoniaSi 0.8605 1.1621 0.7864 0.9417
## DiabetesSi 1.0851 0.9216 0.9909 1.1882
## ObesidadSi 0.8665 1.1540 0.7559 0.9934
##
## Concordance= 0.539 (se = 0.008 )
## Likelihood ratio test= 19.85 on 5 df, p=0.001
## Wald test = 19.85 on 5 df, p=0.001
## Score (logrank) test = 19.88 on 5 df, p=0.001library(forestmodel)
print(forest_model(coxph(Surv(futime, evento) ~ Genero + Edad + Neumonia + Diabetes + Obesidad, data = covidmxsuv)))## Warning in recalculate_width_panels(panel_positions, mapped_text = mapped_text,
## : Unable to resize forest panel to be smaller than its heading; consider a
## smaller text size
##Corroboramos el modelo
#test the proportional hazard model
linelistsurv_ph_test <- cox.zph(linelistsurv_cox)
linelistsurv_ph_test## chisq df p
## Genero 0.136 1 0.7127
## Edad 5.208 1 0.0225
## Neumonia 10.648 1 0.0011
## Diabetes 0.600 1 0.4386
## Obesidad 1.528 1 0.2164
## GLOBAL 17.344 5 0.0039#The graphical verification of this assumption may be performed with the function ggcoxzph() from the survminer package.
survminer::ggcoxzph(linelistsurv_ph_test)