Modelos
Yesith Toloza
2024-01-02
Proceso de modelización de la tasa de filtracion glomerular en pacientes con AR utilizando el paquete lme4
Carge librerias y base de datos
library(readxl) library(tidyverse)
library(lme4)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.3 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.3 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## Loading required package: Matrix
##
##
## Attaching package: 'Matrix'
##
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
##
##
## New names:Curación de la data
# verificamos que la columna de fechas esté en formato de fecha y hora
df$FechaParaclinicos <- as.POSIXct(df$FechaParaclinicos)
# Ordemanos por ID y fecha
df <- df %>% arrange(Identificacion, FechaParaclinicos)
# Calculamos la diferencia de tiempo entre fechas por cada ID (en años)
df <- df %>%
group_by(Identificacion) %>%
mutate(TiempoEntreExamenes = difftime(FechaParaclinicos,
min(FechaParaclinicos),
units = "weeks")/52)
df<- subset(df, df$CKD.EPI<150) #Personas con depuracion menor a 150
df<- subset(df, df$Creatinina<6) #personas con creatininas menores a 6
df <- df[df$TiempoEntreExamenes < 8,] #personas con menos de 9 años de seguimiento
df<- df[df$N..mayor.a.10 =="Si",] # personas que tengan mas de un año de seguimiento
#quitamos los decimales
df$CKD.EPI <- as.integer(df$CKD.EPI)
#df$TiempoEntreExamenes <- as.integer(df$TiempoEntreExamenes)
#Seleccionamos la Data que queremos
df1<- data.frame("ID" = df$Identificacion,
"Time" = df$TiempoEntreExamenes,
"CKDEPI" = df$CKD.EPI,
"HTA_Inicio" = df$V34.Hta.Al.Diagnóstico,
"HTA_Actual" = df$X59..Hta.Actual,
"DM_Inicio" = df$V35.Dm.Al.Diagnótico,
"DM_Actual" = df$V60.Dm.Actual,
"IMC" = df$IMC.Actual..Kg.M2.,
"Sex" = df$Sexo,
"cod" = 1
)
df1$nsex <- as.numeric(df1$Sex=="FEMENINO") # 1 = Femenino, 0 = Masculino
#Curacion los datos para los analisis de regresion
df1 <- df1 %>%
arrange(ID, Time, HTA_Inicio, HTA_Actual,
DM_Inicio, DM_Actual,IMC, nsex )%>%
group_by(ID, Time, HTA_Inicio, HTA_Actual,
DM_Inicio, DM_Actual,IMC, nsex ) %>%
summarise("CKDEPImean" = mean(CKDEPI, na.rm = T))## `summarise()` has grouped output by 'ID', 'Time', 'HTA_Inicio', 'HTA_Actual',
## 'DM_Inicio', 'DM_Actual', 'IMC'. You can override using the `.groups` argument.df1$NuevaHTA_DM <- c(df1$DM_Actual+df1$HTA_Actual)
df1$NuevaHTA_DMcod<- factor(df1$NuevaHTA_DM,
levels = c("0", "1", "2"),
labels = c("Ninguna",
"DM o HTA",
"DM y HTA"))
df1$NuevaHTA_DMcod <- relevel(df1$NuevaHTA_DMcod, ref = "Ninguna")
#Datos para analisis descriptivo
df2<- subset(df, Duplicado==FALSE)Proceso de modelización
Ajuste del modelo lineal mixto (lmer) Se ajusta un modelo
lineal mixto usando la función lmer del paquete
lme4. Este modelo está modelando la variable de respuesta
CKDEPI.
El código es el siguiente:
#codigo de la regresion
fm1 <- lmer(CKDEPImean ~ Time + NuevaHTA_DMcod + (1|ID),
data = df1, control=lmerControl())
iqrvec <- sapply(simulate(fm1, 1000), IQR)
obsval <- IQR(df1$CKDEPImean, na.rm = T)
df1$fitted_values <- fitted(fm1)
vv <- vcov.merMod(fm1, corr=TRUE)
correlacion <- as(vv, "corMatrix")
cor_matrix <- as.matrix(correlacion)
# Nombres para las filas y columnas de la matrix
nombres_variables <- colnames(cor_matrix)
rownames(cor_matrix) <- c("Intercept", "Time",
"DM o HTA", "DM y HTA")
colnames(cor_matrix) <- c("Intercept", "Time",
"DM o HTA", "DM y HTA")Estadisticas del modelo
el modelo lineal mixto ajustado utiliza el método REML (Máxima Verosimilitud Restringida) para analizar cómo la variable dependiente CKDEPImean (tasa de filtración glomerular estimada) se relaciona con el tiempo (Time) y la condición de salud (NuevaHTA_DMcod) de los individuos, con una variación aleatoria en los interceptos para diferentes individuos (ID). A continuación, se detallan las conclusiones y resultados basados en los resultados del modelo:
Resultados de los efectos fijos
Intercepto (101.49854): representa el valor estimado de CKDEPI cuando Time es 0 y para la categoría de referencia de NuevaHTA_DMcod.
Time (-0.97866): Indica que, en promedio, por cada unidad de incremento en el tiempo, CKDEPI disminuye en 0.97866 unidades. El alto valor absoluto de “t” (-38.350) sugiere que este efecto es estadísticamente significativo.
NuevaHTA_DMcod (DM o HTA y DM y HTA): Las categorías “DM o HTA” y “DM y HTA” están asociadas con disminuciones en CKDEPI de 6.43545 y 7.54264 unidades respectivamente, en comparación con personas sin dichas patologias. Ambas categorías muestran efectos significativos en CKDEPI.
Resultados de los efectos aleatorios
Intercepto Aleatorio (Varianza: 188.75): Indica una variabilidad considerable en el intercepto de CKDEPI entre diferentes individuos (ID). La desviación estándar (13.739) muestra cuánto varían los valores base de CKDEPI de un individuo a otro.
Residual (Varianza: 57.86): Representa la variabilidad en CKDEPI que no es explicada por el modelo. La desviación estándar (7.607) da una idea de la dispersión de los residuos alrededor de los valores predichos.
Residuos Los residuos tienen un rango bastante amplio, lo que podría indicar la presencia de valores atípicos o una posible no normalidad en los residuos, la valoracion de los resuduos se mira mas adelante.
Evaluación del modelo El resultado obtenido en la función
mean(obsval >= c(obsval, iqrvec)) indica la proporción
de casos en los que la diferencia entre los IQR observados y los IQR
simulados es mayor o igual a la diferencia entre los IQR observados y
los IQR simulados.
En este caso específico, un valor de 0.7462 significa que aproximadamente el 74% de las simulaciones generadas muestran una diferencia en el IQR igual o mayor a la diferencia entre los IQR observados y los IQR simulados. Esto podría interpretarse como una evaluación del ajuste del modelo a los datos observados: un valor alto sugiere que el modelo simulado produce resultados bastante similares en términos de la variabilidad (medida por el IQR) en comparación con los datos observados, lo que indica un buen ajuste del modelo.
#Resumen del modelo.
coef_df <- as.data.frame(summary(fm1)$coefficients)
coef_df$Term <- rownames(coef_df)
summary(fm1)## Linear mixed model fit by REML ['lmerMod']
## Formula: CKDEPImean ~ Time + NuevaHTA_DMcod + (1 | ID)
## Data: df1
##
## REML criterion at convergence: 186911
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -9.5435 -0.4639 0.0681 0.5292 9.1994
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 188.75 13.739
## Residual 57.86 7.607
## Number of obs: 26077, groups: ID, 1970
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 101.49854 0.36160 280.690
## Time -0.97866 0.02552 -38.350
## NuevaHTA_DMcodDM o HTA -6.43545 0.79489 -8.096
## NuevaHTA_DMcodDM y HTA -7.54264 2.06772 -3.648
##
## Correlation of Fixed Effects:
## (Intr) Time NHTAoH
## Time -0.157
## NHTA_DMDMoH -0.441 -0.014
## NHTA_DMDMyH -0.169 -0.009 0.078mean(obsval >= c(obsval, iqrvec))## [1] 0.7452547Gráficos
# Grafica datos originales y las líneas ajustadas
ggplot(df1, aes(x = Time, y = CKDEPImean, group = ID)) +
geom_line(aes(y = fitted_values), color = "red") + # Líneas ajustadas
geom_point(color = "black", alpha = 0.3) + # Datos originales
labs(x = "Tiempo en Años", y = "CKD-EPI") +
ggtitle("Líneas estimadas por el modelo y puntos utilizados para el modelo")## Don't know how to automatically pick scale for object of type <difftime>.
## Defaulting to continuous.
# Grafica matriz de correlación con información de significancia estadística
corrplot::corrplot(cor_matrix,
p.mat = cor_matrix,
type="lower",
tl.col="black",
tl.srt = 90,
pch.col = "black",
insig = "p-value",
sig.level = -1,
col = heat.colors(11))
#Grafica de Residuales
plot(fm1, type = c("p", "smooth"))
plot(fm1, sqrt(abs(resid(.))) ~ fitted(.), type = c("p", "smooth"))
# Utilizar ggplot para graficar los efectos fijos
ggplot(coef_df, aes(x = Term, y = Estimate,
ymin = Estimate - 2*`Std. Error`,
ymax = Estimate + 2*`Std. Error`)) +
geom_pointrange() +
theme_minimal() +
labs(title = "Efectos Fijos del Modelo", x = "Términos", y = "Estimaciones") +
geom_hline(yintercept = 0, linetype = "dashed", color = "red")
library(lattice)## Warning: package 'lattice' was built under R version 4.3.2qqmath(fm1)
#xyplot(CKDEPImean ~ Time + NuevaHTA_DMcod + (1|ID), df1,
# aspect = "xy",
# #layout = c(9, 2),
# type = c("g", "p", "r"),
# index.cond = function(x, y) coef(lm(y ~ x))[2],
# xlab = "Año de seguimiento",
# ylab = "CKD-EPI",
#as.table = TRUE,
# na.action = na.exclude
# )
pf <- profile(fm1)
xyplot(pf)
densityplot(pf)
splom(pf)
Estadistica Descriptiva
#promedio de edad de los participantes
mean(df2$Edad, na.rm=T)## [1] 47.43928quantile(df2$Edad, c(0.25, 0.75))## 25% 75%
## 39 56conteo <- df1 %>%
group_by(ID) %>%
summarize(numero_examenes = n(),
"min" = min(Time),
"max"= max(Time),
"HTA" = min(HTA_Actual),
"DM" = min(DM_Actual))
# promedio de veces que las personas se calculo la TFG
promedio_examenes <- mean(conteo$numero_examenes)
# Calculo de percentiles 25, 50 (mediana) y 75
percentiles <- quantile(conteo$numero_examenes,
probs = c(0.25, 0.5, 0.75))
print(paste("Promedio de exámenes por persona:", promedio_examenes))## [1] "Promedio de exámenes por persona: 13.2370558375635"print(paste("Percentil 25:", percentiles[1]))## [1] "Percentil 25: 7"print(paste("Mediana:", percentiles[2]))## [1] "Mediana: 13"print(paste("Percentil 75:", percentiles[3]))## [1] "Percentil 75: 19"print(paste("N° de personas con HTA es ", table(conteo$HTA)[2]))## [1] "N° de personas con HTA es 396"print(paste("N° de personas con DM es ", table(conteo$DM)[2]))## [1] "N° de personas con DM es 84"sessionInfo()## R version 4.3.1 (2023-06-16 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 11 x64 (build 22621)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=Spanish_Colombia.utf8 LC_CTYPE=Spanish_Colombia.utf8
## [3] LC_MONETARY=Spanish_Colombia.utf8 LC_NUMERIC=C
## [5] LC_TIME=Spanish_Colombia.utf8
##
## time zone: America/Bogota
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] lattice_0.22-5 lme4_1.1-34 Matrix_1.6-1 lubridate_1.9.2
## [5] forcats_1.0.0 stringr_1.5.0 dplyr_1.1.3 purrr_1.0.2
## [9] readr_2.1.4 tidyr_1.3.0 tibble_3.2.1 ggplot2_3.4.3
## [13] tidyverse_2.0.0 readxl_1.4.3
##
## loaded via a namespace (and not attached):
## [1] sass_0.4.7 utf8_1.2.3 generics_0.1.3 stringi_1.7.12
## [5] hms_1.1.3 digest_0.6.33 magrittr_2.0.3 evaluate_0.21
## [9] grid_4.3.1 timechange_0.2.0 fastmap_1.1.1 cellranger_1.1.0
## [13] jsonlite_1.8.7 fansi_1.0.4 scales_1.2.1 jquerylib_0.1.4
## [17] cli_3.6.1 rlang_1.1.1 munsell_0.5.0 splines_4.3.1
## [21] withr_2.5.0 cachem_1.0.8 yaml_2.3.7 corrplot_0.92
## [25] tools_4.3.1 tzdb_0.4.0 nloptr_2.0.3 minqa_1.2.6
## [29] colorspace_2.1-0 boot_1.3-28.1 vctrs_0.6.3 R6_2.5.1
## [33] lifecycle_1.0.3 MASS_7.3-60 pkgconfig_2.0.3 pillar_1.9.0
## [37] bslib_0.5.1 gtable_0.3.3 Rcpp_1.0.11 glue_1.6.2
## [41] highr_0.10 xfun_0.40 tidyselect_1.2.0 rstudioapi_0.15.0
## [45] knitr_1.43 farver_2.1.1 nlme_3.1-164 htmltools_0.5.6
## [49] labeling_0.4.2 rmarkdown_2.24 compiler_4.3.1