Modelacion
Yesith Toloza
2024-03-14
Proceso de modelización de la EDSS en pacientes con NMO 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$Fechaclinimetria <- as.POSIXct(df$Fecha)
# Ordemanos por ID y fecha
df <- df %>% arrange(Identificacion, Fechaclinimetria)
# Calculamos la diferencia de tiempo entre fechas por cada ID (en años)
library(lubridate)
df <- df%>%
filter(year < 2024) %>%
group_by(Identificacion) %>%
mutate(TiempoEntreExamenes = interval(start = min(Fechaclinimetria),
end = Fechaclinimetria) / months(1))
#Seleccionamos la Data que queremos
df1<- data.frame("ID" = df$Identificacion,
"Time" = df$TiempoEntreExamenes,
"EDSS" = df$Resultado,
"year"= df$year ,
"Edad" = df$Edad.Diagnostico.23,
"Sex" = df$Sexo.23,
"cod" = 1
)
df1$Sex<-factor(df1$Sex)
#df1$Medicamento<-factor(df1$Medicamento)
#df1$Medicamento<-relevel(df1$Medicamento, ref="SIN TRATAMIENTO")
df1$Sex<-relevel(df1$Sex, ref="MASCULINO")
#Curacion los datos para los analisis de regresionProceso 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
EDSS
El codigo es el siguiente:
#codigo de la regresion
fm1 <- lmer(EDSS ~ Time + Edad+Sex+(1|ID)+(1|year),
data = df1, control=lmerControl())
iqrvec <- sapply(simulate(fm1, 1000), IQR)
obsval <- IQR(df1$EDSS, 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", "Edad", "FEMENINO")
colnames(cor_matrix) <- c("Intercept", "Time", "Edad","FEMENINO")Interpretacion del codigo anterior
Modelo lineal mixto La primera líinea esta ajustando el
modelo lineal mixto con la función lmer() del paquete
lme4. El modelo está estimando la relación entre la Escala
de Estado de Discapacidad Expandida (EDSS) y el tiempo (Time), la edad
al diagnostico (Edad), el sexo (Sex), además tiene en cuenta los efectos
aleatorios para cada ID y el año del seguimiento.
Simulación e IQR Luego, se simula el modelo fm1
1000 veces y se calcula el rango intercuartílico IQR para
cada simulación. Esto se hace con el fin de evaluar la estabilidad de la
estimación del IQR del modelo. la fincion obsval es el IQR
calculado de la variable EDSS en los datos observados (df1).
Valores ajustados Se añaden al dataframe df1
los valores ajustados por el modelo fitted(fm1), que son
las predicciones del modelo para cada observación.
Matriz de covarianza y correlación por ultimo, se obtiene la
matriz de covarianza para el modelo
vcov.merMod(fm1, corr=TRUE), y se convierte en una matriz
de correlación corMatrix. Esta matriz de correlación
proporciona información sobre la relación entre los parámetros estimados
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: EDSS ~ Time + Edad + Sex + (1 | ID) + (1 | year)
## Data: df1
##
## REML criterion at convergence: 1564.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.2411 -0.4194 0.0102 0.3943 4.9574
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 5.4450 2.3334
## year (Intercept) 0.0127 0.1127
## Residual 0.7382 0.8592
## Number of obs: 502, groups: ID, 73; year, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.776905 1.284668 2.162
## Time 0.011399 0.003692 3.087
## Edad 0.060860 0.021898 2.779
## SexFEMENINO -2.507390 1.017419 -2.464
##
## Correlation of Fixed Effects:
## (Intr) Time Edad
## Time -0.066
## Edad -0.652 -0.003
## SexFEMENINO -0.646 0.013 -0.112mean(obsval >= c(obsval, iqrvec))## [1] 0.999001Interpretacion del codigo anterior
Extraccion de coeficientes crea un dataframe llamado
coef_df a partir de los coeficientes del resumen del modelo
(summary(fm1)$coefficients) y luego añade una nueva columna
llamada Term que contiene los nombres de las variables o
términos del modelo.
Resumen del modelo La función summary(fm1)
proporciona un resumen detallado del modelo lineal mixto ajustado. El
resumen incluye el criterio REML en la convergencia, los residuos
escalados, los efectos aleatorios y fijos, así como la correlación entre
los efectos fijos.
Interpretación del resumen del modelo
Efectos aleatorios Hay dos efectos aleatorios en el modelo, ID y year, que tienen varianzas de 5.4 y 0.01, respectivamente. Esto indica que hay variabilidad considerable en la puntuación de EDSS entre diferentes sujetos, pero muy poca variabilidad entre años.
Residuos La variabilidad residual (variabilidad no explicada por el modelo) es de 0.7382, con una desviación estándar de 0.8592.
Numero de observaciones El modelo se ajustó con 502 observaciones, que implican 73 personas (ID) y 5 años diferentes (year).
Efectos fijos El tiempo entre las mediciones de las clinimetrias (Time) y la edad al diagnostico (Edad) están positivamente asociados con la puntuación de EDSS, y el sexo femenino (SexFEMENINO) esta asociado negativamente con la puntuación de EDSS, lo que significa que, en promedio, las mujeres tienen una puntuación más baja que los hombres, controlando por las otras variables.
Graficos
# Grafica datos originales y las líneas ajustadas
ggplot(df1, aes(x = Time, y = EDSS, group = ID)) +
geom_line(aes(y = fitted_values), color = "red") + # Líneas ajustadas
geom_point(color = "black", alpha = 0.3) + # Datos originales
labs(x = "Time in months", y = "EDSS") #+
#ggtitle("Líneas estimadas por el modelo y puntos utilizados para el modelo")
# Grafica matriz de correlación
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 General
plot(fm1, type = c("p", "smooth"))
plot(fm1, sqrt(abs(resid(.))) ~ fitted(.), type = c("p", "smooth"))
# grafica 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.2# tamaño de la fuente y los puntos sean adecuados para la visualización
trellis.par.set(theme = col.whitebg())
trellis.par.set(cex = 0.7, pch = 16)
# gráfico
xyplot(EDSS ~ Time + Edad + Sex + (1 | ID), df1,
index.cond = function(x, y)
coef(lm(y ~ x))[2],
xlab = "Mes de seguimiento",
type = c("g", "p", "r"),
ylab = "EDSS",
as.table = TRUE,
na.action = na.exclude,
panel = function(x, y) {
panel.xyplot(x, y)
panel.abline(h = median(y, na.rm = TRUE), col = "red", lwd = 2) # Línea mediana
}
)
#Resuduales
pf <- profile(fm1)
xyplot(pf)
densityplot(pf)
qqmath(fm1)
splom(pf)
Conclusiones El modelo sugiere que tanto el tiempo como la edad tienen un efecto positivo en la puntuación de EDSS, mientras que ser de sexo femenino está asociado con una puntuación más baja. Además, hay una variabilidad significativa asociada con los sujetos (ID). La alta proporción del IQR observado frente a los simulados sugiere que el modelo puede no estar capturando completamente la variabilidad de los datos.
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