### Efecto de la demografia y conteos de celulas T en la recuperacion
### de pacientes con sindrome post-COVID-19
## Lopez Sanchez Jaciel Ivan
# 25/05/2024
# Cargar librerias
library("FactoMineR")
library("factoextra")
## Cargando paquete requerido: ggplot2
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
library("corrplot")
## corrplot 0.92 loaded
library(ggplot2)
library(plyr)
library(scales)
library(grid)
library(devtools)
## Cargando paquete requerido: usethis
library(ggbiplot)
library(car)
## Cargando paquete requerido: carData
library(broom)
library(readr)
##
## Adjuntando el paquete: 'readr'
## The following object is masked from 'package:scales':
##
## col_factor
library(dplyr)
##
## Adjuntando el paquete: 'dplyr'
## The following object is masked from 'package:car':
##
## recode
## The following objects are masked from 'package:plyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tinytex)
getwd()
## [1] "C:/Users/rinco/Downloads"
# Cargar el conjunto de datos
df <- read.csv("C:/Users/rinco/Desktop/01_Base/02_Data/pacientes.csv", header = TRUE, sep = ",")
df2 <- read.csv("C:/Users/rinco/Desktop/01_Base/02_Data/pacientes2.csv", header = TRUE, sep = ",")
### Analisis de Componentes Principales
# PCA Best
#################
# Realiza el PCA en las columnas 4 a 6 del DataFrame
WS.pca <- PCA(df [, 4:6], graph = FALSE)
# Obtiener los valores propios (eigenvalues) del PCA
EigenCOVIDClinic <- get_eig(WS.pca)
# Muestra los valores propios en la consola
EigenCOVIDClinic
## eigenvalue variance.percent cumulative.variance.percent
## Dim.1 1.4135428 47.11809 47.11809
## Dim.2 1.0037889 33.45963 80.57772
## Dim.3 0.5826683 19.42228 100.00000
# Guardar los valores propios en un archivo CSV llamado 'Eigen_pacientes.csv'
write.table(EigenCOVIDClinic, file = "02_Data/Eigen_pacientes.csv",
sep = ",", eol = "\r\n", na = "NA",
col.names = TRUE, row.names = FALSE, fileEncoding = "")
# Genera un scree plot del PCA
fviz_screeplot(WS.pca, addlabels=TRUE, ylim = c(0, 49))
# Obtiener las variables del PCA
var <- get_pca_var(WS.pca)
# Mostrar las variables del PCA en la consola
var
## Principal Component Analysis Results for variables
## ===================================================
## Name Description
## 1 "$coord" "Coordinates for the variables"
## 2 "$cor" "Correlations between variables and dimensions"
## 3 "$cos2" "Cos2 for the variables"
## 4 "$contrib" "contributions of the variables"
# Primeras filas de las coordenadas de las variables en el PCA
head(var$coord)
## Dim.1 Dim.2 Dim.3
## BhLinfocitosAbsolutos -0.8426088 -0.002566662 0.5385200
## edad 0.5048812 0.798508478 0.3278402
## peso 0.6698120 -0.605116937 0.4303317
# Primeras filas de los valores cos2
head(var$cos2)
## Dim.1 Dim.2 Dim.3
## BhLinfocitosAbsolutos 0.7099896 6.587753e-06 0.2900038
## edad 0.2549050 6.376158e-01 0.1074792
## peso 0.4486481 3.661665e-01 0.1851854
# Primeras filas de la contribucion de las variables a los componentes principales
head(var$contrib)
## Dim.1 Dim.2 Dim.3
## BhLinfocitosAbsolutos 50.22767 6.562887e-04 49.77167
## edad 18.03306 6.352091e+01 18.44604
## peso 31.73927 3.647844e+01 31.78229
# Almacenar la contribución de las variables en 'Contrib'
Contrib <- var$contrib
# Guardar la contribucion de las variables
write.table(Contrib, file = "02_Data/contrib_pacientes.csv",
sep = ",", eol = "\r\n", na = "NA",
col.names = TRUE, row.names = FALSE, fileEncoding = "")
# Genera un gráfico de las variables del PCA
fviz_pca_var(WS.pca, col.var = "black")
# Mostrar las primeras 4 filas de los valores cos2
head(var$cos2, 4)
## Dim.1 Dim.2 Dim.3
## BhLinfocitosAbsolutos 0.7099896 6.587753e-06 0.2900038
## edad 0.2549050 6.376158e-01 0.1074792
## peso 0.4486481 3.661665e-01 0.1851854
# Generar un gráfico de correlación de los valores cos2 de las variables
corrplot(var$cos2, is.corr=F, method = "number")
# Generar un gráfico de contribución de las variables a los componentes principales
fviz_contrib(WS.pca, choice = "var", axes = 1, top = 35)
fviz_contrib(WS.pca, choice = "var", axes = 2, top = 35)
# Generar un gráfico de calidad de representación (cos2) de las variables en los componentes principales
fviz_cos2(WS.pca, choice = "var", axes = 1:2)
### Estadistica descriptiva
summary(df)
## Id Recuperacion_Postcovid sexo BhLinfocitosAbsolutos
## Min. : 1.00 Length:102 Min. :1.000 Min. :0.230
## 1st Qu.: 26.25 Class :character 1st Qu.:1.000 1st Qu.:0.780
## Median : 52.50 Mode :character Median :1.000 Median :1.220
## Mean : 52.05 Mean :1.451 Mean :1.414
## 3rd Qu.: 77.75 3rd Qu.:2.000 3rd Qu.:1.802
## Max. :103.00 Max. :2.000 Max. :4.310
## edad peso
## Min. :16.00 Min. : 50.00
## 1st Qu.:41.25 1st Qu.: 65.00
## Median :50.00 Median : 77.00
## Mean :49.40 Mean : 78.85
## 3rd Qu.:58.00 3rd Qu.: 90.00
## Max. :90.00 Max. :122.00
# Realizar un boxplot con los datos para observar la distribucion de los datos
# Recuperacion post-COVID-19 vs conteos de linfocitos absolutos
ggplot(df, aes(Recuperacion_Postcovid, BhLinfocitosAbsolutos))+geom_boxplot()+
labs(x = "Recuperacion post-COVID-19",
y= "Conteo de linfocitos absolutos") +
theme_minimal()
# Recuperacion post-COVID-19 vs Edad
ggplot(df, aes(Recuperacion_Postcovid, edad))+geom_boxplot()+
labs(x = "Recuperacion post-COVID-19",
y= "Edad (años)") +
theme_minimal()
# Recuperacion post-COVID-19 vs peso
ggplot(df, aes(Recuperacion_Postcovid, peso))+geom_boxplot()+
labs(x = "Recuperacion post-COVID-19",
y= "Pesp (kg)") +
theme_minimal()
# Se comparara la magnitid de las tres variables anteriores en ambos grupos
### Supuestos para la prueba de t
## Normalidad
# Se subdiviran los datos para la prueba de normalidad
# Conteos de linfocitos
LNR <-subset(df, Recuperacion_Postcovid=="No recuperado", select = "BhLinfocitosAbsolutos")
LR <-subset(df, Recuperacion_Postcovid=="Recuperado", select = "BhLinfocitosAbsolutos")
# Edad
ENR <-subset(df, Recuperacion_Postcovid=="No recuperado", select = "edad")
ER <-subset(df, Recuperacion_Postcovid=="Recuperado", select = "edad")
# Peso
PNR <-subset(df, Recuperacion_Postcovid=="No recuperado", select = "peso")
PR <-subset(df, Recuperacion_Postcovid=="Recuperado", select = "peso")
# Realizar la prueba de Shapiro-Wilk para estos subconjuntos de datos
# prueba de decicion, se utilizara un valor de confianza de a=0.05
# prueba de hipotesis: H0: los datos siguen una distribucion normal
# HA: La distribucion de los datos no es normal
## Conteos de linfocitos
shapiro.test(LNR$BhLinfocitosAbsolutos) # p = 0.02979, los datos respaldan HA
##
## Shapiro-Wilk normality test
##
## data: LNR$BhLinfocitosAbsolutos
## W = 0.85148, p-value = 0.02979
shapiro.test(LR$BhLinfocitosAbsolutos) # p = 3.424e-06, los datos respaldan HA
##
## Shapiro-Wilk normality test
##
## data: LR$BhLinfocitosAbsolutos
## W = 0.89729, p-value = 3.424e-06
# Ya que ambas subdivisiones del data frame muestran una distribucion distinta a la normal, se utilizaran pruebas no parametricas para realizar las comparaciones.
## Edad
shapiro.test(ENR$edad) # p = 0.07219, los datos respaldan H0
##
## Shapiro-Wilk normality test
##
## data: ENR$edad
## W = 0.88039, p-value = 0.07219
shapiro.test(ER$edad) # p = 0.1943, los datos respaldan H0
##
## Shapiro-Wilk normality test
##
## data: ER$edad
## W = 0.98024, p-value = 0.1934
# Ya que ambas subdivisiones del data frame muestran una distribucion normal, se utilizaran pruebas parametricas para realizar las comparaciones.
## Peso
shapiro.test(PNR$peso) # p = 0.7673, los datos respaldan HA
##
## Shapiro-Wilk normality test
##
## data: PNR$peso
## W = 0.96089, p-value = 0.7673
shapiro.test(PR$peso) # p = 0.01473, los datos respaldan H0
##
## Shapiro-Wilk normality test
##
## data: PR$peso
## W = 0.96417, p-value = 0.01473
# Ya que ambas subdivisiones del data frame muestran distribuciones distintas, se utilizaran pruebas no parametricas para realizar las comparaciones.
## Muestras balanceadas
with(df, table(df$Recuperacion_Postcovid))
##
## No recuperado Recuperado
## 13 89
# Las muestras no estan balanceadas
## homogeneidad de varianzas
# Se utilizara la prueba de Levene dado a que los factores presentan dos niveles
# Prueba de hipotesis: H0: las varianzas son iguales en todos los grupos; HA: Al menos una de las varianzas es diferente.
## Conteos de linfocitos
leveneTest(BhLinfocitosAbsolutos ~ Recuperacion_Postcovid, data = df) # p= 0.5274, los datos respaldan H0
## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.4023 0.5274
## 100
## Edad
leveneTest(edad ~ Recuperacion_Postcovid, data = df) # p= 0.02306, los datos respaldan HA
## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 1.4548 0.2306
## 100
## Peso
leveneTest(peso ~ Recuperacion_Postcovid, data = df) # p= 0.6831, los datos respaldan H0
## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.1676 0.6831
## 100
# Se empleara la prueba de Kruskal-Wallis para comparar los conteos de linfocitos B y la recuperacion
# Prueba de hipotesis, H0: La distribucion de las muestras en todos los grupos es identica
# HA: Al menos una de las muestras en los grupos difiere en su distribucion de las demas
kruskal.test(BhLinfocitosAbsolutos ~ Recuperacion_Postcovid, data = df) # p = 0.3481, los datos respaldan la H0
##
## Kruskal-Wallis rank sum test
##
## data: BhLinfocitosAbsolutos by Recuperacion_Postcovid
## Kruskal-Wallis chi-squared = 0.88042, df = 1, p-value = 0.3481
# Se empleara la prueba de Kruskal-Wallis para comparar el peso y la recuperacion
# Prueba de hipotesis, H0: La distribucion de las muestras en todos los grupos es identica
# HA: Al menos una de las muestras en los grupos difiere en su distribucion de las demas
kruskal.test(edad ~ Recuperacion_Postcovid, data = df) # p = 0.6842, los datos respaldan la H0
##
## Kruskal-Wallis rank sum test
##
## data: edad by Recuperacion_Postcovid
## Kruskal-Wallis chi-squared = 0.16538, df = 1, p-value = 0.6842
# Se empleara la prueba de Kruskal-Wallis para comparar el peso y la recuperacion
# Prueba de hipotesis, H0: La distribucion de las muestras en todos los grupos es identica
# HA: Al menos una de las muestras en los grupos difiere en su distribucion de las demas
kruskal.test(peso ~ Recuperacion_Postcovid, data = df) # p = 0.8057, los datos respaldan la H0
##
## Kruskal-Wallis rank sum test
##
## data: peso by Recuperacion_Postcovid
## Kruskal-Wallis chi-squared = 0.060497, df = 1, p-value = 0.8057
### Evaluar el efecto de las variables predictoras (edad, peso y conteos de linfocitos B) sobre la probabilidad de norecuperacion
### Se establecera al seco como covariable para ajustar su posible efecto
# Convertir la variable de recuperación a factor
df2$Recuperacion_Postcovid <- factor(df2$Recuperacion_Postcovid, levels = c(0,1))
# Convertir la variable de sexo a factor
df2$sexo <- factor(df2$sexo)
# Ajustar el modelo de regresión logística
modelo <- glm(df2$Recuperacion_Postcovid ~ edad + peso + BhLinfocitosAbsolutos + sexo, data = df2, family = binomial)
# Resumen del modelo
summary(modelo)
##
## Call:
## glm(formula = df2$Recuperacion_Postcovid ~ edad + peso + BhLinfocitosAbsolutos +
## sexo, family = binomial, data = df2)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.28272 2.26796 -0.566 0.572
## edad 0.00223 0.02106 0.106 0.916
## peso 0.00228 0.01963 0.116 0.908
## BhLinfocitosAbsolutos -0.43754 0.45495 -0.962 0.336
## sexo2 -0.90079 0.67927 -1.326 0.185
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 77.828 on 101 degrees of freedom
## Residual deviance: 75.125 on 97 degrees of freedom
## AIC: 85.125
##
## Number of Fisher Scoring iterations: 5
# Extraer los coeficientes del modelo y sus intervalos de confianza
coeficientes <- summary(modelo)$coefficients
odds_ratios <- exp(coeficientes[,1]) # Calcular los odds ratios
ci_low <- exp(coeficientes[,1] - 1.96 * coeficientes[,2]) # Calcular el límite inferior del intervalo de confianza
ci_up <- exp(coeficientes[,1] + 1.96 * coeficientes[,2]) # Calcular el límite superior del intervalo de confianza
# Crear un data frame con los resultados
resultados <- data.frame(coeficientes = coeficientes[,1], or = odds_ratios, ci_low = ci_low, ci_up = ci_up, variable = rownames(coeficientes))
# Graficar los odds ratios con intervalos de confianza
ggplot(resultados, aes(x = variable, y = or, ymin = ci_low, ymax = ci_up)) +
geom_pointrange() +
geom_hline(yintercept = 1, linetype = "dashed", color = "red") + # Línea de referencia en OR = 1
labs(x = "",
y = "Odds Ratio",
caption = "Intervalo de Confianza al 95%") +
theme_minimal() + # Estilo del tema del gráfico
coord_cartesian(ylim = c(0, 5))
# Verifica si el directorio existe
dir.exists("02_Data")
## [1] TRUE