0. Preparación de base de datos

Se cargaron y limpiaron los datos de tres sensores PurpleAir (CampMar_1, CampMar_2 y CampMar_3) y una base de referencia (Senamhi). Las principales acciones realizadas fueron:

  1. Homologación de nombres de variables: En la base de Senamhi, la columna CDM fue renombrada como pm25_sensor y date como fecha_hora para facilitar su uso.

  2. Selección de variables relevantes: En las tres bases de CampMar se conservaron únicamente las columnas necesarias: UTCDateTime, current_temp_f, current_humidity, pressure, pm2_5_cf_1, y pm2_5_cf_1_b.

  3. Conversión y filtrado de fechas

1. Promedio horario de Purple Air

RESULTADOS: Los datos brutos de PM₂.₅ fueron registrados a intervalos de 2 minutos. Para garantizar la calidad de los promedios horarios, se evaluó la cobertura dentro de cada hora y se aplicó un criterio de inclusión que requería un mínimo de 23 observaciones válidas por hora (equivalente al 75% de los 30 registros esperados). Las horas que no cumplieron con este umbral fueron excluidas del análisis antes del cálculo de los promedios horarios. En la estación CampMar_1, se eliminaron 7 horas (2.18% del total de 8,269 observaciones); en CampMar_2, se excluyeron 8 horas (2.50% de 8,388 observaciones); y en CampMar_3, se descartaron 11 horas, lo que representa el 4.66% de las 7,086 observaciones registradas.

############
# CampMar_1
############
# Asegurar que UTCDateTime esté bien formateado
CampMar_1 <- CampMar_1 %>%
  mutate(UTCDateTime = ymd_hms(gsub("T", " ", gsub("z", "", UTCDateTime))),
         fecha_hora = floor_date(UTCDateTime, unit = "hour"))

# Todas las horas presentes en la base
todas_las_horas <- CampMar_1 %>%
  distinct(fecha_hora) %>%
  pull(fecha_hora)

# Identificar las horas con al menos 23 mediciones válidas de pm2_5_cf_1
horas_validas <- CampMar_1 %>%
  filter(!is.na(pm2_5_cf_1)) %>%
  group_by(fecha_hora) %>%
  summarise(n_validas = n()) %>%
  filter(n_validas >= 23) %>%
  pull(fecha_hora)

# Filtrar la base original
CampMar_1_filtrado <- CampMar_1 %>%
  filter(fecha_hora %in% horas_validas)

# Calcular promedios horarios de todas las variables numéricas
CampMar_1_h <- CampMar_1_filtrado %>%
  group_by(fecha_hora) %>%
  summarise(across(where(is.numeric), ~ mean(.x, na.rm = TRUE)), .groups = "drop")

# Calcular y mostrar cuántas horas se eliminaron
horas_perdidas <- length(setdiff(todas_las_horas, horas_validas))
cat("Se eliminaron", horas_perdidas, "horas por no tener al menos 23 observaciones válidas de pm2_5_cf_1.\n")
## Se eliminaron 7 horas por no tener al menos 23 observaciones válidas de pm2_5_cf_1.
############
# CampMar_2
############
# Asegurar que UTCDateTime esté bien formateado
CampMar_2 <- CampMar_2 %>%
  mutate(UTCDateTime = ymd_hms(gsub("T", " ", gsub("z", "", UTCDateTime))),
         fecha_hora = floor_date(UTCDateTime, unit = "hour"))

# Todas las horas presentes en la base
todas_las_horas <- CampMar_2 %>%
  distinct(fecha_hora) %>%
  pull(fecha_hora)

# Identificar las horas con al menos 22 mediciones válidas de pm2_5_cf_1
horas_validas <- CampMar_2 %>%
  filter(!is.na(pm2_5_cf_1)) %>%
  group_by(fecha_hora) %>%
  summarise(n_validas = n()) %>%
  filter(n_validas >= 23) %>%
  pull(fecha_hora)

# Filtrar la base original
CampMar_2_filtrado <- CampMar_2 %>%
  filter(fecha_hora %in% horas_validas)

# Calcular promedios horarios de todas las variables numéricas
CampMar_2_h <- CampMar_2_filtrado %>%
  group_by(fecha_hora) %>%
  summarise(across(where(is.numeric), ~ mean(.x, na.rm = TRUE)), .groups = "drop")

# Calcular y mostrar cuántas horas se eliminaron
horas_perdidas <- length(setdiff(todas_las_horas, horas_validas))
cat("Se eliminaron", horas_perdidas, "horas por no tener al menos 23 observaciones válidas de pm2_5_cf_1.\n")
## Se eliminaron 8 horas por no tener al menos 23 observaciones válidas de pm2_5_cf_1.
############
# CampMar_3
############
# Asegurar que UTCDateTime esté bien formateado
CampMar_3 <- CampMar_3 %>%
  mutate(UTCDateTime = ymd_hms(gsub("T", " ", gsub("z", "", UTCDateTime))),
         fecha_hora = floor_date(UTCDateTime, unit = "hour"))

# Todas las horas presentes en la base
todas_las_horas <- CampMar_3 %>%
  distinct(fecha_hora) %>%
  pull(fecha_hora)

# Identificar las horas con al menos 22 mediciones válidas de pm2_5_cf_1
horas_validas <- CampMar_3 %>%
  filter(!is.na(pm2_5_cf_1)) %>%
  group_by(fecha_hora) %>%
  summarise(n_validas = n()) %>%
  filter(n_validas >= 23) %>%
  pull(fecha_hora)

# Filtrar la base original
CampMar_3_filtrado <- CampMar_3 %>%
  filter(fecha_hora %in% horas_validas)

# Calcular promedios horarios de todas las variables numéricas
CampMar_3_h <- CampMar_3_filtrado %>%
  group_by(fecha_hora) %>%
  summarise(across(where(is.numeric), ~ mean(.x, na.rm = TRUE)), .groups = "drop")

# Calcular y mostrar cuántas horas se eliminaron
horas_perdidas <- length(setdiff(todas_las_horas, horas_validas))
cat("Se eliminaron", horas_perdidas, "horas por no tener al menos 23 observaciones válidas de pm2_5_cf_1.\n")
## Se eliminaron 11 horas por no tener al menos 23 observaciones válidas de pm2_5_cf_1.

2. Busqueda de datos no plausibles T, H Y P - Creación de DP

RESULTADOS:

Se revisaron los valores registrados en tres variables ambientales clave: temperatura, presión atmosférica y humedad relativa. Para verificar su plausibilidad, se utilizaron los siguientes rangos de referencia: 0 a 55°C para la temperatura, 900 a 1200 hPa para la presión atmosférica y 0 a 100% para la humedad relativa. No se eliminaron datos en esta etapa. Adicionalmente, se calculó el punto de rocío (°C) a partir de la temperatura y la humedad relativa, utilizando la fórmula propuesta por Campbell Scientific.

## ---- CampMar_1_h ----
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   16.65   19.81   21.98   23.01   25.13   37.17
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   994.7   999.7  1001.1  1001.1  1002.4  1007.0
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   32.50   57.53   65.60   62.39   68.93   76.30
## 
## ---- CampMar_2_h ----
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   17.13   20.06   22.31   23.29   25.56   39.13
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   995.3  1000.3  1001.7  1001.6  1003.0  1007.5
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   31.33   57.53   65.20   62.14   68.37   76.10
## 
## ---- CampMar_3_h ----
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   16.31   18.93   21.19   22.27   24.63   37.17
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   995.2  1000.1  1001.4  1001.3  1002.6  1006.8
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   35.37   59.77   67.20   64.42   70.60   77.57

3. Evaluación de la coherencia entre canales A y B del sensor

Para evaluar la consistencia interna de los sensores PurpleAir, se compararon los valores horarios registrados por los canales A (pm2_5_cf_1) y B (pm2_5_cf_1_b) mediante el cálculo de la diferencia porcentual, utilizando la siguiente fórmula: diferencia horaria = ((A − B) × 2) / (A + B) × 100%. Se definieron dos criterios de calidad para detectar discrepancias sustanciales entre canales: (1) una diferencia porcentual superior al 70%, y (2) una diferencia absoluta mayor a 5 µg/m³. Cada hora fue clasificada como “no válida” si incumplía los dos criterios. Esta estrategia se aplicó de forma independiente a cada estación con el fin de identificar posibles fallos de medición o desviaciones técnicas relevantes.

RESULTADOS: CampMar_1_h: 16 horas no pasaron la validación CampMar_2_h: 991 horas no pasaron la validación CampMar_3_h: 16 horas no pasaron la validación CampMar_1_h: 0.19 % de las horas no pasaron la validación CampMar_2_h: 11.81 % de las horas no pasaron la validación CampMar_3_h: 0.23 % de las horas no pasaron la validación

# CampMar_1_h
CampMar_1_h <- CampMar_1_h %>%
  mutate(
    diferencia_pct = ((pm2_5_cf_1 - pm2_5_cf_1_b) * 2) / (pm2_5_cf_1 + pm2_5_cf_1_b) * 100,
    diferencia_abs = abs(pm2_5_cf_1 - pm2_5_cf_1_b),
    fallo_dif_pct = abs(diferencia_pct) > 70,
    fallo_dif_abs = diferencia_abs > 5,
    no_paso_validacion = ifelse(fallo_dif_pct & fallo_dif_abs, "no pasó", "ok")
  )

