MODELO POTENCIAL

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 las variables

La precipitación acumulada en 7 días se selecciona como la variable independiente o causa \((X)\), porque el aumento de agua acumulada puede contribuir a la saturación del terreno.

El índice de susceptibilidad a deslizamientos se selecciona como la variable dependiente o efecto \((Y)\), porque representa la respuesta que se desea estimar a partir de la precipitación.

Dentro del modelo se establece:

  • \(X=\) precipitación acumulada en 7 días.
  • \(Y=\) índice de susceptibilidad a deslizamientos.
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$susceptibility_index <- limpiar_numerico(
  datos$susceptibility_index
)

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

Y_original <- datos_modelo$susceptibility_index
X_original <- datos_modelo$precipitation_7d_mm

3.- Tabla pares de valores

n_original <- nrow(
  datos_modelo
)

cat(
  "Tamaño muestral =",
  n_original
)
## Tamaño muestral = 11033
tabla_original <- data.frame(
  X = X_original,
  Y = Y_original
)

knitr::kable(
  head(
    tabla_original,
    20
  ),
  digits = 2,
  col.names = c(
    "Precipitación acumulada en 7 días (mm)",
    "Índice de susceptibilidad"
  ),
  caption = paste(
    "Tabla Nro. 1. Pares de valores de la precipitación acumulada",
    "en 7 días y el índice de susceptibilidad a deslizamientos"
  )
)
Tabla Nro. 1. Pares de valores de la precipitación acumulada en 7 días y el índice de susceptibilidad a deslizamientos
Precipitación acumulada en 7 días (mm) Índice de susceptibilidad
60.21 83.27
115.77 93.48
47.14 87.15
52.52 74.52
12.51 37.58
84.97 89.74
121.28 93.52
83.39 80.28
93.87 81.31
105.27 89.49
44.63 61.80
60.71 80.24
9.67 54.34
67.48 58.96
55.84 86.14
31.55 49.00
33.14 49.75
47.54 69.49
23.37 47.01
122.71 92.30

4.- Gráfica de dispersión

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

