MODELO LOGARÍTMICO

0.- Carga de librerías

library(readxl)
library(ggplot2)
library(knitr)

1.- Carga de datos

datos <- read_excel(
  "dataset_landslides.xlsx",
  sheet = "Dataset_landslides"
)

2.- Definición de variables

Causa y efecto:

El incremento de la precipitación acumulada en 7 días puede aumentar el área de terreno afectada por un deslizamiento.

Dentro del modelo se establece:

  • \(X=\) precipitación acumulada en 7 días.
  • \(Y=\) área afectada del deslizamiento.
limpiar_numerico <- function(variable) {
  variable <- as.character(variable)
  variable <- gsub(",", ".", variable)
  as.numeric(variable)
}

datos$precipitation_7d_mm <- limpiar_numerico(
  datos$precipitation_7d_mm
)

datos$area_affected_m2 <- limpiar_numerico(
  datos$area_affected_m2
)

datos_modelo <- datos[
  is.finite(datos$precipitation_7d_mm) &
    is.finite(datos$area_affected_m2) &
    datos$precipitation_7d_mm > 0 &
    datos$area_affected_m2 > 0,
  c(
    "precipitation_7d_mm",
    "area_affected_m2"
  )
]

X_original <- datos_modelo$precipitation_7d_mm
Y_original <- datos_modelo$area_affected_m2

3.- Tabla de pares de valores original

tabla_original <- data.frame(
  X = X_original,
  Y = Y_original
)

n_original <- nrow(
  tabla_original
)

cat(
  "Tamaño muestral original =",
  n_original
)
## Tamaño muestral original = 11033
knitr::kable(
  head(
    tabla_original,
    20
  ),
  digits = 2,
  col.names = c(
    "Precipitación acumulada en 7 días (mm)",
    "Área afectada (m²)"
  ),
  caption = paste(
    "Tabla Nro. 1. Pares de valores originales de la precipitación",
    "acumulada en 7 días y el área afectada del deslizamiento"
  )
)
Tabla Nro. 1. Pares de valores originales de la precipitación acumulada en 7 días y el área afectada del deslizamiento
Precipitación acumulada en 7 días (mm) Área afectada (m²)
60.21 1879.8
115.77 2183.0
47.14 784.5
52.52 1815.4
12.51 1399.4
84.97 2568.7
121.28 2725.3
83.39 2167.9
93.87 1983.9
105.27 1835.7
44.63 413.2
60.71 1609.5
9.67 1917.7
67.48 1624.4
55.84 2450.9
31.55 971.8
33.14 885.1
47.54 881.4
23.37 2208.4
122.71 4071.5

4.- Gráfica original de dispersión

x_max_original <- ceiling(
  max(
    X_original,
    na.rm = TRUE
  ) / 50
) * 50

y_max_original <- ceiling(
  max(
    Y_original,
    na.rm = TRUE
  ) / 5000
) * 5000

ggplot(
  tabla_original,
  aes(
    x = X,
    y = Y
  )
) +
  geom_point(
    color = "#00A6D6",
    alpha = 0.40,
    size = 1.5
  ) +
  labs(
    title = "Gráfica Nro. 1",
    subtitle = paste0(
      "Diagrama original de dispersión entre la precipitación acumulada en 7 días\n",
      "y el área afectada del deslizamiento"
    ),
    x = "Precipitación acumulada en 7 días (mm)",
    y = "Área afectada del deslizamiento (m²)"
  ) +
  scale_x_continuous(
    limits = c(
      0,
      x_max_original
    ),
    expand = expansion(
      mult = 0,
      add = 0
    )
  ) +
  scale_y_continuous(
    limits = c(
      0,
      y_max_original
    ),
    expand = expansion(
      mult = 0,
      add = 0
    )
  ) +
  coord_cartesian(
    xlim = c(
      0,
      x_max_original
    ),
    ylim = c(
      0,
      y_max_original
    ),
    expand = FALSE
  ) +
  theme_bw(
    base_size = 14
  ) +
  theme(
    plot.title = element_text(
      hjust = 0.5,
      face = "bold",
      color = "#6D213C",
      size = 17
    ),
    plot.subtitle = element_text(
      hjust = 0.5,
      face = "bold",
      size = 13,
      lineheight = 1.1
    ),
    axis.title = element_text(
      face = "bold"
    ),
    panel.grid.minor = element_blank()
  )