# CampMar_2_h
CampMar_2_h <- CampMar_2_h %>%
  mutate(
    diferencia_pct = ((pm2_5_cf_1 - pm2_5_cf_1_b) * 2) / (pm2_5_cf_1 + pm2_5_cf_1_b) * 100,
    diferencia_abs = abs(pm2_5_cf_1 - pm2_5_cf_1_b),
    fallo_dif_pct = abs(diferencia_pct) > 70,
    fallo_dif_abs = diferencia_abs > 5,
    no_paso_validacion = ifelse(fallo_dif_pct & fallo_dif_abs, "no pasó", "ok")
  )

# CampMar_3_h
CampMar_3_h <- CampMar_3_h %>%
  mutate(
    diferencia_pct = ((pm2_5_cf_1 - pm2_5_cf_1_b) * 2) / (pm2_5_cf_1 + pm2_5_cf_1_b) * 100,
    diferencia_abs = abs(pm2_5_cf_1 - pm2_5_cf_1_b),
    fallo_dif_pct = abs(diferencia_pct) > 70,
    fallo_dif_abs = diferencia_abs > 5,
    no_paso_validacion = ifelse(fallo_dif_pct & fallo_dif_abs, "no pasó", "ok")
  )

#############################
# % de datos perdidos 
############################
# CampMar_1_h
no_pasan_1 <- sum(CampMar_1_h$no_paso_validacion == "no pasó")
cat("CampMar_1_h: ", no_pasan_1, "horas no pasaron la validación\n")
## CampMar_1_h:  16 horas no pasaron la validación
# CampMar_2_h
no_pasan_2 <- sum(CampMar_2_h$no_paso_validacion == "no pasó")
cat("CampMar_2_h: ", no_pasan_2, "horas no pasaron la validación\n")
## CampMar_2_h:  991 horas no pasaron la validación
# CampMar_3_h
no_pasan_3 <- sum(CampMar_3_h$no_paso_validacion == "no pasó")
cat("CampMar_3_h: ", no_pasan_3, "horas no pasaron la validación\n")
## CampMar_3_h:  16 horas no pasaron la validación
# CampMar_1_h
porc_1 <- no_pasan_1 / nrow(CampMar_1_h) * 100
cat("CampMar_1_h: ", round(porc_1, 2), "% de las horas no pasaron la validación\n")
## CampMar_1_h:  0.19 % de las horas no pasaron la validación
# CampMar_2_h
porc_2 <- no_pasan_2 / nrow(CampMar_2_h) * 100
cat("CampMar_2_h: ", round(porc_2, 2), "% de las horas no pasaron la validación\n")
## CampMar_2_h:  11.81 % de las horas no pasaron la validación
# CampMar_3_h
porc_3 <- no_pasan_3 / nrow(CampMar_3_h) * 100
cat("CampMar_3_h: ", round(porc_3, 2), "% de las horas no pasaron la validación\n")
## CampMar_3_h:  0.23 % de las horas no pasaron la validación
#############################
# Filtrar datos 
############################

CampMar_1_h_valido <- CampMar_1_h %>%
  filter(no_paso_validacion != "no pasó")

CampMar_2_h_valido <- CampMar_2_h %>%
  filter(no_paso_validacion != "no pasó")

CampMar_3_h_valido <- CampMar_3_h %>%
  filter(no_paso_validacion != "no pasó")

Grafica para ver si el horario esta similar o no.

Se hace una revisión de la base de datos previa en tableau y se encuentra que a Senamhi hay que reducirle una hora y a Camp Mar 5 horas.

###########################################
# hacer graficas para ver horas  
###########################################

library(ggplot2)

library(lubridate)

Senamhi <- Senamhi %>%
  mutate(fecha_hora_1 = fecha_hora - hours(1)) 

CampMar_1_h_valido <- CampMar_1_h_valido %>%
  mutate(fecha_hora_1 = fecha_hora - hours(5)) 

CampMar_2_h_valido <- CampMar_2_h_valido %>%
  mutate(fecha_hora_1 = fecha_hora - hours(5)) 

CampMar_3_h_valido <- CampMar_3_h_valido %>%
  mutate(fecha_hora_1 = fecha_hora - hours(5)) 


# Filtrar datos del 31 de diciembre de 2024 al 1 de enero de 2025
CampMar_1_dec31 <- CampMar_1_h_valido %>%
  filter(fecha_hora_1 >= as.POSIXct("2024-12-31 00:00:00") &
         fecha_hora_1 <  as.POSIXct("2025-01-01 00:00:00"))

CampMar_2_dec31 <- CampMar_2_h_valido %>%
  filter(fecha_hora_1 >= as.POSIXct("2024-12-31 00:00:00") &
         fecha_hora_1 <  as.POSIXct("2025-01-01 00:00:00"))

CampMar_3_dec31 <- CampMar_3_h_valido %>%
  filter(fecha_hora_1 >= as.POSIXct("2024-12-31 00:00:00") &
         fecha_hora_1 <  as.POSIXct("2025-01-01 00:00:00"))

Senamhi_dec31 <- Senamhi %>%
  filter(fecha_hora_1 >= as.POSIXct("2024-12-31 00:00:00") &
         fecha_hora_1 <  as.POSIXct("2025-01-01 00:00:00"))

# Graficar
ggplot() +
  geom_line(data = CampMar_1_dec31, aes(x = fecha_hora_1, y = pm2_5_cf_1, color = "CampMar_1"), alpha = 0.8) +
  geom_line(data = CampMar_2_dec31, aes(x = fecha_hora_1, y = pm2_5_cf_1, color = "CampMar_2"), alpha = 0.8) +
  geom_line(data = CampMar_3_dec31, aes(x = fecha_hora_1, y = pm2_5_cf_1, color = "CampMar_3"), alpha = 0.8) +
  geom_line(data = Senamhi_dec31, aes(x = fecha_hora_1, y = pm25_sensor, color = "SENAMHI"), alpha = 0.8) +
  scale_color_manual(values = c("CampMar_1" = "blue", "CampMar_2" = "green",
                                "CampMar_3" = "orange", "SENAMHI" = "red")) +
  labs(title = "Fig 1. Concentración horaria de PM2.5",
       subtitle = "31 de diciembre de 2024",
       x = "Hora",
       y = "PM₂.₅ (µg/m³)",
       color = "Sensor") +
  scale_x_datetime(date_labels = "%H:%M", date_breaks = "2 hours") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

####
# ORDENAR BASE DE DATOS 
####


# CampMar1
CampMar1_Senamhi <- full_join(
  CampMar_1_h_valido,
  Senamhi,
  by = "fecha_hora_1"
)

names(CampMar1_Senamhi)
##  [1] "fecha_hora.x"       "current_temp_f"     "current_humidity"  
##  [4] "pressure"           "pm2_5_cf_1"         "pm2_5_cf_1_b"      
##  [7] "current_temp_c"     "dew_point"          "diferencia_pct"    
## [10] "diferencia_abs"     "fallo_dif_pct"      "fallo_dif_abs"     
## [13] "no_paso_validacion" "fecha_hora_1"       "fecha_hora.y"      
## [16] "pm25_sensor"
# Cambiar nombre 
CampMar1_Senamhi <- CampMar1_Senamhi %>%
  rename(pm25_referencia = pm25_sensor)
# promedio de pm2.5 
CampMar1_Senamhi <- CampMar1_Senamhi %>%
  mutate(pm25_sensor = (pm2_5_cf_1 + pm2_5_cf_1_b) / 2)
# quedarme con las que quiero 
CampMar1_Senamhi <- CampMar1_Senamhi %>%
  select(fecha_hora_1, pm25_sensor, pm2_5_cf_1, pm2_5_cf_1_b,
         pm25_referencia, current_temp_c, dew_point, current_humidity, pressure)

CampMar1_Senamhi <- CampMar1_Senamhi %>%
  rename(
    fecha_hora = fecha_hora_1,
    temperature = current_temp_c,
    humedad = current_humidity,
    presion = pressure
  )

CampMar1_Senamhi <- CampMar1_Senamhi %>%
  rename(temperatura = temperature) 


# CampMar2 
CampMar2_Senamhi <- full_join(
  CampMar_2_h_valido,
  Senamhi,
  by = "fecha_hora_1"
)

names(CampMar2_Senamhi)
##  [1] "fecha_hora.x"       "current_temp_f"     "current_humidity"  
##  [4] "pressure"           "pm2_5_cf_1"         "pm2_5_cf_1_b"      
##  [7] "current_temp_c"     "dew_point"          "diferencia_pct"    
## [10] "diferencia_abs"     "fallo_dif_pct"      "fallo_dif_abs"     
## [13] "no_paso_validacion" "fecha_hora_1"       "fecha_hora.y"      
## [16] "pm25_sensor"
# Cambiar nombre 
CampMar2_Senamhi <- CampMar2_Senamhi %>%
  rename(pm25_referencia = pm25_sensor)
# promedio de pm2.5 
CampMar2_Senamhi <- CampMar2_Senamhi %>%
  mutate(pm25_sensor = (pm2_5_cf_1 + pm2_5_cf_1_b) / 2)
# quedarme con las que quiero 
CampMar2_Senamhi <- CampMar2_Senamhi %>%
  select(fecha_hora_1, pm25_sensor, pm2_5_cf_1, pm2_5_cf_1_b,
         pm25_referencia, current_temp_c, dew_point, current_humidity, pressure)

CampMar2_Senamhi <- CampMar2_Senamhi %>%
  rename(
    fecha_hora = fecha_hora_1,
    temperature = current_temp_c,
    humedad = current_humidity,
    presion = pressure
  )

CampMar2_Senamhi <- CampMar2_Senamhi %>%
  rename(temperatura = temperature) 