ggplot(
  tabla_original,
  aes(
    x = X,
    y = Y
  )
) +
  geom_point(
    color = "#00A6D6",
    alpha = 0.45,
    size = 1.5
  ) +
  labs(
    title = "Gráfica Nro. 1",
    subtitle = paste0(
      "Diagrama de dispersión entre la precipitación acumulada en 7 días\n",
      "y el índice de susceptibilidad a deslizamientos"
    ),
    x = "Precipitación acumulada en 7 días (mm)",
    y = "Índice de susceptibilidad"
  ) +
  scale_x_continuous(
    limits = c(
      0,
      x_max_original
    ),
    expand = expansion(
      mult = 0,
      add = 0
    )
  ) +
  scale_y_continuous(
    limits = c(
      0,
      100
    ),
    breaks = seq(
      0,
      100,
      by = 10
    ),
    expand = expansion(
      mult = 0,
      add = 0
    )
  ) +
  coord_cartesian(
    xlim = c(
      0,
      x_max_original
    ),
    ylim = c(
      0,
      100
    ),
    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

Con el fin de reducir la dispersión de los datos, se calcula el promedio de la precipitación acumulada en 7 días para cada valor del índice de susceptibilidad. Este procedimiento disminuye la variabilidad de las observaciones y facilita la identificación de la tendencia entre ambas variables.

datos_prom <- aggregate(
  precipitation_7d_mm ~ susceptibility_index,
  data = datos_modelo,
  FUN = mean
)

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

Y <- datos_prom$susceptibility_index
X <- datos_prom$precipitation_7d_mm

cat(
  "Tamaño muestral utilizado en el modelo =",
  nrow(
    datos_prom
  )
)
## Tamaño muestral utilizado en el modelo = 4874

5.1.- Tabla pares de valores simplificada

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

knitr::kable(
  head(
    tabla_simplificada,
    20
  ),
  digits = 2,
  col.names = c(
    "Precipitación promedio en 7 días (mm)",
    "Índice de susceptibilidad"
  ),
  caption = paste(
    "Tabla Nro. 2. Pares de valores promediados de la precipitación",
    "acumulada en 7 días y el índice de susceptibilidad"
  )
)
Tabla Nro. 2. Pares de valores promediados de la precipitación acumulada en 7 días y el índice de susceptibilidad
Precipitación promedio en 7 días (mm) Índice de susceptibilidad
0.50 13.62
10.30 18.03
5.79 18.34
7.27 20.97
0.50 21.64
17.46 22.06
14.83 22.46
12.55 23.06
1.47 23.07
17.89 23.24
9.67 23.36
8.41 24.10
2.49 24.22
0.50 24.90
9.47 25.15
7.51 25.19
12.80 25.29
23.68 25.40
13.30 25.42
21.79 25.46

5.2.- Gráfica de dispersión simplificada

x_max_simplificado <- ceiling(
  max(
    X,
    na.rm = TRUE
  ) / 50
) * 50

ggplot(
  tabla_simplificada,
  aes(
    x = X,
    y = Y
  )
) +
  geom_point(
    color = "#00A6D6",
    alpha = 0.65,
    size = 1.8
  ) +
  labs(
    title = "Gráfica Nro. 2",
    subtitle = paste0(
      "Diagrama de dispersión entre la precipitación acumulada en 7 días\n",
      "y el índice de susceptibilidad a deslizamientos"
    ),
    x = "Precipitación acumulada en 7 días (mm)",
    y = "Índice de susceptibilidad"
  ) +
  scale_x_continuous(
    limits = c(
      0,
      x_max_simplificado
    ),
    expand = expansion(
      mult = 0,
      add = 0
    )
  ) +
  scale_y_continuous(
    limits = c(
      0,
      100
    ),
    breaks = seq(
      0,
      100,
      by = 10
    ),
    expand = expansion(
      mult = 0,
      add = 0
    )
  ) +
  coord_cartesian(
    xlim = c(
      0,
      x_max_simplificado
    ),
    ylim = c(
      0,
      100
    ),
    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 en el gráfico presenta una tendencia creciente, lo que indica que un modelo potencial podría describir adecuadamente la relación entre la precipitación acumulada en 7 días y el índice de susceptibilidad. Conforme aumenta la precipitación, la susceptibilidad también se incrementa, reflejando una relación positiva y no lineal entre ambas variables.

Modelo potencial general

\[ Y=aX^b \]

Modelo potencial aplicado al estudio

\[ \text{Susceptibilidad} = a \left( \text{Precipitación}_{7d} \right)^b \]

7.- Parámetros

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

ln_a <- coef(
  modelo_pot
)[1]

a_est <- exp(
  ln_a
)

b_est <- coef(
  modelo_pot
)[2]

ln_a
## (Intercept) 
##    3.036235
a_est
## (Intercept) 
##    20.82668
b_est
##    log(X) 
## 0.3031985

Modelo potencial obtenido:

\[ \widehat{Y} = 20.8267 X^{0.303198} \]

El parámetro \(a\) es la constante del modelo y el parámetro \(b\) es el exponente que determina la forma del crecimiento potencial.

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

Modelo potencial general

\[ Y=aX^b \]

Aplicando logaritmos:

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

Definimos:

\[ Y_1=\ln(Y) \]

\[ X_1=\ln(X) \]

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

Obtenemos un modelo lineal:

\[ Y_1= \beta_0+ \beta_1X_1 \]

Modelo potencial aplicado al estudio

\[ \text{Susceptibilidad} = a \left( \text{Precipitación}_{7d} \right)^b \]

Definiendo:

\[ Y_1= \ln( \text{Susceptibilidad} ) \]

\[ X_1= \ln( \text{Precipitación}_{7d} ) \]

Obtenemos:

\[ \ln( \text{Susceptibilidad} ) = \beta_0+ \beta_1 \ln( \text{Precipitación}_{7d} ) \]

Esta ecuación corresponde a una regresión lineal respecto a los parámetros, por lo que puede ser estimada mediante:

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

Finalmente:

\[ a=e^{\beta_0} \]

\[ b=\beta_1 \]

x_max_modelo <- x_max_simplificado

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

Y_curva <- ifelse(
  X_curva == 0,
  0,
  a_est * X_curva^b_est
)

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

ggplot() +
  geom_point(
    data = tabla_simplificada,
    aes(
      x = X,
      y = Y
    ),
    color = "#00A6D6",
    alpha = 0.65,
    size = 1.8
  ) +
  geom_line(
    data = curva_potencial,
    aes(
      x = X,
      y = Y
    ),
    color = "red",
    linewidth = 1.5
  ) +
  labs(
    title = "Gráfica Nro. 3",
    subtitle = paste0(
      "Modelo potencial entre la precipitación acumulada en 7 días\n",
      "y el índice de susceptibilidad a deslizamientos"
    ),
    x = "Precipitación acumulada en 7 días (mm)",
    y = "Índice de susceptibilidad"
  ) +
  scale_x_continuous(
    limits = c(
      0,
      x_max_modelo
    ),
    expand = expansion(
      mult = 0,
      add = 0
    )
  ) +
  scale_y_continuous(
    limits = c(
      0,
      100
    ),
    breaks = seq(
      0,
      100,
      by = 10
    ),
    expand = expansion(
      mult = 0,
      add = 0
    )
  ) +
  coord_cartesian(
    xlim = c(
      0,
      x_max_modelo
    ),
    ylim = c(
      0,
      100
    ),
    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),
  log(Y)
) * 100

cat(
  "Coeficiente de correlación (r) =",
  round(
    r,
    2
  )
)
## Coeficiente de correlación (r) = 83.5

El coeficiente de correlación es de 83.5 %, lo que indica una relación positiva fuerte entre los logaritmos de la precipitación acumulada y la susceptibilidad.

10.- Restricciones

x_limite_susceptibilidad <- (
  100 / a_est
)^(1 / b_est)

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

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

Dominio de la variable independiente utilizado por el modelo:

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

Dominio físico del índice de susceptibilidad:

\[ D_Y= \left\{ Y\in\mathbb{R}: 0\leq Y\leq100 \right\} \]

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

Respuesta: Sí.

La primera restricción es que la precipitación debe ser mayor que cero, porque el ajuste se realiza mediante \(\ln(X)\), que no está definido para \(X=0\).

Además, al resolver:

\[ aX^b\leq100 \]

se obtiene:

\[ X \leq \left( \frac{100}{a} \right)^{\frac{1}{b}} \]

Con los parámetros obtenidos:

\[ X \leq 176.73 \text{ mm} \]

Por tanto, valores mayores que ese límite generan una susceptibilidad estimada superior a 100 y quedan fuera del dominio físico de \(Y\). El modelo fue calibrado con valores de \(X\) entre 0.5 y 304.16 mm.

11.- Estimación

precipitacion_cantidad <- 80

susceptibilidad_80 <- round(
  a_est *
    precipitacion_cantidad^b_est
)

precipitacion_inicial <- 50
precipitacion_final <- 80

susceptibilidad_50 <- a_est *
  precipitacion_inicial^b_est

susceptibilidad_80_real <- a_est *
  precipitacion_final^b_est

incremento_porcentaje <- (
  (
    susceptibilidad_80_real -
      susceptibilidad_50
  ) /
    susceptibilidad_50
) * 100

Pregunta de cantidad

¿Cuál es el índice de susceptibilidad esperado cuando la precipitación acumulada en 7 días es de 80 mm?

Resultado: 79

Pregunta de porcentaje

¿En qué porcentaje aumenta la susceptibilidad esperada al incrementar la precipitación acumulada en 7 días de 50 a 80 mm?

Resultado: 15.32 %

12.- Conclusión

En conclusión:

Entre la precipitación acumulada en 7 días y el índice de susceptibilidad a deslizamientos existe una relación de tipo potencial, cuyo modelo es \(\widehat{Y}=20.8267X^{0.303198}\). El coeficiente de correlación fue de 83.5 %, lo que evidencia una relación positiva fuerte. La principal restricción es que la precipitación debe ser mayor que cero debido a la transformación logarítmica y, para conservar la susceptibilidad dentro del intervalo de 0 a 100, el modelo debe aplicarse hasta aproximadamente 176.73 mm. Los resultados indican que, a medida que aumenta la precipitación acumulada, la susceptibilidad también se incrementa.