5.- Tratamiento de datos

Para identificar con mayor claridad la tendencia logarítmica, se utiliza un dominio de calibración comprendido entre 50 y 220 mm de precipitación. Luego se aplica el rango intercuartílico IQR al área afectada para reducir la influencia de valores extremos. Finalmente, la precipitación se agrupa en intervalos de 5 mm y se calcula el área afectada promedio de cada grupo.

x_min_modelo <- 50
x_max_modelo <- 220

datos_dominio <- datos_modelo[
  datos_modelo$precipitation_7d_mm >= x_min_modelo &
    datos_modelo$precipitation_7d_mm <= x_max_modelo,
]

Q1_Y <- quantile(
  datos_dominio$area_affected_m2,
  0.25,
  na.rm = TRUE
)

Q3_Y <- quantile(
  datos_dominio$area_affected_m2,
  0.75,
  na.rm = TRUE
)

IQR_Y <- Q3_Y - Q1_Y

lim_inf_Y <- Q1_Y - 1.5 * IQR_Y
lim_sup_Y <- Q3_Y + 1.5 * IQR_Y

datos_filtrados <- datos_dominio[
  datos_dominio$area_affected_m2 >= lim_inf_Y &
    datos_dominio$area_affected_m2 <= lim_sup_Y,
]

datos_filtrados$precipitacion_grupo <- round(
  datos_filtrados$precipitation_7d_mm / 5
) * 5

datos_prom <- aggregate(
  area_affected_m2 ~ precipitacion_grupo,
  data = datos_filtrados,
  FUN = mean
)

datos_prom <- datos_prom[
  order(
    datos_prom$precipitacion_grupo
  ),
]

X <- datos_prom$precipitacion_grupo
Y <- datos_prom$area_affected_m2

cat(
  "Datos originales =",
  nrow(datos_modelo),
  "\n"
)
## Datos originales = 11033
cat(
  "Datos dentro del dominio de 50 a 220 mm =",
  nrow(datos_dominio),
  "\n"
)
## Datos dentro del dominio de 50 a 220 mm = 6805
cat(
  "Datos después del filtrado IQR =",
  nrow(datos_filtrados),
  "\n"
)
## Datos después del filtrado IQR = 6419
cat(
  "Pares promediados utilizados en el modelo =",
  nrow(datos_prom)
)
## Pares promediados utilizados en el modelo = 35

5.1.- Tabla de pares de valores simplificada

tabla_simplificada <- data.frame(
  X = X,
  Y = Y
)

n_modelo <- nrow(
  tabla_simplificada
)

cat(
  "Tamaño muestral del modelo =",
  n_modelo
)
## Tamaño muestral del modelo = 35
knitr::kable(
  tabla_simplificada,
  digits = 2,
  col.names = c(
    "Precipitación agrupada en 7 días (mm)",
    "Área afectada promedio (m²)"
  ),
  caption = paste(
    "Tabla Nro. 2. Pares de valores tratados de la precipitación",
    "acumulada en 7 días y el área afectada del deslizamiento"
  )
)
Tabla Nro. 2. Pares de valores tratados de la precipitación acumulada en 7 días y el área afectada del deslizamiento
Precipitación agrupada en 7 días (mm) Área afectada promedio (m²)
50 1718.65
55 1791.23
60 1930.19
65 1983.03
70 2027.57
75 2073.87
80 2162.30
85 2246.73
90 2637.12
95 2472.60
100 2690.44
105 2663.40
110 2711.86
115 2789.50
120 2965.42
125 3012.48
130 2960.87
135 3283.50
140 3437.81
145 3414.68
150 3523.85
155 3522.32
160 3720.22
165 3717.09
170 4011.88
175 3625.15
180 3875.25
185 4063.09
190 4200.54
195 4583.79
200 4466.79
205 4161.78
210 4862.48
215 3738.19
220 4368.53

El tamaño muestral disminuye debido a la delimitación del dominio, la eliminación de valores extremos y la agrupación de la precipitación.

5.2.- Gráfica simplificada de dispersión

y_max_simplificado <- ceiling(
  max(
    Y,
    na.rm = TRUE
  ) / 1000
) * 1000