# CampMar3 
CampMar3_Senamhi <- full_join(
  CampMar_3_h_valido,
  Senamhi,
  by = "fecha_hora_1"
)

names(CampMar3_Senamhi)
##  [1] "fecha_hora.x"       "current_temp_f"     "current_humidity"  
##  [4] "pressure"           "pm2_5_cf_1"         "pm2_5_cf_1_b"      
##  [7] "current_temp_c"     "dew_point"          "diferencia_pct"    
## [10] "diferencia_abs"     "fallo_dif_pct"      "fallo_dif_abs"     
## [13] "no_paso_validacion" "fecha_hora_1"       "fecha_hora.y"      
## [16] "pm25_sensor"
# Cambiar nombre 
CampMar3_Senamhi <- CampMar3_Senamhi %>%
  rename(pm25_referencia = pm25_sensor)
# promedio de pm2.5 
CampMar3_Senamhi <- CampMar3_Senamhi %>%
  mutate(pm25_sensor = (pm2_5_cf_1 + pm2_5_cf_1_b) / 2)
# quedarme con las que quiero 
CampMar3_Senamhi <- CampMar3_Senamhi %>%
  select(fecha_hora_1, pm25_sensor, pm2_5_cf_1, pm2_5_cf_1_b,
         pm25_referencia, current_temp_c, dew_point, current_humidity, pressure)

CampMar3_Senamhi <- CampMar3_Senamhi %>%
  rename(
    fecha_hora = fecha_hora_1,
    temperature = current_temp_c,
    humedad = current_humidity,
    presion = pressure
  )

CampMar3_Senamhi <- CampMar3_Senamhi %>%
  rename(temperatura = temperature) 

# Crear variable de fecha y hora 
CampMar1_Senamhi <- CampMar1_Senamhi %>%
  mutate(
    fecha = as.Date(fecha_hora),
    hora = format(fecha_hora, "%H:%M:%S")
  )
CampMar2_Senamhi <- CampMar2_Senamhi %>%
  mutate(
    fecha = as.Date(fecha_hora),
    hora = format(fecha_hora, "%H:%M:%S")
  )
CampMar3_Senamhi <- CampMar3_Senamhi %>%
  mutate(
    fecha = as.Date(fecha_hora),
    hora = format(fecha_hora, "%H:%M:%S")
  )


#### Borar datos faltantes 
CampMar1_Senamhi <- CampMar1_Senamhi %>%
  filter(!(is.na(pm25_sensor) & is.na(pm25_referencia)))

CampMar2_Senamhi <- CampMar2_Senamhi %>%
  filter(!(is.na(pm25_sensor) & is.na(pm25_referencia)))

CampMar3_Senamhi <- CampMar3_Senamhi %>%
  filter(!(is.na(pm25_sensor) & is.na(pm25_referencia)))

# solo se tienen datos de senamhi hasta agosto 
CampMar1_Senamhi <- CampMar1_Senamhi %>%
  filter(fecha_hora <= as.POSIXct("2025-04-30 23:59:59"))
CampMar2_Senamhi <- CampMar2_Senamhi %>%
  filter(fecha_hora <= as.POSIXct("2025-04-30 23:59:59"))
CampMar3_Senamhi <- CampMar3_Senamhi %>%
  filter(fecha_hora <= as.POSIXct("2025-04-30 23:59:59"))

## desde esta fecha se puso el monitor 

CampMar1_Senamhi <- CampMar1_Senamhi %>%
  filter(fecha_hora > as.POSIXct("2024-06-18 14:00:00"))

CampMar2_Senamhi <- CampMar2_Senamhi %>%
  filter(fecha_hora > as.POSIXct("2024-06-18 14:00:00"))

CampMar3_Senamhi <- CampMar3_Senamhi %>%
  filter(fecha_hora > as.POSIXct("2024-06-18 14:00:00"))

Grafica para ver meses con datos faltantes, pensando LOWOC

El monitor CampMar_2 presentó una pérdida de datos a partir del mes de febrero de 2025, con un 88% de registros faltantes en ese mes. Por otro lado, el monitor CampMar_3 mostró una disminución en la cobertura desde marzo de 2025, con 46.5% de datos perdidos ese mes. Dado que el procedimiento de validación cruzada Leave-One-Week-Out requiere continuidad temporal en los datos, se decide excluir estos meses con alta proporción de vacíos.

library(dplyr)
library(lubridate)
library(tidyr)
library(ggplot2)


# Función para contar datos faltantes y calcular porcentaje por mes
contar_faltantes <- function(base, nombre_sensor) {
  base %>%
    mutate(
      mes = floor_date(fecha_hora, "month"),
      dias_en_mes = days_in_month(fecha_hora),
      datos_esperados = dias_en_mes * 24
    ) %>%
    group_by(mes) %>%
    summarise(
      faltantes_sensor = sum(is.na(pm25_sensor)),
      faltantes_ref = sum(is.na(pm25_referencia)),
      datos_esperados = first(datos_esperados)
    ) %>%
    pivot_longer(cols = starts_with("faltantes"), names_to = "tipo", values_to = "faltantes") %>%
    mutate(
      sensor = nombre_sensor,
      porcentaje = round((faltantes / datos_esperados) * 100, 1),
      etiqueta = ifelse(porcentaje > 12, paste0(porcentaje, "%"), NA)
    )
}

# Aplicar a cada base
faltantes_1 <- contar_faltantes(CampMar1_Senamhi, "CampMar_1")
faltantes_2 <- contar_faltantes(CampMar2_Senamhi, "CampMar_2")
faltantes_3 <- contar_faltantes(CampMar3_Senamhi, "CampMar_3")

# Combinar resultados
faltantes_total <- bind_rows(faltantes_1, faltantes_2, faltantes_3)

# Renombrar etiquetas
faltantes_total$tipo <- recode(faltantes_total$tipo,
                               "faltantes_sensor" = "Sensor",
                               "faltantes_ref" = "Referencia")

faltantes_total$mes <- as.Date(faltantes_total$mes)

# Gráfico
ggplot(faltantes_total, aes(x = mes, y = faltantes, fill = tipo)) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_text(aes(label = etiqueta),
            position = position_dodge(width = 25),
            vjust = -0.3, size = 3, na.rm = TRUE) +
  facet_wrap(~ sensor, ncol = 1) +
  labs(title = "Fig 2. Datos faltantes por mes y tipo",
       x = "Mes",
       y = "Cantidad de datos faltantes",
       fill = "Tipo de dato") +
  scale_x_date(date_labels = "%b %Y", date_breaks = "1 month") +
  ylim(0, 900) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

CampMar2_Senamhi <- CampMar2_Senamhi %>%
  filter(fecha_hora < ymd("2025-02-01"))

CampMar3_Senamhi <- CampMar3_Senamhi %>%
  filter(fecha_hora < ymd("2025-03-01"))

# Función para contar datos faltantes por mes para cada base
contar_faltantes <- function(base, nombre_sensor) {
  base %>%
    mutate(mes = floor_date(fecha_hora, "month")) %>%
    group_by(mes) %>%
    summarise(
      faltantes_sensor = sum(is.na(pm25_sensor)),
      faltantes_ref = sum(is.na(pm25_referencia))
    ) %>%
    pivot_longer(cols = starts_with("faltantes"), names_to = "tipo", values_to = "faltantes") %>%
    mutate(sensor = nombre_sensor)
}

# Aplicar la función a cada base ya depurada
faltantes_1 <- contar_faltantes(CampMar1_Senamhi, "CampMar_1")
faltantes_2 <- contar_faltantes(CampMar2_Senamhi, "CampMar_2")
faltantes_3 <- contar_faltantes(CampMar3_Senamhi, "CampMar_3")

# Combinar resultados
faltantes_total <- bind_rows(faltantes_1, faltantes_2, faltantes_3)

# Etiquetas limpias
faltantes_total$tipo <- recode(faltantes_total$tipo,
                               "faltantes_sensor" = "Sensor",
                               "faltantes_ref" = "Referencia")

# Convertir mes a Date
faltantes_total$mes <- as.Date(faltantes_total$mes)

# Graficar sin % ni etiquetas
ggplot(faltantes_total, aes(x = mes, y = faltantes, fill = tipo)) +
  geom_bar(stat = "identity", position = "dodge") +
  facet_wrap(~ sensor, ncol = 1) +
  labs(title = "Fig 3. Datos faltantes por mes y tipo",
       x = "Mes",
       y = "Cantidad de datos faltantes",
       fill = "Tipo de dato") +
  scale_x_date(date_labels = "%b %Y", date_breaks = "1 month") +
  scale_y_continuous(limits = c(0, 700)) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

## ---- CampMar_1 ----
## Total registros: 7611 
## Faltantes Sensor: 173  ( 2.27 %)
## Faltantes Referencia: 141  ( 1.85 %)
## ---- CampMar_2 ----
## Total registros: 5469 
## Faltantes Sensor: 30  ( 0.55 %)
## Faltantes Referencia: 46  ( 0.84 %)
## ---- CampMar_3 ----
## Total registros: 6141 
## Faltantes Sensor: 50  ( 0.81 %)
## Faltantes Referencia: 134  ( 2.18 %)

Estacionalidad

4. Definir los modelos que probare

Se definieron un total de 24 modelos lineales con el objetivo de explorar diferentes combinaciones de variables ambientales que puedan contribuir a la calibración de los valores de PM₂.₅ medidos por sensores PurpleAir frente a una referencia oficial. En todos los modelos, la variable dependiente fue pm25_referencia, correspondiente a la medición del instrumento de referencia, mientras que la variable independiente base fue pm25_sensor, correspondiente al canal A del sensor PurpleAir.

