Informe CALUCI
INTRODUCCIÓN:
En este informe vamos a analizar los datos del proyecto CALUCI que registran el estado de salud de pacientes en UCI en el Hospital de Córdoba mediante wearable technologies en forma de variables contínuas junto con variables de ambiente también recogidas de forma contínua.
El objetivo de este estudio estadístico es doble: (1) en primer lugar, examinar cómo influyen las variables ambientales sobre las variables fisiológicas, indicando con detalle qué variables son y cómo influyen sobre la salud de los pacientes en UCI; (2) en segundo lugar, debido a la naturaleza contínua de los datos, nos disponemos a comprobar si la modelización de los datos se ve beneficiada por el uso de técnicas de datos funcionales.
Para realizar este estudio vamos a:
- Purificar los datos brutos de que se disponen.
- Filtrar las variables de interés y gestionar los datos faltantes y datos atípicos.
- Examinar la correlación entre variables para conocer las interacciones más prometedoras.
- Formatear los datos en forma de Dato Funcional.
- Examinar de la correlación funcional y comprobar si las indicaciones coinciden con las del paso 3.
- Proponer una serie de modelos tentativos con diferentes metodologías.
- Realizar combinatorias de variables para proponer modelos de regresión aditiva funcional.
- Comparar dichos modelos con otros posibles: modelos lineares, aditivos generalizados, y lineales funcionales.
- Evaluar dichos modelos mediante una tabla que registre los diferentes R-Sqared y AIC.
- Utilizar dichos modelos para realizar predicciones a corto plazo sobre el estado de salud del paciente.
- Evaluar gráficamente estas predicciones mediante Plots de Violines y BoxPlots.
- Elaborar unas conclusiones sucintas.
Procedemos entonces a la carga y limpieza de datos. El avance del proceso de análisis se puede seguir en la barra desplegable de la izquierda.
PARTE 1: PREPROCESADO
|- Preámbulo:
#* Set working directory: ----
setwd("~/GitHub/FDA-CALUCI")
#* Tune Options: ----
options(scipen = 6, digits = 4) # View outputs in non-scientific notation
memory.limit(30000000) # This is needed on some PCs to increase memory allowance
#* Load up packages: ----
library(fda.usc); library(readxl); library(DataExplorer); library(Hmisc);
library(data.table); library(plyr); library(zoo); library(mgcv); library(ggpubr);
library(ellipse); library(biostatUZH); library(RColorBrewer); library(ggplot2);
library(dplyr); library(forcats)
# install.packages("biostatUZH", repos = "http://R-Forge.R-project.org")
# install.packages("rmdformats")Labores de limpieza y curado de datos. Se trata de una base de datos compleja por lo que la limpieza va a ser muy costosa. Algunas partes de esta limpieza ya vienen hechas en Excel debido a la irregularidad de los datos.
|- Carga de Datos:
Parte de este proceso está facilitado porque se ha trabajado directamente en Excel, fusionando variables como la fecha para que sean más fácilmente manejables en R. También se ha creado la variable ID para cada una de las entradas y se ha eliminado el marcado antiguo que en este software es farragoso y poco práctico.
Datos_Totales <- read_excel("Datos Totales.xlsx",
col_types = c("text", "date", "skip",
"skip", "numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "text", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric"))
save(Datos_Totales, file = "datos_totales.RData")load("~/GitHub/FDA-CALUCI/datos_totales.RData")
str(Datos_Totales[ , 1:48])## tibble [397,623 x 48] (S3: tbl_df/tbl/data.frame)
## $ ID : chr [1:397623] "A1" "A1" "A1" "A1" ...
## $ Fecha y hora : POSIXct[1:397623], format: "2020-10-11 09:08:30" "2020-10-11 09:09:00" ...
## $ Marca : num [1:397623] 0 0 0 0 0 0 0 0 0 0 ...
## $ Tª Amb. : num [1:397623] 24.3 24.5 24.7 24.7 24.8 24.8 24.8 24.8 24.8 24.8 ...
## $ Luz Vis. 1 : num [1:397623] 52.3 73.6 71.8 54.4 52.3 ...
## $ Luz Vis. 2 : num [1:397623] 15.8 21.9 21.1 15.9 15.4 ...
## $ ° X : num [1:397623] -12.3 9.9 10 9.8 9.6 10 9.6 10 9.7 9.8 ...
## $ ° Y : num [1:397623] 75.4 6.6 6.6 6.7 6.7 6.5 6.3 6.6 6.9 6.6 ...
## $ ° Z : num [1:397623] 7.6 78 77.8 78 78.2 77.9 78.4 77.8 78 78.1 ...
## $ Suma ° : num [1:397623] 2798 4200 131 143 132 ...
## $ Suma G : num [1:397623] 13.7 19.5 11.2 11.2 11.3 ...
## $ Nº Mov. : num [1:397623] 50 137 0 0 0 0 0 0 0 0 ...
## $ Tª Sonda : num [1:397623] 23.5 23.8 25.1 26.2 27 27.7 28.2 28.7 29.2 29.6 ...
## $ %HR : num [1:397623] 44.8 45.5 47.5 47.4 47.1 46.8 46.5 46.1 45.8 45.5 ...
## $ Luz IR 1 : num [1:397623] 10.99 12.31 11.43 8.78 8.35 ...
## $ Pico Luz : num [1:397623] 95.9 113.1 114.3 54.9 54.8 ...
## $ Prom. dB : num [1:397623] 65.7 65.7 66.7 64.4 65.3 64.8 69.5 64.2 64.6 64.4 ...
## $ FC : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ TAS : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ TAD : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ SpO2 : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ FR : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ GCSm : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ RASS : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ Adrenalina : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ Noradrenalina : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ Dobutamina : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ Milrirona : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ Levosimendan : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ Midazolam : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ Propofol : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ Fentanilo : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ Remifentanilo : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ Dexmedetomidina: num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ Cl. Morfico : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ SALIVA : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ HDPWC : chr [1:397623] NA NA NA NA ...
## $ Modo VMI : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ ECMO AV : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ Cisatracurio : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ Amiodarona : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ GN : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ Sillon : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ HDI : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ NO3 : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ BCIA : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ VMK : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ ONAF Traqueo : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...|- Cribado de Variables:
Una vez cargados los datos, vamos a realizar un cribado de variables para quedarnos únicamente con aquellas que consideramos útiles para este informe. Estas van a ser aquellas que podrían considerarse como ‘meta-datos’ (ID y fecha/hora), junto con las variables fisiológicas medidas en el paciente y las variables ambientales que son candidatas a afectar a dichas medidas fisiológicas.
sub_datos <- Datos_Totales[, c("ID", "Fecha y hora", "FC", "TAS", "TAD", "SpO2", "FR", "GCSm",
"RASS", "Tª Amb.", "Luz Vis. 1", "Luz Vis. 2", "Luz IR 1",
"Pico Luz", "Prom. dB")]El resultado es un data.frame con la siguiente estructura:
head(sub_datos)## # A tibble: 6 x 15
## ID `Fecha y hora` FC TAS TAD SpO2 FR GCSm RASS `Tª Amb.`
## <chr> <dttm> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 A1 2020-10-11 09:08:30 NA NA NA NA NA NA NA 24.3
## 2 A1 2020-10-11 09:09:00 NA NA NA NA NA NA NA 24.5
## 3 A1 2020-10-11 09:09:30 NA NA NA NA NA NA NA 24.7
## 4 A1 2020-10-11 09:10:00 NA NA NA NA NA NA NA 24.7
## 5 A1 2020-10-11 09:10:30 NA NA NA NA NA NA NA 24.8
## 6 A1 2020-10-11 09:11:00 NA NA NA NA NA NA NA 24.8
## # ... with 5 more variables: `Luz Vis. 1` <dbl>, `Luz Vis. 2` <dbl>,
## # `Luz IR 1` <dbl>, `Pico Luz` <dbl>, `Prom. dB` <dbl>str(sub_datos)## tibble [397,623 x 15] (S3: tbl_df/tbl/data.frame)
## $ ID : chr [1:397623] "A1" "A1" "A1" "A1" ...
## $ Fecha y hora: POSIXct[1:397623], format: "2020-10-11 09:08:30" "2020-10-11 09:09:00" ...
## $ FC : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ TAS : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ TAD : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ SpO2 : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ FR : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ GCSm : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ RASS : num [1:397623] NA NA NA NA NA NA NA NA NA NA ...
## $ Tª Amb. : num [1:397623] 24.3 24.5 24.7 24.7 24.8 24.8 24.8 24.8 24.8 24.8 ...
## $ Luz Vis. 1 : num [1:397623] 52.3 73.6 71.8 54.4 52.3 ...
## $ Luz Vis. 2 : num [1:397623] 15.8 21.9 21.1 15.9 15.4 ...
## $ Luz IR 1 : num [1:397623] 10.99 12.31 11.43 8.78 8.35 ...
## $ Pico Luz : num [1:397623] 95.9 113.1 114.3 54.9 54.8 ...
## $ Prom. dB : num [1:397623] 65.7 65.7 66.7 64.4 65.3 64.8 69.5 64.2 64.6 64.4 ...En este data.frame vemos que los datos fisiológicos de los pacientes son tomados cada 15 minutos en contraste con las variables ambientales que se miden cada medio minuto. Esta estructura va a requerir un procesado más complejo, pero también abre la puerta a pensar en usar técnicas de datos funcionales para el análisis estadístico, ya que las mediciones ambientales son muy frecuentes y con pocas oscilaciones, por lo que podrían ser consideradas como instancias de una función subyacente y resultar en un beneficio a la hora de modelizar y predecir.
PARTE 2: ESTADÍSTICOS DESCRIPTIVOS
Lo primero que llevaremos a cabo es un análisis estadístico descriptivo detallado para conocer en profundidad el dataset y ver qué posibles modificaciones serían necesarias, por ejemplo, en cuanto a la presencia de datos faltantes, atípicos…
|- Datos por paciente:
Los datos ambientales y fisiológicos están asociados a un paciente (variable ID) y, debido a las necesidades de la toma de datos y su casuística, el número de entradas de datos registrados para cada paciente no es idéntico para todos los casos, como sería lo ideal.
ID <- as.vector(t(unique(sub_datos[,1]))) # Contar los IDs existentes
Ncasos <- data.frame(ID = integer(), number = numeric()) # Creamos tabla
for (i in 1:length(ID)) {
temp <- data.frame(ID = as.character(ID[i]),
number = as.numeric(length(sub_datos$ID[sub_datos$ID == ID[i]])))
Ncasos <- rbind(Ncasos, temp)} # Contamos el Nº de casos por paciente y lo mandamos a la tabla
rm(i); rm(temp) # Por higiene de códigoVemos que tenemos un total de 77 pacientes en este estudio, los cuales a su vez tienen un número muy variado de entradas de datos. Con dos líneas extra de código vamos a ver qué casos son los más extremos a la baja (pudiendo ser necesario eliminarlos) y un histograma para tener una imagen de la distribución en cuestión.
sort(Ncasos$number)[1:5] # Pacientes con el menor número de registros## [1] 480 828 1410 1585 2702hist(Ncasos$number, main = "Histograma del Nº de Entradas del Paciente",
xlab = "Nº de Entradas",
ylab = "Frecuencia") # Histograma del número de registros en cada paciente
El histograma nos indica que los pacientes tienen en torno a 5000-6000 entradas cada uno. No obstante, hay varios pacientes que oscilan entre las 2500-4500, y un reducido número de pacientes que, por alguna razón, no alcanzan ni los 2000. Los pacientes C9 con 480, B9 con 828, D17 con 1410, D11 con 1585 casos son los ejemplos más extremos.
Se ha barajado la idea de eliminarlos del análisis, no obstante, tras una reunión con las personas responsables de la toma de datos, se decide conservar estos datos como válidos mientras no sean datos atípicos o ilógicos.
|- Casos Atípicos y Casos Imposibles:
Continuando con el análisis descriptivo, procedemos a marcar qué datos ambientales son directamente ilógicos o imposibles siguiendo recomendaciones de los colaboradores clínicos de este proyecto. Las etiquetas serán las siguientes:
Labels:
Outlier en Tª Amb.
Outlier en Prom. Db.
Outlier en Luz 1
Outlier en Luz 2
Outlier en Luz IR
Outlier en Pico de Luz
Que se corresponden con los límites pre-definidos. Por ejemplo: considerar un caso atípico un Promedio de Decibelios inferior a 0 (por ser imposible) o superior a 120, que sería insostenible en ninguna unidad de cuidados intensivos.
sub_datos$label <- rep(0, nrow(sub_datos))
# Esta variable 'label' sirve para marcar si hay atipicos en esta fila
colnames(sub_datos) <- c("ID", "Fecha y hora", "FC", "TAS", "TAD", "SpO2",
"FR", "GCSm", "RASS", "T Amb.", "Luz Vis. 1",
"Luz Vis. 2", "Luz IR 1", "Pico Luz", "Prom. dB",
"Label")
# Cambiar Tª for T para prevenir erroresfor (i in 1:nrow(sub_datos)) {
if (sub_datos$`T Amb.`[i] < 24 | sub_datos$`T Amb.`[i] > 26) {
sub_datos$Label[i] <- 1 # T Amb.
} else if (sub_datos$`Prom. dB`[i] < 0 | sub_datos$`Prom. dB`[i] > 120) {
sub_datos$Label[i] <- 2 # Prom. Deb.
} else if (sub_datos$`Luz Vis. 1`[i] < 0 | sub_datos$`Luz Vis. 1`[i] > 300) {
sub_datos$Label[i] <- 3 # Luz 1
} else if (sub_datos$`Luz Vis. 2`[i] < 0 | sub_datos$`Luz Vis. 2`[i] > 300) {
sub_datos$Label[i] <- 4 # Luz 2
} else if (sub_datos$`Luz IR 1`[i] < 0 | sub_datos$`Luz IR 1`[i] > 300) {
sub_datos$Label[i] <- 5 # Luz IR
} else if (sub_datos$`Pico Luz`[i] < 0 | sub_datos$`Pico Luz`[i] > 300) {
sub_datos$Label[i] <- 6 # Pico Luz
} else {
#print(i)
next
}
}lable.table <- table(sub_datos$Label)
save(lable.table, file = "lable_table.RData")load("lable_table.RData")
lable.table##
## 0 1 3 4 5 6
## 120138 276028 1315 1 6 135Vemos que una gran cantidad de las entradas para temperatura ambiental son atípicas, por lo que Tª Amb deberá ser usada con cautela como variable explicativa en los modelos que hagamos de aquí en adelante, incluso después de la eliminación de estos casos atípicos, todo apunta a errores sistemáticos en la toma de estos datos.
En menor medida, también Luz 1 es una variable con bastantes casos atípicos, así como la variable Pico Luz. El resto, en principio, parecen contener valores dentro de lo esperable y lógico.
Así como para las variables ambientales se han utilizado los rangos proporcionados por el equipo clínico, para otras variables como GCSm o RASS, se han considerado como valores imposibles aquellos que, directamente, caen fuera de las propias posibilidades descritas para el test.
Es de recalcar que Prom.dB es la única variable ambiental que no tiene ningún dato atípico. Este dato contrasta con las medidas de temperatura.
Eliminación:
Los valores fuera de estos rangos son substituidos por NA (dato faltante). También se convierten en NAs algunos valores de variables fisiológicas que se saldrían de lo normal, aunque se mantiene un margen bastante amplio.
library(naniar); library(magrittr)
datos_logicos <- sub_datos %>%
replace_with_na_at(
.vars = "FC",
condition = ~(.x < 10 | .x > 200)
) %>%
replace_with_na_at(
.vars = "TAS",
condition = ~(.x < 20 | .x > 300)
) %>%
replace_with_na_at(
.vars = "TAD",
condition = ~(.x < 20 | .x > 200)
) %>%
replace_with_na_at(
.vars = "SpO2",
condition = ~(.x < 40 | .x > 100)
) %>%
replace_with_na_at(
.vars = "FR",
condition = ~(.x < 5 | .x > 100)
) %>%
replace_with_na_at(
.vars = "GCSm",
condition = ~(.x < 1 | .x > 6)
) %>%
replace_with_na_at(
.vars = "RASS",
condition = ~(.x < -5 | .x > 6)
) %>%
replace_with_na_at(
.vars = "T Amb.",
condition = ~(.x < 24 | .x > 26)
) %>%
replace_with_na_at(
.vars = 'Luz Vis. 1',
condition = ~(.x > 300)
) %>%
replace_with_na_at(
.vars = 'Luz Vis. 2',
condition = ~(.x > 300)
) %>%
replace_with_na_at(
.vars = 'Luz IR 1',
condition = ~(.x > 300)
) %>%
replace_with_na_at(
.vars = 'Pico Luz',
condition = ~(.x > 300)
)
save(datos_logicos, file = "datos_logicos.RData")load("datos_logicos.RData")|- Datos Faltantes:
En conjunto, una vez analizados los datos atípicos que hemos convertido en datos faltantes, más los datos faltantes per se en el dataset, obtenemos los siguientes estadísticos:
Ambientales:
Para las variables ambientales explicativas vemos unos porcentajes de datos faltantes poco preocupante salvo para la variable Tª Ambiental, que se acerca a un 70% de datos faltantes a lo que hay que sumar que ya anteriormente recomendábamos no utilizar esta variable por la cantidad de datos atípicos que había mostrado.
plot_missing(datos_logicos[ , 10:15],
theme_config = list(legend.position = c("none")))
Fisiológicos:
En cuanto a las variables respuesta fisiológicas, ninguna de ellas muestra una cantidad de datos faltantes anormal o que pueda suponer un problema para el análisis estadístico.
load("explicativas.RData")
plot_missing(df_luz1[, 3:9],
theme_config = list(legend.position = c("none")))
|- Descriptivos Generales:
En esta sección, junto con la siguiente de Análisis de Frecuencias, vamos a valorar la naturaleza de las variables explicativas y respuesta una vez los datos atípicos han sido regularizados y los faltantes analizados. Adjuntamos a continuación una serie de estadísticos obtenidos de forma sistemática para todas las variables.
summary(datos_logicos[ , 3:15], digits = 2)## FC TAS TAD SpO2
## Min. : 19 Min. : 37 Min. : 20 Min. : 40
## 1st Qu.: 78 1st Qu.:105 1st Qu.: 55 1st Qu.: 97
## Median : 92 Median :118 Median : 62 Median : 99
## Mean : 91 Mean :119 Mean : 64 Mean : 98
## 3rd Qu.:104 3rd Qu.:131 3rd Qu.: 70 3rd Qu.:100
## Max. :197 Max. :238 Max. :168 Max. :100
## NA's :384864 NA's :384871 NA's :384868 NA's :384868
## FR GCSm RASS T Amb.
## Min. : 6 Min. :1 Min. :-5 Min. :24
## 1st Qu.: 16 1st Qu.:1 1st Qu.:-5 1st Qu.:24
## Median : 20 Median :1 Median :-3 Median :25
## Mean : 20 Mean :3 Mean :-3 Mean :25
## 3rd Qu.: 24 3rd Qu.:5 3rd Qu.: 0 3rd Qu.:25
## Max. :100 Max. :6 Max. : 2 Max. :26
## NA's :384861 NA's :391924 NA's :387890 NA's :276028
## Luz Vis. 1 Luz Vis. 2 Luz IR 1 Pico Luz Prom. dB
## Min. : 0 Min. : 0.0 Min. : 0.0 Min. : 0 Min. : 25
## 1st Qu.: 1 1st Qu.: 0.3 1st Qu.: 0.0 1st Qu.: 1 1st Qu.: 57
## Median : 11 Median : 2.7 Median : 0.3 Median : 11 Median : 64
## Mean : 33 Mean : 9.4 Mean : 3.8 Mean : 34 Mean : 63
## 3rd Qu.: 44 3rd Qu.: 11.8 3rd Qu.: 1.8 3rd Qu.: 46 3rd Qu.: 66
## Max. :300 Max. :300.0 Max. :300.0 Max. :300 Max. :110
## NA's :3873 NA's :1599 NA's :1831 NA's :4505Hmisc::describe(datos_logicos[ , 3:15])## datos_logicos[, 3:15]
##
## 13 Variables 397623 Observations
## --------------------------------------------------------------------------------
## FC
## n missing distinct Info Mean Gmd .05 .10
## 12759 384864 124 1 90.84 21.68 59 64
## .25 .50 .75 .90 .95
## 78 92 104 115 122
##
## lowest : 19 22 29 38 40, highest: 173 175 176 185 197
## --------------------------------------------------------------------------------
## TAS
## n missing distinct Info Mean Gmd .05 .10
## 12752 384871 152 1 118.6 23.38 87 93
## .25 .50 .75 .90 .95
## 105 118 131 145 156
##
## lowest : 37 48 50 51 52, highest: 221 222 225 229 238
## --------------------------------------------------------------------------------
## TAD
## n missing distinct Info Mean Gmd .05 .10
## 12755 384868 108 0.999 63.6 14.4 46 49
## .25 .50 .75 .90 .95
## 55 62 70 81 88
##
## lowest : 20 24 27 28 29, highest: 129 143 152 159 168
## --------------------------------------------------------------------------------
## SpO2
## n missing distinct Info Mean Gmd .05 .10
## 12755 384868 38 0.937 97.86 2.69 93 95
## .25 .50 .75 .90 .95
## 97 99 100 100 100
##
## lowest : 40 41 42 45 55, highest: 96 97 98 99 100
## --------------------------------------------------------------------------------
## FR
## n missing distinct Info Mean Gmd .05 .10
## 12762 384861 42 0.997 20.3 6.127 12 13
## .25 .50 .75 .90 .95
## 16 20 24 27 29
##
## lowest : 6 7 8 9 10, highest: 46 47 54 96 100
## --------------------------------------------------------------------------------
## GCSm
## n missing distinct Info Mean Gmd
## 5699 391924 6 0.852 2.863 2.242
##
## lowest : 1 2 3 4 5, highest: 2 3 4 5 6
##
## Value 1 2 3 4 5 6
## Frequency 2925 297 299 112 946 1120
## Proportion 0.513 0.052 0.052 0.020 0.166 0.197
## --------------------------------------------------------------------------------
## RASS
## n missing distinct Info Mean Gmd
## 9733 387890 8 0.94 -2.628 2.364
##
## lowest : -5 -4 -3 -2 -1, highest: -2 -1 0 1 2
##
## Value -5 -4 -3 -2 -1 0 1 2
## Frequency 3220 1216 862 601 1029 2651 106 48
## Proportion 0.331 0.125 0.089 0.062 0.106 0.272 0.011 0.005
## --------------------------------------------------------------------------------
## T Amb.
## n missing distinct Info Mean Gmd .05 .10
## 121595 276028 21 0.996 24.67 0.5823 24.0 24.1
## .25 .50 .75 .90 .95
## 24.2 24.6 25.0 25.4 25.6
##
## lowest : 24.0 24.1 24.2 24.3 24.4, highest: 25.6 25.7 25.8 25.9 26.0
## --------------------------------------------------------------------------------
## Luz Vis. 1
## n missing distinct Info Mean Gmd .05 .10
## 393750 3873 10198 1 33.08 45.71 0.05 0.18
## .25 .50 .75 .90 .95
## 1.37 10.75 44.40 100.00 148.00
##
## lowest : 0.00 0.01 0.02 0.03 0.04, highest: 296.00 297.00 298.00 299.00 300.00
## --------------------------------------------------------------------------------
## Luz Vis. 2
## n missing distinct Info Mean Gmd .05 .10
## 396024 1599 8882 1 9.387 13.4 0.01 0.04
## .25 .50 .75 .90 .95
## 0.34 2.67 11.85 27.43 41.40
##
## lowest : 0.00 0.01 0.02 0.03 0.04, highest: 288.00 289.00 292.00 298.00 300.00
## --------------------------------------------------------------------------------
## Luz IR 1
## n missing distinct Info Mean Gmd .05 .10
## 395792 1831 3053 0.986 3.779 6.347 0.00 0.00
## .25 .50 .75 .90 .95
## 0.01 0.35 1.81 11.54 21.00
##
## lowest : 0.00 0.01 0.02 0.03 0.04, highest: 287.00 292.00 297.00 299.00 300.00
## --------------------------------------------------------------------------------
## Pico Luz
## n missing distinct Info Mean Gmd .05 .10
## 393118 4505 16720 1 34.1 47 0.06 0.20
## .25 .50 .75 .90 .95
## 1.42 11.15 46.10 104.16 152.56
##
## lowest : 0.00 0.01 0.03 0.04 0.06, highest: 299.59 299.86 299.93 299.97 300.00
## --------------------------------------------------------------------------------
## Prom. dB
## n missing distinct Info Mean Gmd .05 .10
## 397623 0 569 1 62.8 9.384 48.4 49.7
## .25 .50 .75 .90 .95
## 56.7 64.0 66.4 70.0 80.0
##
## lowest : 25.0 30.6 31.1 31.6 33.9, highest: 106.3 106.5 107.4 107.9 110.0
## --------------------------------------------------------------------------------|- Análisis de Frecuencias:
También incluimos un análisis de frecuencias una vez las variables han sido procesadas y limpiadas. Esta clase de gráficos nos ayudan a hacernos una idea sobre la distribución de los valores obtenidos, si muestran un patrón de normalidad, si siguen existiendo datos extremos que modifican esta distribución…
plot_histogram(datos_logicos[ , 3:15],
ncol = 2,
ggtheme = theme_minimal())
Como norma general, vemos que el trabajo de pre-procesamiento de los datos ha dado su fruto y que se trata ahora de datos razonablemente limpios y lógicos, además de que todos parecen mostrar una normalidad clara. Sigue habiendo una serie de resultados que son dignos de mención y que deben ser tenidos en cuenta como, por ejemplo:
Fisiológicas:
FC: no presenta problemas, aunque hay registros de frecuencias cardíacas anormales como el mínimo de 19 o el máximo de 197. El análisis de frecuencias muestra un claro patrón de normalidad en torno a la media de FC = 91
TAS: Vemos una situación similar a FC. Hay una serie de valores extremos que no serían propios de un paciente estable y saludable, con mínimos de tensión de 37 y máximos de 238. No obstante, estos valores son extremos, la media de 118’6 es plausible, y el análisis de frecuencias indica un claro patrón de normalidad.
TAD: Situación también similar a TAS. Hay valores que mínimos de 20 y máximos de 168 qye son impropios de un paciente saludable. No obstante, son casos extremos, estando la media en TAD = 64 y viéndose en el análisis de frecuencias una distribución que se asemeja mucho a una normal.
SpO2: La saturación de Oxígeno no puede pasar del 100% por la naturaleza de la propia variable. Valores por debajo de 80% ya serían preocupantes, sin embargo, aparecen valores como 40% de saturación que no tienen sentido en este análisis. De nuevo, vemos que son poco frecuentes y que en su mayoría oscilan cerca de la media de 97’8%. El análisis de frecuencias muestra normalidad.
FR: Para la frecuencia respiratoria nos encontramos con una situación similar a FC, TAS, y TAD. Hay valores anormales como el mínimo de FR = 6 o el máximo de FR = 100. No obstante, estos son poco numerosos y en su mayoría (como se ve por los valores de los deciles) se encuentran al rededor de la media de FR = 20. El gráfico de análisis de frecuencias muestra también una distribución claramente normal y con poca incidencia de estos valores anormalmente bajos o altos.
GCSm y RASS: Estas variables son discretas y no se ven valores fuera de lo tipificado como normal.
Ambientales:
Tª Amb.: Esta variable ya venía avisando que no era completamente fiable. Vemos que los valores mínimos de 24ºC son aceptables, pero los máximos presentan temperaturas de 26ºC que no pueden ser realistas. La media, 25ºC, sí parece lógica, sin embargo, el gráfico de frecuencias nos indica que todavía existen valores extremos incluso después de la limpieza de datos. Esto, combinado con la acumulación de datos faltantes, vuelve a indicar la poca fiabilidad de esta variable predictora.
Luz Vis. 1,Luz Vis. 2, Luz IR 1, Pico Luz: Estas cuatro variables relacionadas con la luz presentan todas la misma casuística. Tienen valores mínimos esperables (cercanos a la oscuridad), sin embargo, los valores máximos son abiertamente imposibles incluso después del cribado y eliminación de datos atípicos que realizamos en secciones anteriores. El análisis de frecuencias nos indica, para las cuatro variables, que los valores más comunes entran dentro de lo esperable, sin embargo, valores extremos quizás resultado de una mala medición siguen ensuciando estas variables.
Prom. Db: Por último, el promedio de decibelios es una variable bastante estática y fiable que presenta valores mínimos dentro de lo razonable y valores extremos que, si bien no encajan en el entorno de una Unidad de Cuidados Intensivos, tampoco son descabellados ni muy frecuentes. La media está en 63 decibelios y las gráficas de análisis de frecuencias muestran una distribución que podría acercarse a la normalidad, aunque más bien apunta a una distribución con varias modas, probablemente resultado de diversos patrones de ruido: ruidos nocturnos, diurnos, y horarios de visitas.
|- Correlograma:
El objetivo final de este proyecto es encontrar relaciones de causalidad (si las hubiera) entre las variables ambientales registradas y el ‘outcome’, o sea, las variables fisiológicas que miden de forma indirecta la salud de los pacientes en UCI. Dada la gran cantidad de variables implicadas, lo apropiado sería reducir las posibles combinatorias mediante un análisis preliminar de correlaciones. Para ello, realizamos el siguiente correlograma:
my_colors <- brewer.pal(5, "PRGn")
my_colors <- colorRampPalette(my_colors)(100)
cor <- cor(na.omit(datos_logicos[ , 3:15]))
plotcorr(cor, col = my_colors[((cor + 1)/2) * 100],
lower.panel = "ellipse", upper.panel = "number", diag.panel = NULL,
mar = 0.1 + c(0,0,0,0))
En resumen, esta gráfica nos indica que:
Las variables ambientales de Luz (i.e., Luz1, Luz2, PicoLuz, luzIR) están altamente correlacionadas entre sí.
Las variables fisiológicas GCSm y RASS, que miden movilidad y nivel de consciencia del paciente, están correlacionadas positivamente con las variables cardíacas (i.e., FC, TAD, y TAS). Esto cae dentro de lo esperable dado que los pacientes sedados tendrán una actividad cardíaca más baja, y al revés. También están correlacionadas entre si.
También la Tensión Arteríal Sistólica y Diastólica están altamente correlacionadas entre si, como sería de esperar.
Por lo demás, la única variable prometedora que aparenta tener efecto sobre las variables fisiológicas medidas es el ruido (prom.dB) que sí que muestra una correlación positiva importante con FC, TAS, TAD, SpO2. También la Tª Ambiental, aunque es una variable poco fiable.
PARTE 3: FUNCIONALIZACIÓN DE LOS DATOS
Ahora pasamos a la fase de datos funcionales propiamente dicha. Las líneas de código siguientes son solo orientativas del procesado que requieren los datos y no requieren de mayor análisis por parte del lado Cordobés del proyecto. Es más bien para que tanto yo como Manuel tengamos de primera mano qué se está haciendo.
|- Objetos Dato-Funcional:
Crear los objetos funcionales per se requiere de un pre-procesado de los datos que se realiza a continuación. Comenzamos obteniendo las coordenadas exactas en las que se han registrado los datos fisiológicos.
#|- Coordenadas de filas non-NA ####
coord <- which(!is.na(datos_logicos$FC), arr.ind = TRUE) # donde se que hay datos
coord <- c(1, coord) # Para anhadir el punto inicial
head(coord, 20)## [1] 1 14 44 74 104 134 164 194 224 254 284 314 344 374 404 434 464 494 524
## [20] 554Una vez conocemos dichas coordendas, generamos un bucle recursivo en el que nos quedaremos con los datos fisiológicos registrados en ese punto, junto con los datos ambientales anteriores, junto con (además) los datos anteriores a los anteriores (\(tiempo-2\)), y también los valores promedio de las variables ambientales en ese tiempo. Quedarse con datos anteriores será estratégico para utilizarlos como datos explicativos de las mediciones fisiológicas del paciente, ya que parece lógico pensar que, por ejemplo, un factor importante para conocer el valor de Frecuencia Cardíaca en \(tiempo = 1\) será la propia frecuencia cardíaca anterior, esto es, la frecuencia cardíaca en \(tiempo=-1\).
Además, este bucle recursivo tiene incorporados una serie de estamentos \(IF-ELSE\) que comprobarán que:
- las curvas funcionales generadas tengan el tamaño adecuado
- los datos ambientales y fisiológicos de tiempos anteriores sigan pertenieciendo al mismo paciente.
Se muestra en el informe el proceso de generación de datos funcionales para la variable \(Temperatura\) ya que, para el resto de variables, el proceso es similar.
#|---- Temperatura: ####
df_temperatura <- as.data.frame(matrix(ncol = 72, nrow = 0))
# 72 porque son 10 variables fisiologicas de estudio (+ID+Hora+Label)
# + 30 datos de T-1 (+media2) + 30 datos de T-2 (+media2)
for (i in 3:length(coord)) {
aa <- t(as.matrix(datos_logicos$`T Amb.`[(coord[i - 1]):coord[i]]))
# Temperaturas tiempo-1
aa <- as.data.frame(t(aa[ , -1]))
bb <- datos_logicos[coord[i], c(1:9, 16)]
# Datos fisiologicos
cc <- mean(as.numeric(aa))
dd <- t(as.matrix(datos_logicos$`T Amb.`[(coord[i - 2]):coord[i - 1]]))
# Temperaturas tiempo-2
dd <- as.data.frame(t(dd[ , -1]))
ee <- mean(as.numeric(dd))
temp <- cbind(bb, aa, cc, dd, ee)
# Los empatamos
# Este if/else comprueba la longitud de los datos de ambos batches y que el
# anterior batch sea del mismo paciente. Si se cumplen, los añade al DF madre.
if (length(aa) != 30) {
next
} else if (length(dd) != 30) {
next
} else if (datos_logicos$ID[i] != datos_logicos$ID[i - 1]) {
next
} else {
df_temperatura <- rbind(df_temperatura, temp)
}
print(i)
}
colnames(df_temperatura) <- c("ID", "Fecha y hora", "FC", "TAS", "TAD", "SpO2",
"FR", "GCSm", "RASS", "Label", 1:30, "Promedio-1",
31:60, "Promedio-2")Además, las variables GCSm y RASS son registradas con menos frecuencia que el resto de variables fisiológicas, por lo que se va a asumir que el resultado se mantiene estático (igual al anterior) en los casos en que haya un NA.
# Hay muchos NA's intermedios en GCSm y RASS que se van a rellenar con el valor
# no-NA más cercano hacia arriba, o sea, el último registrado.
df_luz1$GCSm <- na.locf(df_luz1$GCSm, na.rm = F)
df_luz1$GCSm <- na.locf(df_luz1$GCSm, na.rm = F, fromLast = TRUE) # Y al revés
df_luz1$RASS <- na.locf(df_luz1$RASS, na.rm = F)
df_luz1$RASS <- na.locf(df_luz1$RASS, na.rm = F, fromLast = TRUE) # Y al revésEsto se realiza para todas las variables explicativas, aunque nos ahorramos poner aquí las líneas de código ya que son redundantes y, directamente, vamos a cargar los resultados finales.
save(list = c("df_luz1", "df_luz2", "df_luzIR", "df_PicoLuz", "df_promdB", "df_temperatura"),
file = "explicativas.RData")load("explicativas.RData")
load("datos_logicos.RData")Los objetos funcionales con _NA al final todavía conservan los NA y los valores que no son NA pero que se han marcado como outliers, mientras que los que no tienen _NA (e.g. ‘func_temp’) son solo con valores “consolidados”, o sea, que ni son outliers ni tienen NAs.
#|---- Temperatura
func_temp_NA = fdata(df_temperatura[ , 11:40],
names = list(main = "Temperatura con Outliers", xlab = "Tiempo",
ylab = "Cº"))
func_temp = fdata(df_temperatura[df_temperatura$Label == 0 & complete.cases(df_temperatura), 11:40],
names = list(main = "Temperatura sin Outliers", xlab = "Tiempo",
ylab = "Cº"))Ahora ya disponemos de nuestros datos ambientales en forma de dato funcional para proponer modelos de regresión y estimar la influencia de estas variables sobre las variables fisiológicas de interés.
|- Descriptivos Funcionales:
Los descriptivos funcionales podrían arrojar un poco más de luz sobre la naturaleza de los datos con los que trabajamos antes de empezar a proponer modelos.
#|---- Gráfico para sin Outliers: ####
index = 1:20; index2 = 500:520; index3 = 2000:2020
par(mfrow = c(3, 2))
plot.fdata(c(func_luz1[index], func_luz1[index2], func_luz1[index3]), col = "cyan")
plot.fdata(c(func_luz2[index], func_luz2[index2], func_luz2[index3]), col = "blue")
plot.fdata(c(func_luzIR[index], func_luzIR[index2], func_luzIR[index3]), col = "green")
plot.fdata(c(func_PicoLuz[index], func_PicoLuz[index2], func_PicoLuz[index3]), col = "darkgreen")
plot.fdata(c(func_promdB[index], func_promdB[index2], func_promdB[index3]), col = "red")
plot.fdata(c(func_temp[index], func_temp[index2], func_temp[index3]), col = "pink")
El análisis descriptivo funcional no arroja ninguna conclusión destacable que no pudieramos ya predecir con los análisis descriptivos anteriores o con los análisis de correlación. Se ve escasa variabilidad salvo unos pocos pacientes que sí parecen tener altibajos por alguna razón. También se puede ver como, al elegir el mismo índice para los tres primeros indicadores de luz, estos están muy correlacionados y presentan patrones similares.
La temperatura se ve que es sumamente estable, aunque también ha sufrido mucho como variable en cuestión de recorte de outliers y casos atípicos, por lo que volvemos a recordar que no se trata de una medida muy fiable. El único patrón mínimamente destacable es el de la variable de ruido, que parece tener presentar dos grupos distintos: uno de más ruido y otro de menos. Lo más probable es que esto tenga que ver con los horarios de visitas.
|- Pre-procesado ‘ldatos’:
ldatos <- ldata("df" = df_promdB[df_promdB$Label == 0 & complete.cases(df_promdB), ][-(3475:3476),],
temp = func_temp[1:3474], dB = func_promdB[1:3474],
luz1 = func_luz1[1:3474], luz2 = func_luz2[1:3474],
luzIR = func_luzIR[1:3474], pico = func_PicoLuz[1:3474])
# Por alguna razón hay una pequeña disparidad de longitud en dos objetos funcionales
# de modo que limito a 3474 curvas (se pierden 1 o 2 dependiendo del caso)
ldatos_NA <- ldata("df" = df_promdB,
temp = func_temp_NA, dB = func_promdB_NA, luz1 = func_luz1_NA,
luz2 = func_luz2_NA, luzIR = func_luzIR_NA, pico = func_PicoLuz_NA)
names(ldatos$df)[41] = "Promedio1"
names(ldatos$df)[72] = "Promedio2"
names(ldatos_NA$df)[41] = "Promedio1"
names(ldatos_NA$df)[72] = "Promedio2"|- Correlograma Funcional:
A <- cor(na.omit(datos_logicos[ , 3:15]))
names(ldatos_NA)## [1] "df" "temp" "dB" "luz1" "luz2" "luzIR" "pico"names(ldatos_NA$df)## [1] "ID" "Fecha y hora" "FC" "TAS" "TAD"
## [6] "SpO2" "FR" "GCSm" "RASS" "Label"
## [11] "1" "2" "3" "4" "5"
## [16] "6" "7" "8" "9" "10"
## [21] "11" "12" "13" "14" "15"
## [26] "16" "17" "18" "19" "20"
## [31] "21" "22" "23" "24" "25"
## [36] "26" "27" "28" "29" "30"
## [41] "Promedio1" "31" "32" "33" "34"
## [46] "35" "36" "37" "38" "39"
## [51] "40" "41" "42" "43" "44"
## [56] "45" "46" "47" "48" "49"
## [61] "50" "51" "52" "53" "54"
## [66] "55" "56" "57" "58" "59"
## [71] "60" "Promedio2"A1 <- ldatos_NA$temp
A2 <- ldatos_NA$dB
A3 <- ldatos_NA$luz1
A4 <- ldatos_NA$luz2
A5 <- ldatos_NA$luzIR
A6 <- ldatos_NA$pico
A7 <- as.matrix(ldatos_NA$df$FC)
A8 <- as.matrix(ldatos_NA$df$TAS)
A9 <- as.matrix(ldatos_NA$df$TAD)
A10 <- as.matrix(ldatos_NA$df$SpO2)
A11 <- as.matrix(ldatos_NA$df$FR)
A12 <- as.matrix(ldatos_NA$df$GCSm)
A13 <- as.matrix(ldatos_NA$df$RASS)
AA = A# FC:
A[1,2] <- dcor.xy(A7,A8)$estimate
A[1,3] <- dcor.xy(A7,A9)$estimate
A[1,4] <- dcor.xy(A7,A10)$estimate
A[1,5] <- dcor.xy(A7,A11)$estimate
A[1,6] <- dcor.xy(A7,A12)$estimate
A[1,7] <- dcor.xy(A7,A13)$estimate
A[1,8] <- dcor.xy(A7,A1)$estimate
A[1,9] <- dcor.xy(A7,A3)$estimate
A[1,10] <- dcor.xy(A7,A4)$estimate
A[1,11] <- dcor.xy(A7,A5)$estimate
A[1,12] <- dcor.xy(A7,A6)$estimate
A[1,13] <- dcor.xy(A7,A2)$estimate
# TAS:
A[2,1] <- dcor.xy(A8,A7)$estimate
A[2,3] <- dcor.xy(A8,A9)$estimate
A[2,4] <- dcor.xy(A8,A10)$estimate
A[2,5] <- dcor.xy(A8,A11)$estimate
A[2,6] <- dcor.xy(A8,A12)$estimate
A[2,7] <- dcor.xy(A8,A13)$estimate
A[2,8] <- dcor.xy(A8,A1)$estimate
A[2,9] <- dcor.xy(A8,A3)$estimate
A[2,10] <- dcor.xy(A8,A4)$estimate
A[2,11] <- dcor.xy(A8,A5)$estimate
A[2,12] <- dcor.xy(A8,A6)$estimate
A[2,13] <- dcor.xy(A8,A2)$estimate
# TAD:
A[3,1] <- dcor.xy(A9,A7)$estimate
A[3,2] <- dcor.xy(A9,A8)$estimate
A[3,4] <- dcor.xy(A9,A10)$estimate
A[3,5] <- dcor.xy(A9,A11)$estimate
A[3,6] <- dcor.xy(A9,A12)$estimate
A[3,7] <- dcor.xy(A9,A13)$estimate
A[3,8] <- dcor.xy(A9,A1)$estimate
A[3,9] <- dcor.xy(A9,A3)$estimate
A[3,10] <- dcor.xy(A9,A4)$estimate
A[3,11] <- dcor.xy(A9,A5)$estimate
A[3,12] <- dcor.xy(A9,A6)$estimate
A[3,13] <- dcor.xy(A9,A2)$estimate
# SpO2:
A[4,1] <- dcor.xy(A10,A7)$estimate
A[4,2] <- dcor.xy(A10,A8)$estimate
A[4,3] <- dcor.xy(A10,A9)$estimate
A[4,5] <- dcor.xy(A10,A11)$estimate
A[4,6] <- dcor.xy(A10,A12)$estimate
A[4,7] <- dcor.xy(A10,A13)$estimate
A[4,8] <- dcor.xy(A10,A1)$estimate
A[4,9] <- dcor.xy(A10,A3)$estimate
A[4,10] <- dcor.xy(A10,A4)$estimate
A[4,11] <- dcor.xy(A10,A5)$estimate
A[4,12] <- dcor.xy(A10,A6)$estimate
A[4,13] <- dcor.xy(A10,A2)$estimate
# FR:
A[5,1] <- dcor.xy(A11,A7)$estimate
A[5,2] <- dcor.xy(A11,A8)$estimate
A[5,3] <- dcor.xy(A11,A9)$estimate
A[5,4] <- dcor.xy(A11,A10)$estimate
A[5,6] <- dcor.xy(A11,A12)$estimate
A[5,7] <- dcor.xy(A11,A13)$estimate
A[5,8] <- dcor.xy(A11,A1)$estimate
A[5,9] <- dcor.xy(A11,A3)$estimate
A[5,10] <- dcor.xy(A11,A4)$estimate
A[5,11] <- dcor.xy(A11,A5)$estimate
A[5,12] <- dcor.xy(A11,A6)$estimate
A[5,13] <- dcor.xy(A11,A2)$estimate
# GCSm:
A[6,1] <- dcor.xy(A12,A7)$estimate
A[6,2] <- dcor.xy(A12,A8)$estimate
A[6,3] <- dcor.xy(A12,A9)$estimate
A[6,4] <- dcor.xy(A12,A10)$estimate
A[6,5] <- dcor.xy(A12,A12)$estimate
A[6,7] <- dcor.xy(A12,A13)$estimate
A[6,8] <- dcor.xy(A12,A1)$estimate
A[6,9] <- dcor.xy(A12,A3)$estimate
A[6,10] <- dcor.xy(A12,A4)$estimate
A[6,11] <- dcor.xy(A12,A5)$estimate
A[6,12] <- dcor.xy(A12,A6)$estimate
A[6,13] <- dcor.xy(A12,A2)$estimate
# RASS:
A[7,1] <- dcor.xy(A13,A7)$estimate
A[7,2] <- dcor.xy(A13,A8)$estimate
A[7,3] <- dcor.xy(A13,A9)$estimate
A[7,4] <- dcor.xy(A13,A10)$estimate
A[7,5] <- dcor.xy(A13,A11)$estimate
A[7,6] <- dcor.xy(A13,A12)$estimate
A[7,8] <- dcor.xy(A13,A1)$estimate
A[7,9] <- dcor.xy(A13,A3)$estimate
A[7,10] <- dcor.xy(A13,A4)$estimate
A[7,11] <- dcor.xy(A13,A5)$estimate
A[7,12] <- dcor.xy(A13,A6)$estimate
A[7,13] <- dcor.xy(A13,A2)$estimate
# TªAmb:
A[8,1] <- dcor.xy(A1,A7)$estimate
A[8,2] <- dcor.xy(A1,A8)$estimate
A[8,3] <- dcor.xy(A1,A9)$estimate
A[8,4] <- dcor.xy(A1,A10)$estimate
A[8,5] <- dcor.xy(A1,A11)$estimate
A[8,6] <- dcor.xy(A1,A12)$estimate
A[8,7] <- dcor.xy(A1,A13)$estimate
A[8,9] <- dcor.xy(A1,A3)$estimate
A[8,10] <- dcor.xy(A1,A4)$estimate
A[8,11] <- dcor.xy(A1,A5)$estimate
A[8,12] <- dcor.xy(A1,A6)$estimate
A[8,13] <- dcor.xy(A1,A2)$estimate
# Luz Vis. 1:
A[9,1] <- dcor.xy(A3,A7)$estimate
A[9,2] <- dcor.xy(A3,A8)$estimate
A[9,3] <- dcor.xy(A3,A9)$estimate
A[9,4] <- dcor.xy(A3,A10)$estimate
A[9,5] <- dcor.xy(A3,A11)$estimate
A[9,6] <- dcor.xy(A3,A12)$estimate
A[9,7] <- dcor.xy(A3,A13)$estimate
A[9,8] <- dcor.xy(A3,A1)$estimate
A[9,10] <- dcor.xy(A3,A4)$estimate
A[9,11] <- dcor.xy(A3,A5)$estimate
A[9,12] <- dcor.xy(A3,A6)$estimate
A[9,13] <- dcor.xy(A3,A2)$estimate
# Luz Vis. 2:
A[10,1] <- dcor.xy(A4,A7)$estimate
A[10,2] <- dcor.xy(A4,A8)$estimate
A[10,3] <- dcor.xy(A4,A9)$estimate
A[10,4] <- dcor.xy(A4,A10)$estimate
A[10,5] <- dcor.xy(A4,A11)$estimate
A[10,6] <- dcor.xy(A4,A12)$estimate
A[10,7] <- dcor.xy(A4,A13)$estimate
A[10,8] <- dcor.xy(A4,A1)$estimate
A[10,9] <- dcor.xy(A4,A3)$estimate
A[10,11] <- dcor.xy(A4,A5)$estimate
A[10,12] <- dcor.xy(A4,A6)$estimate
A[10,13] <- dcor.xy(A4,A2)$estimate
# Luz IR:
A[11,1] <- dcor.xy(A5,A7)$estimate
A[11,2] <- dcor.xy(A5,A8)$estimate
A[11,3] <- dcor.xy(A5,A9)$estimate
A[11,4] <- dcor.xy(A5,A10)$estimate
A[11,5] <- dcor.xy(A5,A11)$estimate
A[11,6] <- dcor.xy(A5,A12)$estimate
A[11,7] <- dcor.xy(A5,A13)$estimate
A[11,8] <- dcor.xy(A5,A1)$estimate
A[11,9] <- dcor.xy(A5,A3)$estimate
A[11,10] <- dcor.xy(A5,A4)$estimate
A[11,12] <- dcor.xy(A5,A6)$estimate
A[11,13] <- dcor.xy(A5,A2)$estimate
# Pico Luz:
A[12,1] <- dcor.xy(A6,A7)$estimate
A[12,2] <- dcor.xy(A6,A8)$estimate
A[12,3] <- dcor.xy(A6,A9)$estimate
A[12,4] <- dcor.xy(A6,A10)$estimate
A[12,5] <- dcor.xy(A6,A11)$estimate
A[12,6] <- dcor.xy(A6,A12)$estimate
A[12,7] <- dcor.xy(A6,A13)$estimate
A[12,8] <- dcor.xy(A6,A1)$estimate
A[12,9] <- dcor.xy(A6,A3)$estimate
A[12,10] <- dcor.xy(A6,A4)$estimate
A[12,11] <- dcor.xy(A6,A5)$estimate
A[12,13] <- dcor.xy(A6,A2)$estimate
# Prom. dB:
A[13,1] <- dcor.xy(A7,A7)$estimate
A[13,2] <- dcor.xy(A7,A8)$estimate
A[13,3] <- dcor.xy(A7,A9)$estimate
A[13,4] <- dcor.xy(A7,A10)$estimate
A[13,5] <- dcor.xy(A7,A11)$estimate
A[13,6] <- dcor.xy(A7,A12)$estimate
A[13,7] <- dcor.xy(A7,A13)$estimate
A[13,8] <- dcor.xy(A7,A1)$estimate
A[13,9] <- dcor.xy(A7,A3)$estimate
A[13,10] <- dcor.xy(A7,A4)$estimate
A[13,11] <- dcor.xy(A7,A5)$estimate
A[13,12] <- dcor.xy(A7,A6)$estimateround(A,2)
save(A, file = "corr_func.RData")Algo sigue pasando con el correlograma funcional que muchos valores resultan en NA’s. Falta hacerle una revisión aunque en realidad la idea es que sirva de apoyo al correlograma anteriormente mencionado. A falta de verlo entero, parece que apunta hacia las mismas conclusiones.
load("~/GitHub/FDA-CALUCI/corr_func.RData")
plotcorr(A, col = my_colors[((cor + 1)/2) * 100],
lower.panel = "ellipse", upper.panel = "number", diag.panel = NULL,
mar = 0.1 + c(0,0,0,0))
PARTE 4: PROPUESTA DE MODELOS
En este punto no se tiene claro realmente si utilizar modelos de regresión funcional va a ser útil y beneficioso en comparación con, por ejemplo, modelos GAM (Generalized Additive Models), o tratar los datos como series de tiempo… De modo que primero vamos a tantear un poco el terreno.
|- Toma de Contacto
||– Modelo Base Funcional
Este es un modelo de regresión funcional para dos variables que sabemos tienen una alta correlación para el estándar de este estudio. Utilizaremos su R-Sq y Varianza Explicada como referencias para ver si con otros métodos podemos igualar o superar estos resultados.
model.base <- fregre.gsam(FC ~ s(dB), data = ldatos)
summary(model.base)
# R-sq.(adj) = 0.212 Deviance explained = 21.9%||– Dependencia de T-1
Una idea que surgió es que a menudo el mejor predictor para la FC en un momento dado sería, probablemente, la FC anterior de ese paciente. Testeamos esta hipótesis y, efectivamente, si ajustamos un modelo de regresión funcional para FC en base a la FC anterior (FC0) y el promedio de decibelios, vemos que el R-sq es muy alto, pero que también lo es si eliminas los decibelios de la fórmula, lo cual indica que realmente el efecto total de una variable ambiental aquí es reducido.
Esto tiene mucho sentido ya que hablamos de una UCI en la que los pacientes tendrán un prognóstico mayormente influenciado por su propio estado de salud, no tanto por las variables ambientales. No obstante, lo que buscamos es precisamente los efectos de las variables ambientales en los pacientes ingresados, por lo que vamos a intentar girar hacia otros modelos que, aunque obtengan resultados más pobres, se centren en nuestras variables de interés. De todos modos, podría ser interesante comentar esto en un artículo ya que, además, se me ha comentado por varios investigadores en las conferencias de BIOSTATNET’23.
id = "A9"
ldatos0 = ldatos[ldatos$df$ID == id, row = T]
ldatos0 = ldatos0[2:132, row = T]
ldatos0$df$FC0 = ldatos[ldatos$df$ID == id, row = T]$df$FC[1:131]
modelT0.1 <- fregre.gsam(FC ~ s(FC0) + s(dB), data = ldatos0)
summary(modelT0.1)
# R-sq.(adj) = 0.868 Deviance explained = 89.8%
modelT0.2 <- fregre.gsam(FC ~ s(FC0), data = ldatos0)
summary(modelT0.2)
# R-sq.(adj) = 0.809 Deviance explained = 81.1%||– Base BSpline
Probamos creando unas bases de BSpline para ver si hay una mejoría. Las bases de B-spline se suelen usar para funciones no periódicas. Las funciones de base de B-spline son segmentos polinómicos unidos extremo a extremo en valores llamados nudos.
La idea era buena, pero los resultados indican que no existe ninguna mejoría y solamente incrementa el tiempo de computación.
lbsp <- list("dB" = create.bspline.basis(ldatos$dB$rangeval, 7))
model.spline <- fregre.gsam(FC ~ s(dB), data = ldatos, basis.b = lbsp)
summary(model.spline)
# R-sq.(adj) = 0.212 Deviance explained = 21.9%||– Base Componentes Principales
Con componentes principales (técnica de reducción de la dimensionalidad de un conjunto de datos con muchas variables en un conjunto de variables linealmente correlacionadas conocidas como componentes principales) los resultados son buenos pero no tan buenos como con una regresión funcional como las que hemos estado viendo hasta ahora. No parece que vaya a mejorar el análisis.
lpc <- list("dB" = create.pc.basis(ldatos$dB, 1:7))
model.PC <- fregre.gsam(FC ~ s(dB), data = ldatos, basis.x = lpc)
summary(model.PC)
# R-sq.(adj) = 0.199 Deviance explained = 20.5%||– Series de Tiempo
Probamos una aproximación con series de tiempo testeando para varios pacientes. En general vemos que no tiene una mejor utilidad y que lo único que hace es reforzar nuestra observación de que en gran medida las variables fisiológicas vienen muy determinadas por la propia variable en tiempos anteriores.
ID = "B5"
k = df_luz1$FC[df_luz1$Label == 0 & df_luz1$ID == ID]
ts.model <- gam(k[-1] ~ s(k[-length(k)]))
summary(ts.model)plot.ts(k, main = paste0("Frecuencia Cardiaca paciente ", as.character(ID)))
ar(k, aic = TRUE, order.max = NULL) ##
## Call:
## ar(x = k, aic = TRUE, order.max = NULL)
##
## Coefficients:
## 1
## 0.863
##
## Order selected 1 sigma^2 estimated as 5.04# A1: R-sq.(adj) = 0.57 / Deviance explained = 58.1% / Coefficients = 0.745
# B5: R-sq.(adj) = 0.79 / Deviance explained = 79.3% / Coefficients = 0.863
# C5: R-sq.(adj) = 0.589 / Deviance explained = 59.7% / Coefficients = 0.737
# D10: R-sq.(adj)= 0.888 / Deviance explained = 88.9% / Coefficients = 0.942|- Modelos de Regresión Funcional
||– FC ~
Para la frecuencia cardíaca vemos que el único modelo que parece tener interés es el de FC~dB. El resto tienen un R-Sq muy bajo. Añadiremos este modelo de interés a la fase final de comparación entre metodologías.
model01.1 <- fregre.gsam(FC ~ s(temp), data = ldatos)
summary(model01.1)$r.sq; AIC(model01.1)
# R-sq.(adj) = 0.01421 AIC = 30895
model01.2 <- fregre.gsam(FC ~ s(dB), data = ldatos)
summary(model01.2)$r.sq; AIC(model01.2)
# R-sq.(adj) = 0.2124 AIC = 30129 (1º)
model01.3 <- fregre.gsam(FC ~ s(luz1), data = ldatos)
summary(model01.3)$r.sq; AIC(model01.3)
# R-sq.(adj) = 0.04429 AIC = 30793 (8º)
model01.4 <- fregre.gsam(FC ~ s(luz2), data = ldatos)
summary(model01.4)$r.sq; AIC(model01.4)
# R-sq.(adj) = 0.04355 AIC = 30794
model01.5 <- fregre.gsam(FC ~ s(luzIR), data = ldatos)
summary(model01.5)$r.sq; AIC(model01.5)
# R-sq.(adj) = 0.02137 AIC = 30866
model01.6 <- fregre.gsam(FC ~ s(pico), data = ldatos)
summary(model01.6)$r.sq; AIC(model01.6)
# R-sq.(adj) = 0.04262 AIC = 30799||– TAS ~
Para la Tensión Arterial Sistólica de nuevo la variable ambiental de ruido es la que más influencia tiene, aunque esta vez menor que para la FC (TAS~dB). El resto de variables no llegan ni de lejos a un 0’1 de R-Sq por lo que no serán valoradas para la sección final de este informe.
model02.1 <- fregre.gsam(TAS ~ s(temp), data = ldatos)
summary(model02.1)$r.sq; AIC(model02.1)
# R-sq.(adj) = 0.0204 AIC = 31207
model02.2 <- fregre.gsam(TAS ~ s(dB), data = ldatos)
summary(model02.2)$r.sq; AIC(model02.2)
# R-sq.(adj) = 0.1636 AIC = 30664 (3º)
model02.3 <- fregre.gsam(TAS ~ s(luz1), data = ldatos)
summary(model02.3)$r.sq; AIC(model02.3)
# R-sq.(adj) = 0.02045 AIC = 31205
model02.4 <- fregre.gsam(TAS ~ s(luz2), data = ldatos)
summary(model02.4)$r.sq; AIC(model02.4)
# R-sq.(adj) = 0.02153 AIC = 31201 (12º)
model02.5 <- fregre.gsam(TAS ~ s(luzIR), data = ldatos)
summary(model02.5)$r.sq; AIC(model02.5)
# R-sq.(adj) = 0.007196 AIC = 31250
model02.6 <- fregre.gsam(TAS ~ s(pico), data = ldatos)
summary(model02.6)$r.sq; AIC(model02.6)
# R-sq.(adj) = 0.01854 AIC = 31215||– TAD ~
De nuevo es el ruido lo que más efecto tiene, aunque menos todavía que con TAS. El modelo que valoraremos - si acaso - para la fase final del estudio será TAD~dB. Es curioso y quizás requiera un poco de meditación el ver que el ruido tiene un gran efecto sobre la FC, menos sobre la TAS, y menos todavía sobre la TAD. Por qué podría ser esto?
model03.1 <- fregre.gsam(TAD ~ s(temp), data = ldatos)
summary(model03.1)$r.sq; AIC(model03.1)
# R-sq.(adj) = 0.02274 AIC = 26845 (11º)
model03.2 <- fregre.gsam(TAD ~ s(dB), data = ldatos)
summary(model03.2)$r.sq; AIC(model03.2)
# R-sq.(adj) = 0.08999 AIC = 26593 (5º)
model03.3 <- fregre.gsam(TAD ~ s(luz1), data = ldatos)
summary(model03.3)$r.sq; AIC(model03.3)
# R-sq.(adj) = 0.01926 AIC = 26838
model03.4 <- fregre.gsam(TAD ~ s(luz2), data = ldatos)
summary(model03.4)$r.sq; AIC(model03.4)
# R-sq.(adj) = 0.01437 AIC = 26857
model03.5 <- fregre.gsam(TAD ~ s(luzIR), data = ldatos)
summary(model03.5)$r.sq; AIC(model03.5)
# R-sq.(adj) = 0.01553 AIC = 26850
model03.6 <- fregre.gsam(TAD ~ s(pico), data = ldatos)
summary(model03.6)$r.sq; AIC(model03.6)
# R-sq.(adj) = 0.01757 AIC = 26841||– SpO2 ~
De nuevo vemos que la única variable que podría considerarse de interés es la implicada en el modelo SpO2~dB. Tampoco tiene un R-Sq destacable pero es la única que podría entrar en la fase final del informe.
model04.1 <- fregre.gsam(SpO2 ~ s(temp), data = ldatos)
summary(model04.1)$r.sq; AIC(model04.1)
# R-sq.(adj) = 0.003161 AIC = 19385
model04.2 <- fregre.gsam(SpO2 ~ s(dB), data = ldatos)
summary(model04.2)$r.sq; AIC(model04.2)
# R-sq.(adj) = 0.0495 AIC = 19236 (7º)
model04.3 <- fregre.gsam(SpO2 ~ s(luz1), data = ldatos)
summary(model04.3)$r.sq; AIC(model04.3)
# R-sq.(adj) = 0.02735 AIC = 19304
model04.4 <- fregre.gsam(SpO2 ~ s(luz2), data = ldatos)
summary(model04.4)$r.sq; AIC(model04.4)
# R-sq.(adj) = 0.03055 AIC = 19291 (10º)
model04.5 <- fregre.gsam(SpO2 ~ s(luzIR), data = ldatos)
summary(model04.5)$r.sq; AIC(model04.5)
# R-sq.(adj) = 0.01831 AIC = 19331
model04.6 <- fregre.gsam(SpO2 ~ s(pico), data = ldatos)
summary(model04.6)$r.sq; AIC(model04.6)
# R-sq.(adj) = 0.02869 AIC = 19300||– FR ~
De nuevo, es el ruido la variable que destaca. En el modelo FR~dB se alcanza un R-Sq de 0’073 y será el modelo que utilizaremos en la última sección del informe.
model05.1 <- fregre.gsam(FR ~ s(temp), data = ldatos)
summary(model05.1)$r.sq; AIC(model05.1)
# R-sq.(adj) = 0.005235 AIC = 21301
model05.2 <- fregre.gsam(FR ~ s(dB), data = ldatos)
summary(model05.2)$r.sq; AIC(model05.2)
# R-sq.(adj) = 0.07371 AIC = 21057 (6º)
model05.3 <- fregre.gsam(FR ~ s(luz1), data = ldatos)
summary(model05.3)$r.sq; AIC(model05.3)
# R-sq.(adj) = 0.03965 AIC = 21196 (9º)
model05.4 <- fregre.gsam(FR ~ s(luz2), data = ldatos)
summary(model05.4)$r.sq; AIC(model05.4)
# R-sq.(adj) = 0.03659 AIC = 21205
model05.5 <- fregre.gsam(FR ~ s(luzIR), data = ldatos)
summary(model05.5)$r.sq; AIC(model05.5)
# R-sq.(adj) = 0.002133 AIC = 21308
model05.6 <- fregre.gsam(FR ~ s(pico), data = ldatos)
summary(model05.6)$r.sq; AIC(model05.6)
# R-sq.(adj) = 0.03081 AIC = 21223||– GCSm ~
Para GCSm volvemos a encontrar que es GCSm~dB el modelo con los mejores resultados y que pasará a la fase de predicción. Empiezo a pensar que este patrón va a ser fuerte y que deberíamos quizás incluir también la 2ª variable con más impacto en los análisis o nos terminaremos quedando solamente con el ruido.
model06.1 <- fregre.gsam(GCSm ~ s(temp), data = ldatos)
summary(model06.1)$r.sq; AIC(model06.1)
# R-sq.(adj) = 0.0145 AIC = 14235 (14º)
model06.2 <- fregre.gsam(GCSm ~ s(dB), data = ldatos)
summary(model06.2)$r.sq; AIC(model06.2)
# R-sq.(adj) = 0.128 AIC = 13818 (3º)
model06.3 <- fregre.gsam(GCSm ~ s(luz1), data = ldatos)
summary(model06.3)$r.sq; AIC(model06.3)
# R-sq.(adj) = 0.0107 AIC = 14253
model06.4 <- fregre.gsam(GCSm ~ s(luz2), data = ldatos)
summary(model06.4)$r.sq; AIC(model06.4)
# R-sq.(adj) = 0.005607 AIC = 14267
model06.5 <- fregre.gsam(GCSm ~ s(luzIR), data = ldatos)
summary(model06.5)$r.sq; AIC(model06.5)
# R-sq.(adj) = 0.0008525 AIC = 14282
model06.6 <- fregre.gsam(GCSm ~ s(pico), data = ldatos)
summary(model06.6)$r.sq; AIC(model06.6)
# R-sq.(adj) = 0.01344 AIC = 14251||– RASS ~
De nuevo es el ruido. RASS~dB es el modelo con un R-Sq superior.
model07.1 <- fregre.gsam(RASS ~ s(temp), data = ldatos)
summary(model07.1)$r.sq; AIC(model07.1)
# R-sq.(adj) = 0.006167 AIC = 14174
model07.2 <- fregre.gsam(RASS ~ s(dB), data = ldatos)
summary(model07.2)$r.sq; AIC(model07.2)
# R-sq.(adj) = 0.1213 AIC = 13757 (4º)
model07.3 <- fregre.gsam(RASS ~ s(luz1), data = ldatos)
summary(model07.3)$r.sq; AIC(model07.3)
# R-sq.(adj) = 0.01561 AIC = 14145
model07.4 <- fregre.gsam(RASS ~ s(luz2), data = ldatos)
summary(model07.4)$r.sq; AIC(model07.4)
# R-sq.(adj) = 0.01944 AIC = 14128 (13º)
model07.5 <- fregre.gsam(RASS ~ s(luzIR), data = ldatos)
summary(model07.5)$r.sq; AIC(model07.5)
# R-sq.(adj) = 0.0194 AIC = 14130
model07.6 <- fregre.gsam(RASS ~ s(pico), data = ldatos)
summary(model07.6)$r.sq; AIC(model07.6)
# R-sq.(adj) = 0.01558 AIC = 14141Ahora hemos obtenido una combinatoria de modelos y sabemos de forma orientativa el nivel de ajuste y adecuación que estos tienen para nuestros datos. Para cada variable fisiológica hemos elegido los dos mejores modelos y los hemos puesto en una tabla para comparar si efectivamente usar regresión funcional es el mejor método o si habríamos obtenido resultados similares o mejores con modelos más simples.
|- Comparativa de Modelos Seleccionados:
La tabla que se adjunta actúa de resumen de toda esta sección. Vemos que efectivamente los modelos de regresión funcional GSAM son (salvo para un modelo) los que obtienen el mejor R-sq y en todos salvo en tres casos el mejor AIC. No obstante, vemos que son mejorías de modelización marginales aunque constantes. La mejoría que implica utilizar datos funcionales no es muy grande, sin embargo parece evidente que sí que la hay y no existen razones para rechazar este pequeño avance en la capacidad de modelización.
||– Tabla Resumen:
||– Lineal Models:
lm.01 <- lm(ldatos$df$FC ~ ldatos$dB$data)
summary(lm.01)$r.sq ; AIC(lm.01)
lm.02 <- lm(ldatos$df$TAS ~ ldatos$dB$data)
summary(lm.02)$r.sq ; AIC(lm.02)
lm.03 <- lm(ldatos$df$GCSm ~ ldatos$dB$data)
summary(lm.03)$r.sq ; AIC(lm.03)
lm.04 <- lm(ldatos$df$RASS ~ ldatos$dB$data)
summary(lm.04)$r.sq ; AIC(lm.04)
lm.05 <- lm(ldatos$df$TAD ~ ldatos$dB$data)
summary(lm.05)$r.sq ; AIC(lm.05)
lm.06 <- lm(ldatos$df$FR ~ ldatos$dB$data)
summary(lm.06)$r.sq ; AIC(lm.06)
lm.07 <- lm(ldatos$df$SpO2 ~ ldatos$dB$data)
summary(lm.07)$r.sq ; AIC(lm.07)
lm.08 <- lm(ldatos$df$FC ~ ldatos$luz1$data)
summary(lm.08)$r.sq ; AIC(lm.08)
lm.09 <- lm(ldatos$df$FR ~ ldatos$luz1$data)
summary(lm.09)$r.sq ; AIC(lm.09)
lm.10 <- lm(ldatos$df$SpO2 ~ ldatos$luz2$data)
summary(lm.10)$r.sq ; AIC(lm.10)
lm.11 <- lm(ldatos$df$TAD ~ ldatos$temp$data)
summary(lm.11)$r.sq ; AIC(lm.11)
lm.12 <- lm(ldatos$df$TAS ~ ldatos$luz2$data)
summary(lm.12)$r.sq ; AIC(lm.12)
lm.13 <- lm(ldatos$df$RASS ~ ldatos$luz2$data)
summary(lm.13)$r.sq ; AIC(lm.13)
lm.14 <- lm(ldatos$df$GCSm ~ ldatos$temp$data)
summary(lm.14)$r.sq ; AIC(lm.14)||– Functional Linear Models:
func.lin.01 <- fregre.lm(FC ~ dB, data = ldatos)
summary(func.lin.01)$r.sq ; AIC(func.lin.01)
func.lin.02 <- fregre.lm(TAS ~ dB, data = ldatos)
summary(func.lin.02)$r.sq ; AIC(func.lin.02)
func.lin.03 <- fregre.lm(GCSm ~ dB, data = ldatos)
summary(func.lin.03)$r.sq ; AIC(func.lin.03)
func.lin.04 <- fregre.lm(RASS ~ dB, data = ldatos)
summary(func.lin.04)$r.sq ; AIC(func.lin.04)
func.lin.05 <- fregre.lm(TAD ~ dB, data = ldatos)
summary(func.lin.05)$r.sq ; AIC(func.lin.05)
func.lin.06 <- fregre.lm(FR ~ dB, data = ldatos)
summary(func.lin.06)$r.sq ; AIC(func.lin.06)
func.lin.07 <- fregre.lm(SpO2 ~ dB, data = ldatos)
summary(func.lin.07)$r.sq ; AIC(func.lin.07)
func.lin.08 <- fregre.lm(FC ~ luz1, data = ldatos)
summary(func.lin.08)$r.sq ; AIC(func.lin.08)
func.lin.09 <- fregre.lm(FR ~ luz1, data = ldatos)
summary(func.lin.09)$r.sq ; AIC(func.lin.09)
func.lin.10 <- fregre.lm(SpO2 ~ luz2, data = ldatos)
summary(func.lin.10)$r.sq ; AIC(func.lin.10)
func.lin.11 <- fregre.lm(TAD ~ temp, data = ldatos)
summary(func.lin.11)$r.sq ; AIC(func.lin.11)
func.lin.12 <- fregre.lm(TAS ~ luz2, data = ldatos)
summary(func.lin.12)$r.sq ; AIC(func.lin.12)
func.lin.13 <- fregre.lm(RASS ~ luz2, data = ldatos)
summary(func.lin.13)$r.sq ; AIC(func.lin.13)
func.lin.14 <- fregre.lm(GCSm ~ temp, data = ldatos)
summary(func.lin.14)$r.sq ; AIC(func.lin.14)||– Generalized Additive Models:
gen.add.01 <- gam(ldatos$df$FC ~ s(ldatos$dB$data))
summary(gen.add.01)$r.sq ; AIC(gen.add.01)
gen.add.02 <- gam(ldatos$df$TAS ~ s(ldatos$dB$data))
summary(gen.add.02)$r.sq ; AIC(gen.add.02)
gen.add.03 <- gam(ldatos$df$GCSm ~ s(ldatos$dB$data))
summary(gen.add.03)$r.sq ; AIC(gen.add.03)
gen.add.04 <- gam(ldatos$df$RASS ~ s(ldatos$dB$data))
summary(gen.add.04)$r.sq ; AIC(gen.add.04)
gen.add.05 <- gam(ldatos$df$TAD ~ s(ldatos$dB$data))
summary(gen.add.05)$r.sq ; AIC(gen.add.05)
gen.add.06 <- gam(ldatos$df$FR ~ s(ldatos$dB$data))
summary(gen.add.06)$r.sq ; AIC(gen.add.06)
gen.add.07 <- gam(ldatos$df$SpO2 ~ s(ldatos$dB$data))
summary(gen.add.07)$r.sq ; AIC(gen.add.07)
gen.add.08 <- gam(ldatos$df$FC ~ s(ldatos$luz1$data))
summary(gen.add.08)$r.sq ; AIC(gen.add.08)
gen.add.09 <- gam(ldatos$df$FR ~ s(ldatos$luz1$data))
summary(gen.add.09)$r.sq ; AIC(gen.add.09)
gen.add.10 <- gam(ldatos$df$SpO2 ~ s(ldatos$luz2$data))
summary(gen.add.10)$r.sq ; AIC(gen.add.10)
gen.add.11 <- gam(ldatos$df$TAD ~ s(ldatos$temp$data))
summary(gen.add.11)$r.sq ; AIC(gen.add.11)
gen.add.12 <- gam(ldatos$df$TAS ~ s(ldatos$luz2$data))
summary(gen.add.12)$r.sq ; AIC(gen.add.12)
gen.add.13 <- gam(ldatos$df$RASS ~ s(ldatos$luz2$data))
summary(gen.add.13)$r.sq ; AIC(gen.add.13)
gen.add.14 <- gam(ldatos$df$GCSm ~ s(ldatos$temp$data))
summary(gen.add.14)$r.sq ; AIC(gen.add.14)||– Functional Generalized Spectral Additive Models:
func.gen.add.01 <- fregre.gsam(FC ~ s(dB), data = ldatos)
summary(func.gen.add.01)$r.sq ; AIC(func.gen.add.01)
func.gen.add.02 <- fregre.gsam(TAS ~ s(dB), data = ldatos)
summary(func.gen.add.02)$r.sq ; AIC(func.gen.add.02)
func.gen.add.03 <- fregre.gsam(GCSm ~ s(dB), data = ldatos)
summary(func.gen.add.03)$r.sq ; AIC(func.gen.add.03)
func.gen.add.04 <- fregre.gsam(RASS ~ s(dB), data = ldatos)
summary(func.gen.add.04)$r.sq ; AIC(func.gen.add.04)
func.gen.add.05 <- fregre.gsam(TAD ~ s(dB), data = ldatos)
summary(func.gen.add.05)$r.sq ; AIC(func.gen.add.05)
func.gen.add.06 <- fregre.gsam(FR ~ s(dB), data = ldatos)
summary(func.gen.add.06)$r.sq ; AIC(func.gen.add.06)
func.gen.add.07 <- fregre.gsam(SpO2 ~ s(dB), data = ldatos)
summary(func.gen.add.07)$r.sq ; AIC(func.gen.add.07)
func.gen.add.08 <- fregre.gsam(FC ~ s(luz1), data = ldatos)
summary(func.gen.add.08)$r.sq ; AIC(func.gen.add.08)
func.gen.add.09 <- fregre.gsam(FR ~ s(luz1), data = ldatos)
summary(func.gen.add.09)$r.sq ; AIC(func.gen.add.09)
func.gen.add.10 <- fregre.gsam(SpO2 ~ s(luz2), data = ldatos)
summary(func.gen.add.10)$r.sq ; AIC(func.gen.add.10)
func.gen.add.11 <- fregre.gsam(TAD ~ s(temp), data = ldatos)
summary(func.gen.add.11)$r.sq ; AIC(func.gen.add.11)
func.gen.add.12 <- fregre.gsam(TAS ~ s(luz2), data = ldatos)
summary(func.gen.add.12)$r.sq ; AIC(func.gen.add.12)
func.gen.add.13 <- fregre.gsam(RASS ~ s(luz2), data = ldatos)
summary(func.gen.add.13)$r.sq ; AIC(func.gen.add.13)
func.gen.add.14 <- fregre.gsam(GCSm ~ s(temp), data = ldatos)
summary(func.gen.add.14)$r.sq ; AIC(func.gen.add.14)PARTE 5: PREDICCIÓN EN BASE A MODELOS
Dado que el objetivo de este estudio se centraba en el uso de datos funcionales y que en 13 de los 14 modelos los datos funcionales superan a sus modelos competidores, continuaré a la fase de predicción en base a modelos únicamente utilizando estos modelos de regresión funcional. Existe un único caso en que marginalmente los GAM tuvieron un mejor comportamiento, pero esto complicaría la comprensión de los resultados generales.
# Filtrar ldatos0 por tiempos:
datos_pred <- df_promdB[year(df_promdB$`Fecha y hora`) == 2021 & df_promdB$Label == 0,
1:9]
datos_pred$predict <- FALSE # Variable booleana para marcar qué se va a predecir
ID <- as.vector(t(unique(datos_pred[, 1])))
coord = 0 # marcapasos
for (i in 1:length(ID)) {
n <- nrow(datos_pred[datos_pred$ID == ID[i],]) # coordenadas de cada ID
datos_pred[coord + n , 10] = TRUE # el último y el
datos_pred[coord + n - 1, 10] = TRUE # penúltimo valor pasan a ser para predecir
coord <- coord + n # y se actualiza el marcapasos
}
# Crear nuevos objetos funcionales:
func_temp_pred = fdata(df_temperatura[year(df_temperatura$`Fecha y hora`) == 2021 & df_temperatura$Label == 0, 11:40],
names = list(main = "Temperatura con Outliers", xlab = "Tiempo",
ylab = "Cº"))
func_luz1_pred = fdata(df_luz1[year(df_luz1$`Fecha y hora`) == 2021 & df_luz1$Label == 0, 11:40],
names = list(main = "Luz 1 con Outliers", xlab = "Tiempo",
ylab = "Luz 1 (lumens)"))
func_luz2_pred = fdata(df_luz2[year(df_luz2$`Fecha y hora`) == 2021 & df_luz2$Label == 0, 11:40],
names = list(main = "Luz 2 con Outliers", xlab = "Tiempo",
ylab = "Luz 2 (lumens)"))
func_luzIR_pred = fdata(df_luzIR[year(df_luzIR$`Fecha y hora`) == 2021 & df_luzIR$Label == 0, 11:40],
names = list(main = "Luz IR con Outliers", xlab = "Tiempo",
ylab = "Luz IR (lumens)"))
func_PicoLuz_pred = fdata(df_PicoLuz[year(df_PicoLuz$`Fecha y hora`) == 2021 & df_PicoLuz$Label == 0, 11:40],
names = list(main = "Pico Luz con Outliers", xlab = "Tiempo",
ylab = "Pico Luz (lumens)"))
func_promdB_pred = fdata(df_promdB[year(df_promdB$`Fecha y hora`) == 2021 & df_promdB$Label == 0, 11:40],
names = list(main = "Ruido con Outliers", xlab = "Tiempo",
ylab = "Decibelios (dB)"))
# Crear ldatos_pred:
ldatos_pred <- ldata("df" = datos_pred,
temp = func_temp_pred, dB = func_promdB_pred,
luz1 = func_luz1_pred, luz2 = func_luz2_pred,
luzIR = func_luzIR_pred, pico = func_PicoLuz_pred)# Ajustar modelos:
new.coord <- which(datos_pred$predict == TRUE, arr.ind = TRUE)
ldatos_pred_FALSE <- ldata("df" = datos_pred[datos_pred$predict == FALSE ,],
temp = func_temp_pred[-new.coord],
dB = func_promdB_pred[-new.coord],
luz1 = func_luz1_pred[-new.coord],
luz2 = func_luz2_pred[-new.coord],
luzIR = func_luzIR_pred[-new.coord],
pico = func_PicoLuz_pred[-new.coord])
func.gen.add.01 <- fregre.gsam(FC ~ s(dB), data = ldatos_pred_FALSE)
func.gen.add.02 <- fregre.gsam(TAS ~ s(dB), data = ldatos_pred_FALSE)
func.gen.add.03 <- fregre.gsam(GCSm ~ s(dB), data = ldatos_pred_FALSE)
func.gen.add.04 <- fregre.gsam(RASS ~ s(dB), data = ldatos_pred_FALSE)
func.gen.add.05 <- fregre.gsam(TAD ~ s(dB), data = ldatos_pred_FALSE)
func.gen.add.06 <- fregre.gsam(FR ~ s(dB), data = ldatos_pred_FALSE)
func.gen.add.07 <- fregre.gsam(SpO2 ~ s(dB), data = ldatos_pred_FALSE)
func.gen.add.08 <- fregre.gsam(FC ~ s(luz1), data = ldatos_pred_FALSE)
func.gen.add.09 <- fregre.gsam(FR ~ s(luz1), data = ldatos_pred_FALSE)
func.gen.add.10 <- fregre.gsam(SpO2 ~ s(luz2), data = ldatos_pred_FALSE)
func.gen.add.11 <- fregre.gsam(TAD ~ s(temp), data = ldatos_pred_FALSE)
func.gen.add.12 <- fregre.gsam(TAS ~ s(luz2), data = ldatos_pred_FALSE)
func.gen.add.13 <- fregre.gsam(RASS ~ s(luz2), data = ldatos_pred_FALSE)
func.gen.add.14 <- fregre.gsam(GCSm ~ s(temp), data = ldatos_pred_FALSE)ID <- as.vector(t(unique(datos_pred$ID)))
# Predicción por pacientes + RMSE:
ldatos_pred_TRUE <- ldata("df" = datos_pred[datos_pred$predict == TRUE ,],
temp = func_temp_pred, dB = func_promdB_pred,
luz1 = func_luz1_pred, luz2 = func_luz2_pred,
luzIR = func_luzIR_pred, pico = func_PicoLuz_pred)
# Creamos tabla para anidar resultados:
RMSE <- data.frame(model = integer(), RMSE = numeric(), NRMSE = numeric())
# length(ID) = 46, por lo que al final RMSE tiene que tener 46 valores
# multiplicado por 14 modelos = 644 valores.
model = "FC~Prom.dB"
for (i in 1:length(ID)) {
newdata <- ldatos_pred_TRUE[(i*2 - 1):(i*2), row = T] # Dos datos por paciente
pred <- predict(func.gen.add.01, newdata)
rmse <- pred2meas(ldatos_pred_TRUE$df$FC[(i*2 - 1):(i*2)], as.numeric(pred), measure = "RMSE")
nrmse = sqrt(mean(((as.numeric(pred)) - (ldatos_pred_TRUE$df$FC[(i*2 - 1):(i*2)]))^{2})/var(ldatos_pred$df$FC))
xx <- cbind(model, rmse, nrmse)
RMSE <- rbind(RMSE, xx)
}
model = "TAS~Prom.dB"
for (i in 1:length(ID)) {
newdata <- ldatos_pred_TRUE[(i*2 - 1):(i*2), row = T]
pred <- predict(func.gen.add.02, newdata)
rmse <- pred2meas(ldatos_pred_TRUE$df$TAS[(i*2 - 1):(i*2)], pred, measure = "RMSE")
nrmse = sqrt(mean(((as.numeric(pred)) - (ldatos_pred_TRUE$df$TAS[(i*2 - 1):(i*2)]))^{2})/var(ldatos_pred$df$TAS))
xx <- cbind(model, rmse, nrmse)
RMSE <- rbind(RMSE, xx)
}
model = "GCSm~Prom.dB"
for (i in 1:length(ID)) {
newdata <- ldatos_pred_TRUE[(i*2 - 1):(i*2), row = T]
pred <- predict(func.gen.add.03, newdata)
rmse <- pred2meas(ldatos_pred_TRUE$df$GCSm[(i*2 - 1):(i*2)], pred, measure = "RMSE")
nrmse = sqrt(mean(((as.numeric(pred)) - (ldatos_pred_TRUE$df$GCSm[(i*2 - 1):(i*2)]))^{2})/var(ldatos_pred$df$GCSm))
xx <- cbind(model, rmse, nrmse)
RMSE <- rbind(RMSE, xx)
}
model = "RASS~Prom.dB"
for (i in 1:length(ID)) {
newdata <- ldatos_pred_TRUE[(i*2 - 1):(i*2), row = T]
pred <- predict(func.gen.add.04, newdata)
rmse <- pred2meas(ldatos_pred_TRUE$df$RASS[(i*2 - 1):(i*2)], pred, measure = "RMSE")
nrmse = sqrt(mean(((as.numeric(pred)) - (ldatos_pred_TRUE$df$RASS[(i*2 - 1):(i*2)]))^{2})/var(ldatos_pred$df$RASS))
xx <- cbind(model, rmse, nrmse)
RMSE <- rbind(RMSE, xx)
}
model = "TAD~Prom.dB"
for (i in 1:length(ID)) {
newdata <- ldatos_pred_TRUE[(i*2 - 1):(i*2), row = T]
pred <- predict(func.gen.add.05, newdata)
rmse <- pred2meas(ldatos_pred_TRUE$df$TAD[(i*2 - 1):(i*2)], pred, measure = "RMSE")
nrmse = sqrt(mean(((as.numeric(pred)) - (ldatos_pred_TRUE$df$TAD[(i*2 - 1):(i*2)]))^{2})/var(ldatos_pred$df$TAD, na.rm = TRUE))
xx <- cbind(model, rmse, nrmse)
RMSE <- rbind(RMSE, xx)
}
model = "FR~Prom.dB"
for (i in 1:length(ID)) {
newdata <- ldatos_pred_TRUE[(i*2 - 1):(i*2), row = T]
pred <- predict(func.gen.add.06, newdata)
rmse <- pred2meas(ldatos_pred_TRUE$df$FR[(i*2 - 1):(i*2)], pred, measure = "RMSE")
nrmse = sqrt(mean(((as.numeric(pred)) - (ldatos_pred_TRUE$df$FR[(i*2 - 1):(i*2)]))^{2})/var(ldatos_pred$df$FR, na.rm = TRUE))
xx <- cbind(model, rmse, nrmse)
RMSE <- rbind(RMSE, xx)
}
model = "SpO2~Prom.dB"
for (i in 1:length(ID)) {
newdata <- ldatos_pred_TRUE[(i*2 - 1):(i*2), row = T]
pred <- predict(func.gen.add.07, newdata)
rmse <- pred2meas(ldatos_pred_TRUE$df$SpO2[(i*2 - 1):(i*2)], pred, measure = "RMSE")
nrmse = sqrt(mean(((as.numeric(pred)) - (ldatos_pred_TRUE$df$SpO2[(i*2 - 1):(i*2)]))^{2})/var(ldatos_pred$df$SpO2, na.rm = TRUE))
xx <- cbind(model, rmse, nrmse)
RMSE <- rbind(RMSE, xx)
}
model = "FC~Luz1"
for (i in 1:length(ID)) {
newdata <- ldatos_pred_TRUE[(i*2 - 1):(i*2), row = T]
pred <- predict(func.gen.add.08, newdata)
rmse <- pred2meas(ldatos_pred_TRUE$df$FC[(i*2 - 1):(i*2)], pred, measure = "RMSE")
nrmse = sqrt(mean(((as.numeric(pred)) - (ldatos_pred_TRUE$df$FC[(i*2 - 1):(i*2)]))^{2})/var(ldatos_pred$df$FC, na.rm = TRUE))
xx <- cbind(model, rmse, nrmse)
RMSE <- rbind(RMSE, xx)
}
model = "FR~Luz1"
for (i in 1:length(ID)) {
newdata <- ldatos_pred_TRUE[(i*2 - 1):(i*2), row = T]
pred <- predict(func.gen.add.09, newdata)
rmse <- pred2meas(ldatos_pred_TRUE$df$FR[(i*2 - 1):(i*2)], pred, measure = "RMSE")
nrmse = sqrt(mean(((as.numeric(pred)) - (ldatos_pred_TRUE$df$FR[(i*2 - 1):(i*2)]))^{2})/var(ldatos_pred$df$FR, na.rm = TRUE))
xx <- cbind(model, rmse, nrmse)
RMSE <- rbind(RMSE, xx)
}
model = "SpO2~Luz2"
for (i in 1:length(ID)) {
newdata <- ldatos_pred_TRUE[(i*2 - 1):(i*2), row = T]
pred <- predict(func.gen.add.10, newdata)
rmse <- pred2meas(ldatos_pred_TRUE$df$SpO2[(i*2 - 1):(i*2)], pred, measure = "RMSE")
nrmse = sqrt(mean(((as.numeric(pred)) - (ldatos_pred_TRUE$df$SpO2[(i*2 - 1):(i*2)]))^{2})/var(ldatos_pred$df$SpO2, na.rm = TRUE))
xx <- cbind(model, rmse, nrmse)
RMSE <- rbind(RMSE, xx)
}
model = "TAD~Temp"
for (i in 1:length(ID)) {
newdata <- ldatos_pred_TRUE[(i*2 - 1):(i*2), row = T]
pred <- predict(func.gen.add.11, newdata)
rmse <- pred2meas(ldatos_pred_TRUE$df$TAD[(i*2 - 1):(i*2)], pred, measure = "RMSE")
nrmse = sqrt(mean(((as.numeric(pred)) - (ldatos_pred_TRUE$df$TAD[(i*2 - 1):(i*2)]))^{2})/var(ldatos_pred$df$TAD, na.rm = TRUE))
xx <- cbind(model, rmse, nrmse)
RMSE <- rbind(RMSE, xx)
}
model = "TAS~Luz2"
for (i in 1:length(ID)) {
newdata <- ldatos_pred_TRUE[(i*2 - 1):(i*2), row = T]
pred <- predict(func.gen.add.12, newdata)
rmse <- pred2meas(ldatos_pred_TRUE$df$TAS[(i*2 - 1):(i*2)], pred, measure = "RMSE")
nrmse = sqrt(mean(((as.numeric(pred)) - (ldatos_pred_TRUE$df$TAS[(i*2 - 1):(i*2)]))^{2})/var(ldatos_pred$df$TAS, na.rm = TRUE))
xx <- cbind(model, rmse, nrmse)
RMSE <- rbind(RMSE, xx)
}
model = "RASS~Luz2"
for (i in 1:length(ID)) {
newdata <- ldatos_pred_TRUE[(i*2 - 1):(i*2), row = T]
pred <- predict(func.gen.add.13, newdata)
rmse <- pred2meas(ldatos_pred_TRUE$df$RASS[(i*2 - 1):(i*2)], pred, measure = "RMSE")
nrmse = sqrt(mean(((as.numeric(pred)) - (ldatos_pred_TRUE$df$RASS[(i*2 - 1):(i*2)]))^{2})/var(ldatos_pred$df$RASS, na.rm = TRUE))
xx <- cbind(model, rmse, nrmse)
RMSE <- rbind(RMSE, xx)
}
model = "GCSm~Temp"
for (i in 1:length(ID)) {
newdata <- ldatos_pred_TRUE[(i*2 - 1):(i*2), row = T]
pred <- predict(func.gen.add.14, newdata)
rmse <- pred2meas(ldatos_pred_TRUE$df$GCSm[(i*2 - 1):(i*2)], pred, measure = "RMSE")
nrmse = sqrt(mean(((as.numeric(pred)) - (ldatos_pred_TRUE$df$GCSm[(i*2 - 1):(i*2)]))^{2})/var(ldatos_pred$df$GCSm, na.rm = TRUE))
xx <- cbind(model, rmse, nrmse)
RMSE <- rbind(RMSE, xx)
}
RMSE <- na.omit(RMSE)
RMSE$model <- as.factor(RMSE$model)
RMSE$rmse <- round(as.numeric(RMSE$rmse), 4)
RMSE$nrmse <- round(as.numeric(RMSE$nrmse), 4)
RMSE$nrmse[RMSE$nrmse > 10] <- mean(RMSE$nrmse)
RMSE$rmse[RMSE$model == "SpO2~Prom.dB" & RMSE$rmse > 50 | RMSE$model == "SpO2~Luz2" & RMSE$rmse > 50] <- mean(RMSE$rmse[RMSE$model == "SpO2~Prom.dB" | RMSE$model == "SpO2~Luz2"])# Gráficos Violines con Boxplot incorporado:
# Violin basico
RMSE <- RMSE %>%
mutate(model = factor(model, levels = c("FC~Prom.dB", "TAS~Prom.dB",
"GCSm~Prom.dB", "RASS~Prom.dB",
"TAD~Prom.dB", "FR~Prom.dB",
"SpO2~Prom.dB", "FC~Luz1",
"FR~Luz1", "SpO2~Luz2",
"TAD~Temp", "TAS~Luz2",
"RASS~Luz2", "GCSm~Temp"))) %>%
mutate(model = fct_reorder(model, desc(model))) # Para invertir el orden
violin <- ggplot(RMSE, aes(x = model, y = nrmse, fill = model)) +
geom_violin(trim = T) +
# Rotate
coord_flip() +
# Añadir boxplot
geom_boxplot(width = 0.1, fill = "white") +
theme_minimal(base_size = 15) +
theme(legend.position = "none") +
labs(title = "",x = "", y = "NRMSE")
ggsave(
"violin.png",
plot = violin,
scale = 1,
width = NA,
height = NA,
dpi = 300,
bg = "white"
)En primer lugar se muestra un violin plot que junta a todos los modelos propuestos (por orden decreciente de R-Sq) y muestra su comportamiento en funcion del NRMSE (Normalized Root Mean Squared Error) incluyendo en su interior un boxplot.
violin
violin.FC <- ggplot(RMSE[RMSE$model == "FC~Prom.dB" | RMSE$model == "FC~Luz1",], aes(x = model, y = rmse, fill = model)) +
geom_violin(trim = T, fill = 2) +
# Rotate
coord_flip() +
# Añadir boxplot
geom_boxplot(width = 0.1, fill = "white") +
theme_minimal(base_size = 15) +
theme(legend.position = "none") +
labs(title = "",x = "", y = "")
#
violin.TAS <- ggplot(RMSE[RMSE$model == "TAS~Prom.dB" | RMSE$model == "TAS~Luz2",], aes(x = model, y = rmse, fill = model)) +
geom_violin(trim = T, fill = 3) +
# Rotate
coord_flip() +
# Añadir boxplot
geom_boxplot(width = 0.1, fill = "white") +
theme_minimal(base_size = 15) +
theme(legend.position = "none") +
labs(title = "",x = "", y = "")
#
violin.TAD <- ggplot(RMSE[RMSE$model == "TAD~Prom.dB" | RMSE$model == "TAD~Temp",], aes(x = model, y = rmse, fill = model)) +
geom_violin(trim = T, fill = 4) +
# Rotate
coord_flip() +
# Añadir boxplot
geom_boxplot(width = 0.1, fill = "white") +
theme_minimal(base_size = 15) +
theme(legend.position = "none") +
labs(title = "",x = "", y = "")
#
violin.FR <- ggplot(RMSE[RMSE$model == "FR~Prom.dB" | RMSE$model == "FR~Luz1",], aes(x = model, y = rmse, fill = model)) +
geom_violin(trim = T, fill = 5) +
# Rotate
coord_flip() +
# Añadir boxplot
geom_boxplot(width = 0.1, fill = "white") +
theme_minimal(base_size = 15) +
theme(legend.position = "none") +
labs(title = "",x = "", y = "")
#
violin.SpO2 <- ggplot(RMSE[RMSE$model == "SpO2~Prom.dB" | RMSE$model == "SpO2~Luz2",], aes(x = model, y = rmse, fill = model)) +
geom_violin(trim = T, fill = 6) +
# Rotate
coord_flip() +
# Añadir boxplot
geom_boxplot(width = 0.1, fill = "white") +
theme_minimal(base_size = 15) +
theme(legend.position = "none") +
labs(title = "",x = "", y = "")
#
violin.RASS <- ggplot(RMSE[RMSE$model == "RASS~Prom.dB" | RMSE$model == "RASS~Luz2",], aes(x = model, y = rmse, fill = model)) +
geom_violin(trim = T, fill = 7) +
# Rotate
coord_flip() +
# Añadir boxplot
geom_boxplot(width = 0.1, fill = "white") +
theme_minimal(base_size = 15) +
theme(legend.position = "none") +
labs(title = "",x = "", y = "")
#
violin.GCSm <- ggplot(RMSE[RMSE$model == "GCSm~Prom.dB" | RMSE$model == "GCSm~Temp",], aes(x = model, y = rmse, fill = model)) +
geom_violin(trim = T, fill = 8) +
# Rotate
coord_flip() +
# Añadir boxplot
geom_boxplot(width = 0.1, fill = "white") +
theme_minimal(base_size = 15) +
theme(legend.position = "none") +
labs(title = "",x = "", y = "")
#El problema del NRMSE es que, pese a que por un lado convierte en comparables las predicciones de diversos modelos sobre diversas variables, también puede inducir a ciertos errores debido a la normalización. Es por ello que en la siguiente figura incluimos violin plots desglosados, dos por cada variable, para que así compartan la escala del eje X (RMSE) y poder valorar con detalle la capacidad predictiva.
Por ejemplo, no será lo mismo un error de 10 puntos en la Frecuencia Cardíaca, que tiene una \(SD=23.26\), que un error de 10 puntos en la saturación de oxígeno, que supondría la diferencia - prácticamente - entre la vida y la muerte.
violin.multi <- ggarrange(violin.FC, violin.TAS, violin.TAD, violin.FR, violin.SpO2,
violin.RASS, violin.GCSm,
font.label = list(size = 10, color = "black", face = "plain", family = NULL),
labels = c(paste0("FC ", "mean = ",
round(mean(ldatos_pred[["df"]]$FC),2), " +/- ",
round(sd(ldatos_pred[["df"]]$FC),2)),
paste0("TAS ", "mean = ",
round(mean(ldatos_pred[["df"]]$TAS),2), " +/- ",
round(sd(ldatos_pred[["df"]]$TAS),2)),
paste0("TAD ", "mean = ",
round(mean(ldatos_pred[["df"]]$TAD, na.rm = T),2), " +/- ",
round(sd(ldatos_pred[["df"]]$TAD),2)),
paste0("FR ", "mean = ",
round(mean(ldatos_pred[["df"]]$FR),2), " +/- ",
round(sd(ldatos_pred[["df"]]$FR),2)),
paste0("SpO2 ", "mean = ",
round(mean(ldatos_pred[["df"]]$SpO2,na.rm = T),2), " +/- ",
round(sd(ldatos_pred[["df"]]$SpO2),2)),
paste0("RASS ", "mean = ",
round(mean(ldatos_pred[["df"]]$RASS),2), " +/- ",
round(sd(ldatos_pred[["df"]]$RASS),2)),
paste0("GCSm ", "mean = ",
round(mean(ldatos_pred[["df"]]$GCSm),2), " +/- ",
round(sd(ldatos_pred[["df"]]$GCSm),2))),
ncol = 2, nrow = 4)
plot(violin.multi)
# ggsave(
# "violin_multi.png",
# plot = violin.multi,
# scale = 1,
# width = 30,
# height = 25,
# units = c("cm"),
# dpi = 300,
# bg = "white"
# )PARTE 6: CONCLUSIONES
Tras analizar estos datos podemos obtener una serie de conclusiones iniciales:
Algunas de las variables han sufrido muchos problemas de registro que deberían considerarse a la hora de tomar decisiones o, por lo menos, en futuros estudios. Especialmente la variable de Tª Ambiental muestra una gran cantidad de datos atípicos. Tras una limpieza minuciosa y exigente, consideramos que los resultados en base a esta variable pueden considerarse válidos, aunque originalmente su utilidad se vió muy comprometida.
La variable ambiental con el mayor impacto y más ubicuo a lo largo de las diferentes variables fisiológicas ha sido, con diferencia, el Promedio de Decibelios. Esta variable ha sido la que obtuvo mejores R-Sqared para todas las combinatorias con variables fisiológicas, por lo que se concluye que es la que más efecto está teniendo sobre la salud de los pacientes de UCI. Es interesante valorar (como se ve en los descriptivos generales) la naturaleza multi-modal de la distribución de esta variable, lo que indica alguna clase de patrón que podría ser corregido en las UCIs para minimizar su impacto.
En menor medida, también las variables relacionadas con la luz (Luz 1 y Luz 2) aparecen a menudo como buenas candidatas a ser evaluadas en las UCIs, aunque con R-Squared values muy inferiores a los que involucran al Promedio de Decibelios.
En cuanto al uso de técnicas de Datos Funcionales, se ve que en la inmensa mayoría de los casos estudiados los modelos son más precisos cuando se utilizan DF que con sus alternativas. Este beneficio es marginal pero también es constante, por lo que se puede concluir que efectivamente suponen una mejoría a la hora de modelar esta clase de datos.
En especial, son la Frecuencia Cardíaca y la Tensión Arterial Sistólica (que ya de por sí están muy correlacionadas) las variables fisiológicas que más afectadas se ven por la variable ambiental Promedio de Decibelios.
PARTE 7: NOTAS
¿Podría ser interesante testear también los GAM y sacar unos violinplot para comparar?
La correlación funcional sigue dando problemas y aparecen NAs de la nada. No parece comprometer el análisis ni aportar nada nuevo, sin embargo trataré de hacer que funcione para los artículos finales.
¿Sería conveniente un workflow del análisis para explicar de forma visual en un artículo?