ggplot(
  tabla_simplificada,
  aes(
    x = X,
    y = Y
  )
) +
  geom_point(
    color = "#00A6D6",
    alpha = 0.80,
    size = 2.5
  ) +
  labs(
    title = "Gráfica Nro. 2",
    subtitle = paste0(
      "Diagrama simplificado entre la precipitación acumulada en 7 días\n",
      "y el área afectada promedio del deslizamiento"
    ),
    x = "Precipitación acumulada en 7 días (mm)",
    y = "Área afectada promedio del deslizamiento (m²)"
  ) +
  scale_x_continuous(
    limits = c(
      x_min_modelo,
      x_max_modelo
    ),
    breaks = seq(
      x_min_modelo,
      x_max_modelo,
      by = 20
    ),
    expand = expansion(
      mult = 0,
      add = 0
    )
  ) +
  scale_y_continuous(
    limits = c(
      0,
      y_max_simplificado
    ),
    breaks = pretty(
      c(
        0,
        y_max_simplificado
      ),
      n = 8
    ),
    expand = expansion(
      mult = 0,
      add = 0
    )
  ) +
  coord_cartesian(
    xlim = c(
      x_min_modelo,
      x_max_modelo
    ),
    ylim = c(
      0,
      y_max_simplificado
    ),
    expand = FALSE
  ) +
  theme_bw(
    base_size = 14
  ) +
  theme(
    plot.title = element_text(
      hjust = 0.5,
      face = "bold",
      color = "#6D213C",
      size = 17
    ),
    plot.subtitle = element_text(
      hjust = 0.5,
      face = "bold",
      size = 13,
      lineheight = 1.1
    ),
    axis.title = element_text(
      face = "bold"
    ),
    panel.grid.minor = element_blank()
  )

6.- Conjetura

La distribución de los puntos tratados presenta una tendencia ascendente que aumenta con mayor rapidez al inicio y reduce progresivamente su pendiente. Por ello, se propone un modelo logarítmico para representar la relación entre la precipitación acumulada en 7 días y el área afectada del deslizamiento.

Modelo logarítmico general:

\[ Y=a+b\ln(X) \]

Modelo logarítmico aplicado al estudio:

\[ \text{Área afectada} = a+ b\ln( \text{Precipitación}_{7d} ) \]

7.- Parámetros

modelo_log <- lm(
  Y ~ log(X)
)

a_est <- unname(
  coef(modelo_log)[1]
)

b_est <- unname(
  coef(modelo_log)[2]
)

a_est
## [1] -6464.977
b_est
## [1] 1999.981

Modelo logarítmico obtenido:

\[ \widehat{Y} = -6464.9767 + 1999.9813 \ln(X) \]

El intercepto \(a\) representa el valor teórico del área afectada y el coeficiente \(b\) determina el cambio asociado al logaritmo natural de la precipitación.

Justificación del uso de la regresión lineal (lm)

Aunque el modelo presenta una relación logarítmica, se utiliza la función lm() porque puede expresarse como un modelo lineal respecto a sus parámetros.

Modelo logarítmico general:

\[ Y=a+b\ln(X) \]

Aplicando la transformación:

\[ X_1=\ln(X), \qquad \beta_0=a, \qquad \beta_1=b \]

Se obtiene el modelo lineal:

\[ Y= \beta_0+ \beta_1X_1 \]

Esta ecuación puede estimarse mediante:

modelo_log <- lm(
  Y ~ log(X)
)

Finalmente, los coeficientes permiten reconstruir el modelo logarítmico:

\[ Y=a+b\ln(X) \]

8.- Comparación de la realidad con el modelo

X_curva <- seq(
  x_min_modelo,
  x_max_modelo,
  length.out = 500
)

Y_curva <- a_est +
  b_est * log(
    X_curva
  )

curva_logaritmica <- data.frame(
  X = X_curva,
  Y = Y_curva
)

y_max_modelo <- ceiling(
  max(
    c(
      Y,
      Y_curva
    ),
    na.rm = TRUE
  ) / 1000
) * 1000