A partir de este modelo base, se incluyeron progresivamente variables ambientales (temperatura, humedad relativa, punto de rocío y presión atmosférica), así como sus interacciones y combinaciones no lineales.

lista_modelos <- list(
  modelo_1  = pm25_referencia ~ pm25_sensor,
  modelo_2  = pm25_referencia ~ pm25_sensor + temperatura,
  modelo_3  = pm25_referencia ~ pm25_sensor + humedad,
  modelo_4  = pm25_referencia ~ pm25_sensor + dew_point,
  modelo_5  = pm25_referencia ~ pm25_sensor + presion,
  modelo_6  = pm25_referencia ~ pm25_sensor + temperatura + humedad,
  modelo_7  = pm25_referencia ~ pm25_sensor + temperatura * humedad,
  modelo_8  = pm25_referencia ~ pm25_sensor + temperatura * dew_point,
  modelo_9  = pm25_referencia ~ pm25_sensor + temperatura * presion,
  modelo_10 = pm25_referencia ~ pm25_sensor + humedad * dew_point,
  modelo_11 = pm25_referencia ~ pm25_sensor + humedad * presion,
  modelo_12 = pm25_referencia ~ pm25_sensor + dew_point * presion,
  modelo_13 = pm25_referencia ~ pm25_sensor * temperatura,
  modelo_14 = pm25_referencia ~ pm25_sensor * humedad,
  modelo_15 = pm25_referencia ~ pm25_sensor * dew_point,
  modelo_16 = pm25_referencia ~ pm25_sensor * presion,
  modelo_17 = pm25_referencia ~ pm25_sensor + humedad * temperatura + dew_point,
  modelo_18 = pm25_referencia ~ pm25_sensor + humedad * temperatura + presion,
  modelo_19 = pm25_referencia ~ pm25_sensor + humedad * dew_point + presion,
  modelo_20 = pm25_referencia ~ pm25_sensor + dew_point * presion + dew_point,
  modelo_con_estacion = pm25_referencia ~ pm25_sensor + temperatura * dew_point + epoca_humedad,
  modelo_22 = pm25_referencia ~ pm25_sensor + temperatura * humedad + epoca_humedad,
  modelo_23 = pm25_referencia ~ pm25_sensor + humedad * presion + epoca_humedad,
  modelo_24 = pm25_referencia ~ pm25_sensor + temperatura + humedad + epoca_humedad,
  modelo_25 = pm25_referencia ~ pm25_sensor + temperatura * presion + epoca_humedad,
  modelo_26 = pm25_referencia ~ pm25_sensor * temperatura + epoca_humedad,
  modelo_27 = pm25_referencia ~ pm25_sensor * humedad + epoca_humedad
)

# Cargar tibble
library(tibble)

# Crear tabla con nombres de modelos y fórmulas
tabla_modelos <- tibble(
  Modelo = names(lista_modelos),
  Formula = sapply(lista_modelos, function(f) deparse(f))
)

# Mostrar la tabla en el notebook
tabla_modelos
## # A tibble: 27 × 2
##    Modelo    Formula                                                
##    <chr>     <chr>                                                  
##  1 modelo_1  pm25_referencia ~ pm25_sensor                          
##  2 modelo_2  pm25_referencia ~ pm25_sensor + temperatura            
##  3 modelo_3  pm25_referencia ~ pm25_sensor + humedad                
##  4 modelo_4  pm25_referencia ~ pm25_sensor + dew_point              
##  5 modelo_5  pm25_referencia ~ pm25_sensor + presion                
##  6 modelo_6  pm25_referencia ~ pm25_sensor + temperatura + humedad  
##  7 modelo_7  pm25_referencia ~ pm25_sensor + temperatura * humedad  
##  8 modelo_8  pm25_referencia ~ pm25_sensor + temperatura * dew_point
##  9 modelo_9  pm25_referencia ~ pm25_sensor + temperatura * presion  
## 10 modelo_10 pm25_referencia ~ pm25_sensor + humedad * dew_point    
## # ℹ 17 more rows

5. Proceso de Leave-One-Week-Out Cross

Este enfoque consiste en dividir el conjunto de datos por semanas. En cada iteración, se selecciona una semana específica (semana s) como conjunto de prueba, mientras que todas las demás semanas se utilizan como conjunto de entrenamiento.

El procedimiento sigue estos pasos: 1. Selección del conjunto de prueba: Se deja fuera una semana (semana s) del conjunto de datos para usarla como conjunto de prueba. 2. Entrenamiento del modelo: El modelo se entrena utilizando los datos de las semanas restantes (todas excepto la semana s). 3. Predicción y evaluación: El modelo entrenado realiza predicciones sobre los datos de la semana s (conjunto de prueba), y se calculan los errores de predicción (como RMSE, MAE, Bias, R²). 4. Repetición del proceso: Se repite este procedimiento dejando fuera una semana diferente en cada iteración, hasta que todas las semanas hayan sido utilizadas una vez como conjunto de prueba.

library(Metrics)

# SOLO PARA CAMPO DE MARTE 1 

# Crear variable de semana en la base
CampMar1_Senamhi <- CampMar1_Senamhi %>%
  mutate(id_semana = isoweek(fecha_hora))

# Inicializar resultados
semanas <- unique(CampMar1_Senamhi$id_semana)
resultados <- data.frame()

# Bucle LOWO-CV
for (modelo_nombre in names(lista_modelos)) {
  formula_actual <- lista_modelos[[modelo_nombre]]
  
  for (s in semanas) {
    datos_prueba <- subset(CampMar1_Senamhi, id_semana == s)
    datos_entrenamiento <- subset(CampMar1_Senamhi, id_semana != s)
    
    modelo <- lm(formula_actual, data = datos_entrenamiento)
    pred <- predict(modelo, newdata = datos_prueba)
    real <- datos_prueba$pm25_referencia
    
    metrica <- data.frame(
      modelo = modelo_nombre,
      semana = s,
      RMSE = rmse(real, pred),
      MAE = mae(real, pred),
      Bias = mean(pred - real),
      R2 = 1 - sum((pred - real)^2) / sum((real - mean(real))^2)
    )
    
    resultados <- rbind(resultados, metrica)
  }
}
# SOLO PARA CAMPO DE MARTE 2 

# Crear variable de semana en la base
CampMar2_Senamhi <- CampMar2_Senamhi %>%
  mutate(id_semana = isoweek(fecha_hora))

# Inicializar resultados
semanas <- unique(CampMar2_Senamhi$id_semana)
resultados2 <- data.frame()

# Bucle LOWO-CV
for (modelo_nombre in names(lista_modelos)) {
  formula_actual2 <- lista_modelos[[modelo_nombre]]
  
  for (s in semanas) {
    datos_prueba_2 <- subset(CampMar2_Senamhi, id_semana == s)
    datos_entrenamiento_2 <- subset(CampMar2_Senamhi, id_semana != s)
    
    modelo2 <- lm(formula_actual2, data = datos_entrenamiento_2)
    pred2 <- predict(modelo2, newdata = datos_prueba_2)
    real2 <- datos_prueba_2$pm25_referencia
    
    metrica2 <- data.frame(
      modelo2 = modelo_nombre,
      semana = s,
      RMSE = rmse(real2, pred2),
      MAE = mae(real2, pred2),
      Bias = mean(pred2 - real2),
      R2 = 1 - sum((pred2 - real2)^2) / sum((real2 - mean(real2))^2)
    )
    
    resultados2 <- rbind(resultados2, metrica2)
  }
}
# SOLA CAMPO DE MARTE 3 

# Crear variable de semana
CampMar3_Senamhi <- CampMar3_Senamhi %>%
  mutate(id_semana = isoweek(fecha_hora))

# Inicializar resultados
semanas_3 <- unique(CampMar3_Senamhi$id_semana)
resultados3 <- data.frame()

# Bucle LOWO-CV
for (modelo_nombre_3 in names(lista_modelos)) {
  formula_actual_3 <- lista_modelos[[modelo_nombre_3]]
  
  for (s in semanas_3) {
    datos_prueba_3 <- subset(CampMar3_Senamhi, id_semana == s)
    datos_entrenamiento_3 <- subset(CampMar3_Senamhi, id_semana != s)
    
    modelo_3 <- lm(formula_actual_3, data = datos_entrenamiento_3)
    pred_3 <- predict(modelo_3, newdata = datos_prueba_3)
    real_3 <- datos_prueba_3$pm25_referencia
    
    metrica_3 <- data.frame(
      modelo3 = modelo_nombre_3,
      semana = s,
      RMSE = rmse(real_3, pred_3),
      MAE = mae(real_3, pred_3),
      Bias = mean(pred_3 - real_3),
      R2 = 1 - sum((pred_3 - real_3)^2) / sum((real_3 - mean(real_3))^2)
    )
    
    resultados3 <- rbind(resultados3, metrica_3)
  }
}

6. Criterios de evaluación del modelo:

Criterios de evaluación del modelo:

Los resultados obtenidos a partir del proceso de validación cruzada Leave-One-Week-Out se evaluaron utilizando los siguientes indicadores:

•   Coeficiente de determinación (R²): Se consideró aceptable un R² igual o superior a 0.75, lo cual indica una correlación casi perfecta entre los datos calibrados del sensor PurpleAir y los valores del monitor de referencia.

•   Errores de predicción (RMSE, MAE y MAPE): Se buscaron valores lo más bajos posible para estos indicadores, ya que reflejan una mayor precisión del modelo al estimar los niveles reales de PM₂.₅.

•   Sesgo (Bias): Se aceptó un sesgo dentro del rango de –5 a 5 μg/m³. Valores fuera de este intervalo indican una sobreestimación o subestimación sistemática que compromete la validez del modelo.

En caso de que varios modelos cumplieran con estos criterios, se priorizó la selección de aquel que maximizara el R² y minimizara simultáneamente el RMSE y el MAE.

# Calcular métricas promedio por modelo para Campo de Marte 1
resumen_modelos_1 <- aggregate(cbind(RMSE, MAE, Bias, R2) ~ modelo, data = resultados, mean)

# Ordenar por R2 descendente, luego por RMSE y MAE ascendentes
resumen_modelos_ordenado_1 <- resumen_modelos_1 %>%
  arrange(desc(R2), RMSE, MAE)

# Mostrar resultados ordenados
print(resumen_modelos_ordenado_1)
##                 modelo     RMSE      MAE       Bias        R2
## 1            modelo_11 2.871283 2.168658 -0.4650165 0.8729463
## 2            modelo_22 2.873077 2.167780 -0.3762396 0.8717648
## 3             modelo_7 2.877728 2.180512 -0.3477333 0.8716420
## 4            modelo_17 2.879048 2.182674 -0.3483576 0.8716020
## 5             modelo_3 2.887125 2.177250 -0.4967850 0.8713737
## 6            modelo_14 2.888826 2.177807 -0.4982354 0.8713645
## 7            modelo_24 2.895104 2.182209 -0.3711624 0.8707383
## 8            modelo_18 2.888468 2.184962 -0.3582503 0.8698969
## 9            modelo_23 2.893437 2.192640 -0.4317115 0.8698837
## 10            modelo_6 2.916035 2.214816 -0.3336784 0.8692151
## 11           modelo_10 2.896114 2.192493 -0.3481613 0.8691885
## 12           modelo_27 2.911914 2.209681 -0.4405717 0.8686570
## 13            modelo_9 2.907120 2.184982 -0.4672534 0.8679352
## 14            modelo_8 2.853518 2.157507 -0.3306374 0.8679093
## 15           modelo_19 2.906893 2.197400 -0.3578301 0.8675162
## 16 modelo_con_estacion 2.858469 2.161562 -0.3348895 0.8673083
## 17            modelo_5 2.940235 2.215331 -0.4635012 0.8672787
## 18            modelo_1 2.930191 2.202046 -0.4839195 0.8671598
## 19           modelo_25 2.911290 2.195589 -0.4313118 0.8671200
## 20           modelo_16 2.949366 2.221190 -0.4686736 0.8668162
## 21           modelo_12 2.940796 2.204726 -0.3666454 0.8656705
## 22           modelo_20 2.940796 2.204726 -0.3666454 0.8656705
## 23            modelo_2 2.948089 2.215346 -0.4677738 0.8652353
## 24           modelo_13 2.972243 2.231184 -0.4687531 0.8641203
## 25           modelo_26 2.975086 2.242828 -0.4487636 0.8634398
## 26            modelo_4 2.982454 2.249593 -0.3604781 0.8625874
## 27           modelo_15 3.005881 2.264486 -0.3587895 0.8616852
# Calcular métricas promedio por modelo para Campo de Marte 2
resumen_modelos_2 <- aggregate(cbind(RMSE, MAE, Bias, R2) ~ modelo2, data = resultados2, mean)

# Ordenar por R2 descendente, luego por RMSE y MAE ascendentes
resumen_modelos_ordenado_2 <- resumen_modelos_2 %>%
  arrange(desc(R2), RMSE, MAE)

# Mostrar resultados ordenados
print(resumen_modelos_ordenado_2)
##                modelo2     RMSE      MAE       Bias        R2
## 1             modelo_8 3.021703 2.119506 -0.3106278 0.8601177
## 2  modelo_con_estacion 3.027684 2.125113 -0.3081807 0.8593338
## 3            modelo_18 3.051548 2.139912 -0.4277438 0.8521240
## 4            modelo_17 3.046906 2.138206 -0.4154127 0.8518990
## 5             modelo_7 3.047427 2.141282 -0.4141086 0.8517471
## 6            modelo_10 3.058891 2.146732 -0.4137933 0.8509974
## 7            modelo_22 3.053145 2.145807 -0.4206595 0.8509789
## 8            modelo_19 3.064619 2.147551 -0.4261553 0.8509493
## 9            modelo_24 3.097058 2.181979 -0.4316055 0.8492248
## 10            modelo_1 3.133786 2.216812 -0.4036962 0.8469882
## 11            modelo_6 3.108331 2.200986 -0.4246341 0.8465834
## 12            modelo_5 3.138034 2.215542 -0.4127500 0.8465499
## 13           modelo_16 3.153676 2.223208 -0.4332269 0.8459182
## 14            modelo_3 3.097881 2.192840 -0.4124900 0.8448102
## 15           modelo_14 3.107062 2.193363 -0.4218957 0.8443839
## 16            modelo_4 3.153677 2.233014 -0.3958596 0.8439925
## 17           modelo_27 3.130008 2.215263 -0.4305103 0.8412392
## 18            modelo_2 3.132455 2.219444 -0.3860640 0.8405003
## 19           modelo_11 3.100732 2.183417 -0.4141719 0.8404942
## 20           modelo_23 3.121415 2.202049 -0.4252255 0.8383206
## 21           modelo_13 3.162205 2.215065 -0.3846781 0.8348139
## 22           modelo_12 3.166616 2.237154 -0.3365209 0.8344003
## 23           modelo_20 3.166616 2.237154 -0.3365209 0.8344003
## 24           modelo_26 3.180316 2.236997 -0.4031840 0.8321132
## 25           modelo_15 3.237443 2.238927 -0.3931561 0.8273555
## 26            modelo_9 3.149446 2.227217 -0.3378373 0.8245153
## 27           modelo_25 3.163600 2.242934 -0.3568522 0.8231056
# Calcular métricas promedio por modelo para Campo de Marte 3
resumen_modelos_3 <- aggregate(cbind(RMSE, MAE, Bias, R2) ~ modelo3, data = resultados3, mean)

# Ordenar por R2 descendente, luego por RMSE y MAE ascendentes
resumen_modelos_ordenado_3 <- resumen_modelos_3 %>%
  arrange(desc(R2), RMSE, MAE)


# Mostrar resultados ordenados
print(resumen_modelos_ordenado_3)
##                modelo3     RMSE      MAE       Bias        R2
## 1  modelo_con_estacion 2.670585 2.032783 -0.3532635 0.8929898
## 2             modelo_8 2.675570 2.037205 -0.3560837 0.8923137
## 3            modelo_22 2.717883 2.051122 -0.5284363 0.8874300
## 4            modelo_24 2.768805 2.087608 -0.4846241 0.8861961
## 5            modelo_17 2.757322 2.093494 -0.5604026 0.8782078
## 6            modelo_10 2.767119 2.101827 -0.5485990 0.8776263
## 7             modelo_7 2.760958 2.097383 -0.5647030 0.8775899
## 8            modelo_19 2.775553 2.108810 -0.5574130 0.8763712
## 9            modelo_18 2.770150 2.104351 -0.5744561 0.8762753
## 10           modelo_12 2.855747 2.157038 -0.4454490 0.8748637
## 11           modelo_20 2.855747 2.157038 -0.4454490 0.8748637
## 12            modelo_9 2.813645 2.133281 -0.5761677 0.8741919
## 13           modelo_11 2.806646 2.145628 -0.6841325 0.8714762
## 14           modelo_14 2.846018 2.167845 -0.7511578 0.8700880
## 15            modelo_3 2.849910 2.177532 -0.7450740 0.8695073
## 16            modelo_6 2.856649 2.182618 -0.5369670 0.8691542
## 17           modelo_25 2.844187 2.164981 -0.6008829 0.8683198
## 18            modelo_1 2.881202 2.200569 -0.7211397 0.8667404
## 19            modelo_5 2.891580 2.222136 -0.6665048 0.8648465
## 20           modelo_15 2.907559 2.224491 -0.5257804 0.8648399
## 21           modelo_13 2.899770 2.213142 -0.6458618 0.8644763
## 22           modelo_23 2.846970 2.183318 -0.7014933 0.8642756
## 23           modelo_16 2.902880 2.226283 -0.6708787 0.8641603
## 24            modelo_4 2.925058 2.237367 -0.5453549 0.8633718
## 25            modelo_2 2.919584 2.230061 -0.6538289 0.8631514
## 26           modelo_26 2.919661 2.233600 -0.6735134 0.8609457
## 27           modelo_27 2.896811 2.221957 -0.7472537 0.8602383

Selección de modelo 1

Camp1 = modelo_11 = pm25_referencia ~ pm25_sensor + humedad * presion

# Camp1 - modelo_11 = pm25_referencia ~ pm25_sensor + humedad * presion,
modelo_final_1 <- lm(pm25_referencia ~ pm25_sensor + humedad * presion, data = CampMar1_Senamhi)
summary(modelo_final_1)
## 
## Call:
## lm(formula = pm25_referencia ~ pm25_sensor + humedad * presion, 
##     data = CampMar1_Senamhi)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -61.224  -1.657  -0.170   1.433  31.750 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -9.930e+02  1.360e+02  -7.301 3.15e-13 ***
## pm25_sensor      6.069e-01  2.241e-03 270.765  < 2e-16 ***
## humedad          1.420e+01  2.171e+00   6.542 6.46e-11 ***
## presion          9.995e-01  1.360e-01   7.351 2.19e-13 ***
## humedad:presion -1.424e-02  2.171e-03  -6.561 5.71e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.21 on 7292 degrees of freedom
##   (314 observations deleted due to missingness)
## Multiple R-squared:  0.9242, Adjusted R-squared:  0.9242 
## F-statistic: 2.223e+04 on 4 and 7292 DF,  p-value: < 2.2e-16
CampMar1_Senamhi <- CampMar1_Senamhi %>%
  mutate(pm25_ajustado = 
           -993.0 +
           0.6069 * pm25_sensor +
           14.20 * humedad +
           0.9995 * presion +
           (-0.01424) * humedad * presion)