ggplot() +
  geom_point(
    data = tabla_simplificada,
    aes(
      x = X,
      y = Y
    ),
    color = "#00A6D6",
    alpha = 0.80,
    size = 2.5
  ) +
  geom_line(
    data = curva_logaritmica,
    aes(
      x = X,
      y = Y
    ),
    color = "red",
    linewidth = 1.5
  ) +
  labs(
    title = "Gráfica Nro. 3",
    subtitle = paste0(
      "Modelo logarítmico entre la precipitación acumulada en 7 días\n",
      "y el área afectada del deslizamiento"
    ),
    x = "Precipitación acumulada en 7 días (mm)",
    y = "Área afectada del deslizamiento (m²)"
  ) +
  scale_x_continuous(
    limits = c(
      x_min_modelo,
      x_max_modelo
    ),
    breaks = seq(
      x_min_modelo,
      x_max_modelo,
      by = 20
    ),
    expand = expansion(
      mult = 0,
      add = 0
    )
  ) +
  scale_y_continuous(
    limits = c(
      0,
      y_max_modelo
    ),
    breaks = pretty(
      c(
        0,
        y_max_modelo
      ),
      n = 8
    ),
    expand = expansion(
      mult = 0,
      add = 0
    )
  ) +
  coord_cartesian(
    xlim = c(
      x_min_modelo,
      x_max_modelo
    ),
    ylim = c(
      0,
      y_max_modelo
    ),
    expand = FALSE
  ) +
  theme_bw(
    base_size = 14
  ) +
  theme(
    plot.title = element_text(
      hjust = 0.5,
      face = "bold",
      color = "#6D213C",
      size = 17
    ),
    plot.subtitle = element_text(
      hjust = 0.5,
      face = "bold",
      size = 13,
      lineheight = 1.1
    ),
    axis.title = element_text(
      face = "bold"
    ),
    panel.grid.minor = element_blank(),
    legend.position = "none"
  )

9.- Test de bondad

r <- cor(
  log(X),
  Y,
  method = "pearson"
) * 100

R2 <- summary(
  modelo_log
)$r.squared * 100

cat(
  "Coeficiente de correlación de Pearson =",
  round(
    r,
    2
  ),
  "%\n"
)
## Coeficiente de correlación de Pearson = 96.47 %
cat(
  "Coeficiente de determinación =",
  round(
    R2,
    2
  ),
  "%"
)
## Coeficiente de determinación = 93.06 %

El coeficiente de correlación entre el logaritmo de la precipitación y el área afectada es de 96.47 %, lo que evidencia una relación positiva fuerte.

El coeficiente de determinación es de 93.06 %.

10.- Restricciones

x_limite_fisico <- exp(
  -a_est / b_est
)

x_min_calibracion <- min(
  X,
  na.rm = TRUE
)

x_max_calibracion <- max(
  X,
  na.rm = TRUE
)

modelo_fuera_dominio <- any(
  Y_curva < 0
)

Dominios matemáticos:

\[ D_X= \left\{ X\in\mathbb{R}: X>0 \right\} \]

\[ D_Y= \left\{ Y\in\mathbb{R}: Y\geq0 \right\} \]

Pregunta: ¿Existe algún valor de \(X\) que, reemplazado en el modelo matemático, genere un valor fuera del dominio de \(Y\)?

Respuesta: No.

El modelo presenta una restricción matemática porque el logaritmo natural solo está definido para valores de precipitación mayores que cero.

El límite físico se obtiene resolviendo:

\[ a+b\ln(X)\geq0 \]

Por tanto:

\[ X \geq e^{-\frac{a}{b}} \]

Con los parámetros obtenidos:

\[ X \geq 25.34 \text{ mm} \]

El modelo debe utilizarse preferentemente dentro del intervalo de calibración:

\[ 50 \leq X \leq 220 \]

11.- Estimación

precipitacion_cantidad <- median(
  X,
  na.rm = TRUE
)

area_cantidad <- a_est +
  b_est * log(
    precipitacion_cantidad
  )

precipitacion_inicial <- 100
precipitacion_final <- 150

area_inicial <- a_est +
  b_est * log(
    precipitacion_inicial
  )

area_final <- a_est +
  b_est * log(
    precipitacion_final
  )

incremento_area_porcentaje <- (
  (
    area_final -
      area_inicial
  ) /
    area_inicial
) * 100

Pregunta de cantidad

¿Cuál es el área afectada esperada cuando la precipitación acumulada en 7 días es de 135 mm?

Resultado: 3345.48 m²

12.- Conclusión

En conclusión:

Entre la precipitación acumulada en 7 días y el área afectada del deslizamiento existe una relación de tipo logarítmico, representada por el modelo Y = -6464.9767 + 1999.9813ln(X) . La precipitación acumulada debe ser mayor que cero para la aplicación del modelo. El modelo indica que el área afectada del deslizamiento aumenta a medida que se incrementa la precipitación acumulada en 7 días.