ggplot(CampMar1_Senamhi, aes(x = pm25_referencia, y = pm25_ajustado)) +
  geom_point(alpha = 0.4, color = "black") +
  geom_smooth(method = "lm", se = FALSE, color = "blue", linetype = "solid") +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "red") +
  labs(
    title = "Fig 4. Comparación entre PM2.5 ajustado y valor de referencia",
    x = "PM2.5 de referencia (µg/m³)",
    y = "PM2.5 ajustado (µg/m³)"
  ) +
  theme_minimal() +
  theme(axis.text = element_text(size = 10))
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 314 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 314 rows containing missing values or values outside the scale range
## (`geom_point()`).

# Grafica 2

CampMar1_Senamhi  %>%
  select(fecha_hora, pm25_sensor, pm25_referencia, pm25_ajustado) %>%
  pivot_longer(cols = starts_with("pm25_"),
               names_to = "tipo",
               values_to = "pm25") %>%
  ggplot(aes(x = fecha_hora, y = pm25, color = tipo, linetype = tipo, size = tipo)) +
  geom_line() +
  scale_color_manual(
    values = c(
      "pm25_sensor" = "gray80",
      "pm25_referencia" = "steelblue",       # azul más suave
      "pm25_ajustado" = "darkgreen"
    )
  ) +
  scale_linetype_manual(
    values = c(
      "pm25_sensor" = "dotted",
      "pm25_referencia" = "solid",
      "pm25_ajustado" = "solid"
    )
  ) +
  scale_size_manual(
    values = c(
      "pm25_sensor" = 0.3,
      "pm25_referencia" = 0.8,
      "pm25_ajustado" = 0.8
    )
  ) +
  scale_x_datetime(date_breaks = "1 month", date_labels = "%b-%Y") +
  labs(
    title = "Fig 5.Serie de tiempo de PM2.5: Sensor, Referencia y Ajustado",
    x = "Fecha",
    y = "PM2.5 (µg/m³)",
    color = "Tipo de medición",
    linetype = "Tipo de medición"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`geom_line()`).

## 7. Bland-Altman:

Una vez seleccionada la mejor ecuación de calibración mediante el proceso de validación cruzada Leave-One-Week-Out, se procedió a ajustar el modelo sobre el conjunto completo de datos disponibles. El modelo seleccionado fue el modelo 11 para CampMar_1, cuya especificación fue:

pm25_referencia ~ pm25_sensor + humedad * presion

A partir de esta ecuación, se estimaron las concentraciones calibradas de PM₂.₅ para cada observación en el conjunto de datos, generando una nueva variable llamada pm25_ajustado.

Para evaluar la concordancia entre los valores calibrados y los datos del monitor de referencia, se elaboró un gráfico de Bland-Altman. En dicho gráfico, la diferencia entre el valor ajustado y el valor de referencia se graficó contra la media de ambos.

La línea azul representa una tendencia ajustada linealmente (modelo lm), la cual muestra una ligera pendiente negativa. Esto sugiere que, a medida que aumentan los niveles de PM₂.₅, el sensor calibrado tiende a subestimar ligeramente las concentraciones con respecto al monitor de referencia.

Los límites de acuerdo (líneas rojas discontinuas) fueron aproximadamente ±6.3 µg/m³, mientras que la media de las diferencias fue cercana a cero (-0.27 µg/m³), lo cual indica buena concordancia general, con una pequeña variabilidad esperada dentro de los rangos aceptables.

# Filtrar datos completos

bland_altman1 <- CampMar1_Senamhi %>%
  filter(!is.na(pm25_ajustado), !is.na(pm25_referencia),
         is.finite(pm25_ajustado), is.finite(pm25_referencia)) %>%
  mutate(
    media = (pm25_ajustado + pm25_referencia) / 2,
    diferencia = pm25_ajustado - pm25_referencia
  ) %>%
  filter(is.finite(media), is.finite(diferencia))

# Calcular límites de acuerdo
media_diff <- mean(bland_altman1$diferencia)
sd_diff <- sd(bland_altman1$diferencia)
limite_superior <- media_diff + 1.96 * sd_diff
limite_inferior <- media_diff - 1.96 * sd_diff

# Graficar Bland-Altman

ggplot(bland_altman1, aes(x = media, y = diferencia)) +
  geom_point(alpha = 0.4, color = "black") +
  geom_smooth(method = "lm", se = FALSE, color = "blue", linewidth = 1) +
  geom_hline(yintercept = media_diff, color = "black", linetype = "solid") +
  geom_hline(yintercept = limite_superior, color = "red", linetype = "dashed") +
  geom_hline(yintercept = limite_inferior, color = "red", linetype = "dashed") +
  annotate("text", x = max(bland_altman1$media, na.rm = TRUE) * 0.95, y = media_diff,
           label = paste0("Mean: ", round(media_diff, 2)), hjust = 1) +
  annotate("text", x = max(bland_altman1$media, na.rm = TRUE) * 0.95, y = limite_superior,
           label = paste0("Upper limit: ", round(limite_superior, 2)), hjust = 1) +
  annotate("text", x = max(bland_altman1$media, na.rm = TRUE) * 0.95, y = limite_inferior,
           label = paste0("Lower limit: ", round(limite_inferior, 2)), hjust = 1) +
  labs(
    title = "Fig 6. Gráfico de Bland-Altman con tendencia (PurpleAir vs. referencia)",
    x = "Media entre ajustado y referencia (µg/m³)",
    y = "Diferencia (ajustado - Referencia) (µg/m³)"
  ) +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

# Graficar Bland-Altman 2
ggplot(bland_altman1, aes(x = media, y = diferencia)) +
  geom_point(alpha = 0.4, color = "black") +
  geom_smooth(method = "loess", se = FALSE, color = "blue", linewidth = 1, span = 0.75) +
  geom_hline(yintercept = media_diff, color = "black", linetype = "solid") +
  geom_hline(yintercept = limite_superior, color = "red", linetype = "dashed") +
  geom_hline(yintercept = limite_inferior, color = "red", linetype = "dashed") +
  annotate("text", x = max(bland_altman1$media, na.rm = TRUE) * 0.95, y = media_diff,
           label = paste0("Mean: ", round(media_diff, 2)), hjust = 1) +
  annotate("text", x = max(bland_altman1$media, na.rm = TRUE) * 0.95, y = limite_superior,
           label = paste0("Upper limit: ", round(limite_superior, 2)), hjust = 1) +
  annotate("text", x = max(bland_altman1$media, na.rm = TRUE) * 0.95, y = limite_inferior,
           label = paste0("Lower limit: ", round(limite_inferior, 2)), hjust = 1) +
  labs(
    title = "Fig 7. Gráfico de Bland-Altman con tendencia (PurpleAir vs. referencia)",
    x = "Media entre ajustado y referencia (µg/m³)",
    y = "Diferencia (ajustado - Referencia) (µg/m³)"
  ) +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

Línea azul con method = “lm” (modelo lineal):

La línea azul representa un ajuste lineal simple (lm) entre la diferencia (sensor calibrado – referencia) y la media de ambas mediciones. Esta línea permite identificar si existe una tendencia sistemática a sobreestimar o subestimar en todo el rango de concentraciones. En este caso, la pendiente negativa sugiere una leve subestimación del sensor a medida que aumentan los niveles de PM₂.₅.

Línea azul con method = “loess” (suavizado local):

La línea azul corresponde a un suavizado local (loess) que captura variaciones no lineales en la diferencia entre el sensor calibrado y el monitor de referencia a lo largo del rango de concentración. Aunque refleja mejor la variabilidad local, puede ser más sensible a pocos datos extremos, como se observa en el aumento de la curva en niveles altos de PM₂.₅, donde hay menor densidad de observaciones.

Selección de modelo 2

Camp2 = modelo_7 pm25_referencia ~ pm25_sensor + temperatura * humedad

# Camp2 = modelo_7  pm25_referencia ~ pm25_sensor + temperatura * humedad

# Ajustar el modelo
modelo_final_2 <- lm(pm25_referencia ~ pm25_sensor + temperatura * humedad, data = CampMar2_Senamhi)
summary(modelo_final_2)
## 
## Call:
## lm(formula = pm25_referencia ~ pm25_sensor + temperatura * humedad, 
##     data = CampMar2_Senamhi)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -56.429  -1.651  -0.246   1.293  31.510 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         35.380394   2.217134   15.96   <2e-16 ***
## pm25_sensor          0.604833   0.002374  254.78   <2e-16 ***
## temperatura         -1.035430   0.081991  -12.63   <2e-16 ***
## humedad             -0.470913   0.033815  -13.93   <2e-16 ***
## temperatura:humedad  0.015720   0.001363   11.53   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.146 on 5388 degrees of freedom
##   (76 observations deleted due to missingness)
## Multiple R-squared:  0.935,  Adjusted R-squared:  0.935 
## F-statistic: 1.939e+04 on 4 and 5388 DF,  p-value: < 2.2e-16
# Aplicar la fórmula del modelo 8 para calcular pm25_ajustado en CampMar2_Senamhi
CampMar2_Senamhi <- CampMar2_Senamhi %>%
  mutate(pm25_ajustado = 
           33.142193 +
           0.602629 * pm25_sensor +
          (-0.892808) * temperatura +
          (-0.418474) * humedad +
           0.012680 * temperatura * humedad)

# Gráfico de comparación
ggplot(CampMar2_Senamhi, aes(x = pm25_referencia, y = pm25_ajustado)) +
  geom_point(alpha = 0.4, color = "black") +
  geom_smooth(method = "lm", se = FALSE, color = "blue", linetype = "solid") +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "red") +
  labs(
    title = "Fig 8. Comparación entre PM2.5 ajustado y valor de referencia (Camp 2)",
    x = "PM2.5 de referencia (µg/m³)",
    y = "PM2.5 ajustado (µg/m³)"
  ) +
  theme_minimal() +
  theme(axis.text = element_text(size = 10))
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 76 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 76 rows containing missing values or values outside the scale range
## (`geom_point()`).

# Gráfica de serie de tiempo para CampMar2_Senamhi
CampMar2_Senamhi %>%
  select(fecha_hora, pm25_sensor, pm25_referencia, pm25_ajustado) %>%
  pivot_longer(cols = starts_with("pm25_"),
               names_to = "tipo",
               values_to = "pm25") %>%
  ggplot(aes(x = fecha_hora, y = pm25, color = tipo, linetype = tipo, size = tipo)) +
  geom_line() +
  scale_color_manual(
    values = c(
      "pm25_sensor" = "gray80",
      "pm25_referencia" = "steelblue",
      "pm25_ajustado" = "darkgreen"
    )
  ) +
  scale_linetype_manual(
    values = c(
      "pm25_sensor" = "dotted",
      "pm25_referencia" = "solid",
      "pm25_ajustado" = "solid"
    )
  ) +
  scale_size_manual(
    values = c(
      "pm25_sensor" = 0.3,
      "pm25_referencia" = 0.8,
      "pm25_ajustado" = 0.8
    )
  ) +
  scale_x_datetime(date_breaks = "1 month", date_labels = "%b-%Y") +
  labs(
    title = "Fig 9. Serie de tiempo de PM2.5: Sensor, Referencia y Ajustado (Camp 2)",
    x = "Fecha",
    y = "PM2.5 (µg/m³)",
    color = "Tipo de medición",
    linetype = "Tipo de medición"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Filtrar datos completos

# Crear tabla para Bland-Altman
bland_altman2 <- CampMar2_Senamhi %>%
  filter(!is.na(pm25_ajustado), !is.na(pm25_referencia),
         is.finite(pm25_ajustado), is.finite(pm25_referencia)) %>%
  mutate(
    media = (pm25_ajustado + pm25_referencia) / 2,
    diferencia = pm25_ajustado - pm25_referencia
  ) %>%
  filter(is.finite(media), is.finite(diferencia))

# Calcular límites de acuerdo
media_diff <- mean(bland_altman2$diferencia)
sd_diff <- sd(bland_altman2$diferencia)
limite_superior <- media_diff + 1.96 * sd_diff
limite_inferior <- media_diff - 1.96 * sd_diff

# Verificar media real antes de graficar
cat("Media de la diferencia:", round(media_diff, 4), "\n")
## Media de la diferencia: -0.0142
# Bland-Altman con línea de tendencia lineal
ggplot(bland_altman2, aes(x = media, y = diferencia)) +
  geom_point(alpha = 0.4, color = "black") +
  geom_smooth(method = "lm", se = FALSE, color = "blue", linewidth = 1) +
  geom_hline(yintercept = media_diff, color = "black", linetype = "solid") +
  geom_hline(yintercept = limite_superior, color = "red", linetype = "dashed") +
  geom_hline(yintercept = limite_inferior, color = "red", linetype = "dashed") +
  annotate("text", x = max(bland_altman2$media, na.rm = TRUE) * 0.95, y = media_diff,
           label = paste0("Mean: ", round(media_diff, 2)), hjust = 1) +
  annotate("text", x = max(bland_altman2$media, na.rm = TRUE) * 0.95, y = limite_superior,
           label = paste0("Upper limit: ", round(limite_superior, 2)), hjust = 1) +
  annotate("text", x = max(bland_altman2$media, na.rm = TRUE) * 0.95, y = limite_inferior,
           label = paste0("Lower limit: ", round(limite_inferior, 2)), hjust = 1) +
  labs(
    title = "Fig 10. Gráfico de Bland-Altman con tendencia lineal (PurpleAir vs. referencia)",
    x = "Media entre ajustado y referencia (µg/m³)",
    y = "Diferencia (ajustado - Referencia) (µg/m³)"
  ) +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

# Bland-Altman con suavizado LOESS
ggplot(bland_altman2, aes(x = media, y = diferencia)) +
  geom_point(alpha = 0.4, color = "black") +
  geom_smooth(method = "loess", se = FALSE, color = "blue", linewidth = 1, span = 0.75) +
  geom_hline(yintercept = media_diff, color = "black", linetype = "solid") +
  geom_hline(yintercept = limite_superior, color = "red", linetype = "dashed") +
  geom_hline(yintercept = limite_inferior, color = "red", linetype = "dashed") +
  annotate("text", x = max(bland_altman2$media, na.rm = TRUE) * 0.95, y = media_diff,
           label = paste0("Mean: ", round(media_diff, 2)), hjust = 1) +
  annotate("text", x = max(bland_altman2$media, na.rm = TRUE) * 0.95, y = limite_superior,
           label = paste0("Upper limit: ", round(limite_superior, 2)), hjust = 1) +
  annotate("text", x = max(bland_altman2$media, na.rm = TRUE) * 0.95, y = limite_inferior,
           label = paste0("Lower limit: ", round(limite_inferior, 2)), hjust = 1) +
  labs(
    title = "Fig 11. Gráfico de Bland-Altman con tendencia suavizada (PurpleAir vs. referencia)",
    x = "Media entre ajustado y referencia (µg/m³)",
    y = "Diferencia (ajustado - Referencia) (µg/m³)"
  ) +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

# Selección de modelo 3

modelo_con_estacion pm25_referencia ~ pm25_sensor + temperatura * dew_point + epoca_humedad

# Camp3 = modelo_con_estacion   pm25_referencia ~ pm25_sensor + temperatura * dew_point + epoca_humedad

# Ajustar el modelo
modelo_final_3 <- lm(pm25_referencia ~ pm25_sensor + temperatura * dew_point + epoca_humedad, 
                          data = CampMar3_Senamhi)
summary(modelo_final_3)
## 
## Call:
## lm(formula = pm25_referencia ~ pm25_sensor + temperatura * dew_point + 
##     epoca_humedad, data = CampMar3_Senamhi)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -45.857  -1.631  -0.222   1.320  28.588 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           -31.503084   1.948269 -16.170   <2e-16 ***
## pm25_sensor             0.553928   0.002146 258.156   <2e-16 ***
## temperatura             2.052835   0.082262  24.955   <2e-16 ***
## dew_point               1.935712   0.121927  15.876   <2e-16 ***
## epoca_humedad          -0.405732   0.157786  -2.571   0.0102 *  
## temperatura:dew_point  -0.113852   0.004840 -23.523   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.038 on 5951 degrees of freedom
##   (184 observations deleted due to missingness)
## Multiple R-squared:  0.9388, Adjusted R-squared:  0.9387 
## F-statistic: 1.825e+04 on 5 and 5951 DF,  p-value: < 2.2e-16
# Aplicar la fórmula del modelo 8 para calcular pm25_ajustado en CampMar2_Senamhi
CampMar3_Senamhi <- CampMar3_Senamhi %>%
  mutate(pm25_ajustado = 
           -28.510520 +
           0.553317 * pm25_sensor +
           1.942324 * temperatura +
           1.714124 * dew_point +
           (-0.105699) * temperatura * dew_point +
           ifelse(epoca_humedad == "humeda", -0.412816, 0)
         )

# Gráfico de comparación
ggplot(CampMar3_Senamhi, aes(x = pm25_referencia, y = pm25_ajustado)) +
  geom_point(alpha = 0.4, color = "black") +
  geom_smooth(method = "lm", se = FALSE, color = "blue", linetype = "solid") +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "red") +
  labs(
    title = "Fig 12. Comparación entre PM2.5 ajustado y valor de referencia (Camp 3)",
    x = "PM2.5 de referencia (µg/m³)",
    y = "PM2.5 ajustado (µg/m³)"
  ) +
  theme_minimal() +
  theme(axis.text = element_text(size = 10))
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 184 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 184 rows containing missing values or values outside the scale range
## (`geom_point()`).

# Gráfica de serie de tiempo para CampMar3_Senamhi
CampMar3_Senamhi %>%
  select(fecha_hora, pm25_sensor, pm25_referencia, pm25_ajustado) %>%
  pivot_longer(cols = starts_with("pm25_"),
               names_to = "tipo",
               values_to = "pm25") %>%
  ggplot(aes(x = fecha_hora, y = pm25, color = tipo, linetype = tipo, size = tipo)) +
  geom_line() +
  scale_color_manual(
    values = c(
      "pm25_sensor" = "gray80",
      "pm25_referencia" = "steelblue",
      "pm25_ajustado" = "darkgreen"
    )
  ) +
  scale_linetype_manual(
    values = c(
      "pm25_sensor" = "dotted",
      "pm25_referencia" = "solid",
      "pm25_ajustado" = "solid"
    )
  ) +
  scale_size_manual(
    values = c(
      "pm25_sensor" = 0.3,
      "pm25_referencia" = 0.8,
      "pm25_ajustado" = 0.8
    )
  ) +
  scale_x_datetime(date_breaks = "1 month", date_labels = "%b-%Y") +
  labs(
    title = "Fig 9. Serie de tiempo de PM2.5: Sensor, Referencia y Ajustado (Camp 2)",
    x = "Fecha",
    y = "PM2.5 (µg/m³)",
    color = "Tipo de medición",
    linetype = "Tipo de medición"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).

# Filtrar datos completos

# Crear tabla para Bland-Altman
bland_altman3 <- CampMar3_Senamhi %>%
  filter(!is.na(pm25_ajustado), !is.na(pm25_referencia),
         is.finite(pm25_ajustado), is.finite(pm25_referencia)) %>%
  mutate(
    media = (pm25_ajustado + pm25_referencia) / 2,
    diferencia = pm25_ajustado - pm25_referencia
  ) %>%
  filter(is.finite(media), is.finite(diferencia))

# Calcular límites de acuerdo
media_diff <- mean(bland_altman3$diferencia)
sd_diff <- sd(bland_altman3$diferencia)
limite_superior <- media_diff + 1.96 * sd_diff
limite_inferior <- media_diff - 1.96 * sd_diff

# Verificar media real antes de graficar
cat("Media de la diferencia:", round(media_diff, 4), "\n")
## Media de la diferencia: 0.1967
# Bland-Altman con línea de tendencia lineal
ggplot(bland_altman3, aes(x = media, y = diferencia)) +
  geom_point(alpha = 0.4, color = "black") +
  geom_smooth(method = "lm", se = FALSE, color = "blue", linewidth = 1) +
  geom_hline(yintercept = media_diff, color = "black", linetype = "solid") +
  geom_hline(yintercept = limite_superior, color = "red", linetype = "dashed") +
  geom_hline(yintercept = limite_inferior, color = "red", linetype = "dashed") +
  annotate("text", x = max(bland_altman3$media, na.rm = TRUE) * 0.95, y = media_diff,
           label = paste0("Mean: ", round(media_diff, 2)), hjust = 1) +
  annotate("text", x = max(bland_altman3$media, na.rm = TRUE) * 0.95, y = limite_superior,
           label = paste0("Upper limit: ", round(limite_superior, 2)), hjust = 1) +
  annotate("text", x = max(bland_altman3$media, na.rm = TRUE) * 0.95, y = limite_inferior,
           label = paste0("Lower limit: ", round(limite_inferior, 2)), hjust = 1) +
  labs(
    title = "Fig 13. Gráfico de Bland-Altman con tendencia lineal (PurpleAir vs. referencia)",
    x = "Media entre ajustado y referencia (µg/m³)",
    y = "Diferencia (ajustado - Referencia) (µg/m³)"
  ) +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

# Bland-Altman con suavizado LOESS
ggplot(bland_altman3, aes(x = media, y = diferencia)) +
  geom_point(alpha = 0.4, color = "black") +
  geom_smooth(method = "loess", se = FALSE, color = "blue", linewidth = 1, span = 0.75) +
  geom_hline(yintercept = media_diff, color = "black", linetype = "solid") +
  geom_hline(yintercept = limite_superior, color = "red", linetype = "dashed") +
  geom_hline(yintercept = limite_inferior, color = "red", linetype = "dashed") +
  annotate("text", x = max(bland_altman3$media, na.rm = TRUE) * 0.95, y = media_diff,
           label = paste0("Mean: ", round(media_diff, 2)), hjust = 1) +
  annotate("text", x = max(bland_altman3$media, na.rm = TRUE) * 0.95, y = limite_superior,
           label = paste0("Upper limit: ", round(limite_superior, 2)), hjust = 1) +
  annotate("text", x = max(bland_altman3$media, na.rm = TRUE) * 0.95, y = limite_inferior,
           label = paste0("Lower limit: ", round(limite_inferior, 2)), hjust = 1) +
  labs(
    title = "Fig 14. Gráfico de Bland-Altman con tendencia suavizada (PurpleAir vs. referencia)",
    x = "Media entre ajustado y referencia (µg/m³)",
    y = "Diferencia (ajustado - Referencia) (µg/m³)"
  ) +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

bland_altman_limpio1 <- bland_altman1 %>%
  filter(media <= 90)  # Elimina los dos puntos más extremos

# Recalcular límites
media_diff <- mean(bland_altman_limpio1$diferencia, na.rm = TRUE)
sd_diff <- sd(bland_altman_limpio1$diferencia, na.rm = TRUE)
limite_superior <- media_diff + 1.96 * sd_diff
limite_inferior <- media_diff - 1.96 * sd_diff

# Gráfico actualizado
ggplot(bland_altman_limpio1, aes(x = media, y = diferencia)) +
  geom_point(alpha = 0.4, color = "black") +
  geom_smooth(method = "lm", se = FALSE, color = "blue", linewidth = 1) +
  geom_hline(yintercept = media_diff, color = "black") +
  geom_hline(yintercept = limite_superior, color = "red", linetype = "dashed") +
  geom_hline(yintercept = limite_inferior, color = "red", linetype = "dashed") +
  annotate("text", x = max(bland_altman_limpio1$media, na.rm = TRUE) * 0.95, y = media_diff,
           label = paste0("Mean: ", round(media_diff, 2)), hjust = 1) +
  annotate("text", x = max(bland_altman_limpio1$media, na.rm = TRUE) * 0.95, y = limite_superior,
           label = paste0("Upper limit: ", round(limite_superior, 2)), hjust = 1) +
  annotate("text", x = max(bland_altman_limpio1$media, na.rm = TRUE) * 0.95, y = limite_inferior,
           label = paste0("Lower limit: ", round(limite_inferior, 2)), hjust = 1) +
  labs(
    title = "Fig 15.Gráfico de Bland-Altman sin valores extremos",
    x = "Media entre ajustado y referencia (µg/m³)",
    y = "Diferencia (ajustado - Referencia) (µg/m³)"
  ) +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

bland_altman_limpio2 <- bland_altman2 %>%
  filter(media <= 90)  # Elimina los dos puntos más extremos

# Recalcular límites
media_diff <- mean(bland_altman_limpio2$diferencia, na.rm = TRUE)
sd_diff <- sd(bland_altman_limpio2$diferencia, na.rm = TRUE)
limite_superior <- media_diff + 1.96 * sd_diff
limite_inferior <- media_diff - 1.96 * sd_diff

# Gráfico actualizado
ggplot(bland_altman_limpio2, aes(x = media, y = diferencia)) +
  geom_point(alpha = 0.4, color = "black") +
  geom_smooth(method = "lm", se = FALSE, color = "blue", linewidth = 1) +
  geom_hline(yintercept = media_diff, color = "black") +
  geom_hline(yintercept = limite_superior, color = "red", linetype = "dashed") +
  geom_hline(yintercept = limite_inferior, color = "red", linetype = "dashed") +
  annotate("text", x = max(bland_altman_limpio2$media, na.rm = TRUE) * 0.95, y = media_diff,
           label = paste0("Mean: ", round(media_diff, 2)), hjust = 1) +
  annotate("text", x = max(bland_altman_limpio2$media, na.rm = TRUE) * 0.95, y = limite_superior,
           label = paste0("Upper limit: ", round(limite_superior, 2)), hjust = 1) +
  annotate("text", x = max(bland_altman_limpio2$media, na.rm = TRUE) * 0.95, y = limite_inferior,
           label = paste0("Lower limit: ", round(limite_inferior, 2)), hjust = 1) +
  labs(
    title = "Fig 16.Gráfico de Bland-Altman sin valores extremos",
    x = "Media entre ajustado y referencia (µg/m³)",
    y = "Diferencia (ajustado - Referencia) (µg/m³)"
  ) +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

bland_altman_limpio3 <- bland_altman3 %>%
  filter(media <= 90)  # Elimina los dos puntos más extremos

# Recalcular límites
media_diff <- mean(bland_altman_limpio3$diferencia, na.rm = TRUE)
sd_diff <- sd(bland_altman_limpio3$diferencia, na.rm = TRUE)
limite_superior <- media_diff + 1.96 * sd_diff
limite_inferior <- media_diff - 1.96 * sd_diff

# Gráfico actualizado
ggplot(bland_altman_limpio3, aes(x = media, y = diferencia)) +
  geom_point(alpha = 0.4, color = "black") +
  geom_smooth(method = "lm", se = FALSE, color = "blue", linewidth = 1) +
  geom_hline(yintercept = media_diff, color = "black") +
  geom_hline(yintercept = limite_superior, color = "red", linetype = "dashed") +
  geom_hline(yintercept = limite_inferior, color = "red", linetype = "dashed") +
  annotate("text", x = max(bland_altman_limpio3$media, na.rm = TRUE) * 0.95, y = media_diff,
           label = paste0("Mean: ", round(media_diff, 2)), hjust = 1) +
  annotate("text", x = max(bland_altman_limpio3$media, na.rm = TRUE) * 0.95, y = limite_superior,
           label = paste0("Upper limit: ", round(limite_superior, 2)), hjust = 1) +
  annotate("text", x = max(bland_altman_limpio3$media, na.rm = TRUE) * 0.95, y = limite_inferior,
           label = paste0("Lower limit: ", round(limite_inferior, 2)), hjust = 1) +
  labs(
    title = "Fig 17.Gráfico de Bland-Altman sin valores extremos",
    x = "Media entre ajustado y referencia (µg/m³)",
    y = "Diferencia (ajustado - Referencia) (µg/m³)"
  ) +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'