MODELO POLINÓMICO

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 la variable

Causa y efecto:

El índice de humedad del suelo se considera la variable independiente \((X)\), porque representa la cantidad relativa de humedad presente en el terreno y puede influir en la saturación y pérdida de estabilidad.

El área afectada del deslizamiento se considera la variable dependiente \((Y)\), porque representa el tamaño de la superficie impactada por el evento.

Dentro del modelo se establece:

\(X=\) índice de humedad del suelo.
\(Y=\) área afectada del deslizamiento.

limpiar_numerico <- function(variable) {
  variable <- as.character(variable)
  variable <- gsub(",", ".", variable)
  as.numeric(variable)
}

datos$soil_moisture_index <- limpiar_numerico(
  datos$soil_moisture_index
)

datos$area_affected_m2 <- limpiar_numerico(
  datos$area_affected_m2
)

datos_originales <- datos[
  is.finite(datos$soil_moisture_index) &
    is.finite(datos$area_affected_m2) &
    datos$soil_moisture_index >= 0 &
    datos$soil_moisture_index <= 1 &
    datos$area_affected_m2 > 0,
  c(
    "soil_moisture_index",
    "area_affected_m2"
  )
]

datos_originales <- datos_originales[
  order(
    datos_originales$soil_moisture_index
  ),
]

X_original <- datos_originales$soil_moisture_index
Y_original <- datos_originales$area_affected_m2

3.- Tabla de pares de valores original

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

cat(
  "Tamaño muestral original =",
  nrow(
    tabla_original
  )
)

## Tamaño muestral original = 11033

knitr::kable(
  head(
    tabla_original,
    20
  ),
  digits = 4,
  col.names = c(
    "Índice de humedad del suelo",
    "Área afectada del deslizamiento (m²)"
  ),
  caption = paste(
    "Tabla Nro. 1. Pares de valores originales del índice",
    "de humedad del suelo y el área afectada del deslizamiento"
  )
)

Tabla Nro. 1. Pares de valores originales del índice de humedad del suelo y el área afectada del deslizamiento
Índice de humedad del suelo	Área afectada del deslizamiento (m²)
0.0978	604.2
0.1229	920.3
0.1236	410.0
0.1240	1680.0
0.1254	982.7
0.1401	669.1
0.1494	2695.8
0.1494	1544.9
0.1499	1885.1
0.1502	1919.7
0.1510	688.2
0.1520	998.1
0.1523	1092.7
0.1524	1043.4
0.1539	573.5
0.1559	663.4
0.1576	570.9
0.1578	198.6
0.1586	3253.5
0.1634	4557.4

4.- Gráfica original de dispersión

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 el índice de humedad del suelo\n",
      "y el área afectada del deslizamiento"
    ),
    x = "Índice de humedad del suelo",
    y = "Área afectada del deslizamiento (m²)"
  ) +
  scale_x_continuous(
    limits = c(
      0,
      1
    ),
    breaks = seq(
      0,
      1,
      by = 0.1
    ),
    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,
      1
    ),
    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.- Tratamientos de datos

Simplificación de los datos:

Los registros se ordenan según el índice de humedad del suelo y se dividen en 15 grupos con cantidades similares de observaciones. En cada grupo se calcula la mediana de la humedad y la mediana del área afectada.

Este tratamiento reduce la dispersión de los datos originales y permite observar con mayor claridad la curvatura de la relación sin inventar ni modificar observaciones individuales.

numero_grupos <- 15

cortes <- quantile(
  datos_originales$soil_moisture_index,
  probs = seq(
    0,
    1,
    length.out = numero_grupos + 1
  ),
  na.rm = TRUE
)

cortes <- unique(
  cortes
)

datos_originales$grupo <- cut(
  datos_originales$soil_moisture_index,
  breaks = cortes,
  include.lowest = TRUE,
  labels = FALSE
)

datos_prom <- aggregate(
  cbind(
    soil_moisture_index,
    area_affected_m2
  ) ~ grupo,
  data = datos_originales,
  FUN = median
)

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

X <- datos_prom$soil_moisture_index
Y <- datos_prom$area_affected_m2

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

## Tamaño muestral utilizado en el modelo = 15

5.1.- Tabla pares de valores simplificada

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

knitr::kable(
  tabla_simplificada,
  digits = 4,
  col.names = c(
    "Índice de humedad representativo",
    "Área afectada representativa (m²)"
  ),
  caption = paste(
    "Tabla Nro. 2. Pares representativos del índice",
    "de humedad del suelo y el área afectada"
  )
)

Tabla Nro. 2. Pares representativos del índice de humedad del suelo y el área afectada
Índice de humedad representativo	Área afectada representativa (m²)
0.2955	1091.10
0.4054	1202.10
0.4838	1336.40
0.5533	1413.10
0.6163	1538.15
0.6748	1603.60
0.7238	1647.40
0.7727	1665.60
0.8183	1793.60
0.8610	2029.50
0.8990	2162.90
0.9321	2289.70
0.9606	2634.10
0.9834	3110.70
0.9971	4809.85

5.2.- Gráfica simplificada de dispersión

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

ggplot(
  tabla_simplificada,
  aes(
    x = X,
    y = Y
  )
) +
  geom_point(
    color = "#00A6D6",
    alpha = 1,
    size = 3.2
  ) +
  labs(
    title = "Gráfica Nro. 2",
    subtitle = paste0(
      "Diagrama simplificado de dispersión entre el índice de humedad del suelo\n",
      "y el área afectada del deslizamiento"
    ),
    x = "Índice de humedad del suelo",
    y = "Área afectada del deslizamiento (m²)"
  ) +
  scale_x_continuous(
    limits = c(
      0,
      1
    ),
    breaks = seq(
      0,
      1,
      by = 0.1
    ),
    expand = expansion(
      mult = 0,
      add = 0
    )
  ) +
  scale_y_continuous(
    limits = c(
      0,
      y_max_simplificado
    ),
    expand = expansion(
      mult = 0,
      add = 0
    )
  ) +
  coord_cartesian(
    xlim = c(
      0,
      1
    ),
    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 presenta una tendencia creciente con cambios en la pendiente. El área afectada aumenta lentamente en niveles medios de humedad y con mayor rapidez cuando el índice se aproxima a la saturación. Por esta razón se propone un modelo polinómico de grado 4.

7.- Parámetros

Modelo polinómico general de grado 4:

\[ Y = a+ bX+ cX^2+ dX^3+ eX^4 \]

Modelo polinómico aplicado al estudio:

\[ \text{Área afectada} = a+ b(\text{Humedad})+ c(\text{Humedad})^2+ d(\text{Humedad})^3+ e(\text{Humedad})^4 \]

Ajuste del modelo:

modelo_pol <- lm(
  Y ~ poly(
    X,
    4,
    raw = TRUE
  )
)

param <- coef(
  modelo_pol
)

a_est <- unname(
  param[1]
)

b_est <- unname(
  param[2]
)

c_est <- unname(
  param[3]
)

d_est <- unname(
  param[4]
)

e_est <- unname(
  param[5]
)

a_est

## [1] 15836.51

b_est

## [1] -117708.7

c_est

## [1] 330237.6

d_est

## [1] -386458.9

e_est

## [1] 162219.3

Ecuación polinómica obtenida:

\[ \widehat{Y} = 15836.5083 - 1.177e+05X + 3.302e+05X^2 - 3.865e+05X^3 + 1.622e+05X^4 \]

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

Aunque el modelo presenta una relación polinómica, se utiliza la función lm() porque la ecuación es lineal respecto a sus parámetros.

El modelo de grado 4 es:

\[ Y = a+ bX+ cX^2+ dX^3+ eX^4 \]

Se definen los términos:

\[ X_1=X, \qquad X_2=X^2, \qquad X_3=X^3, \qquad X_4=X^4 \]

y los parámetros:

\[ \beta_0=a, \qquad \beta_1=b, \qquad \beta_2=c, \qquad \beta_3=d, \qquad \beta_4=e \]

Por tanto:

\[ Y = \beta_0+ \beta_1X_1+ \beta_2X_2+ \beta_3X_3+ \beta_4X_4 \]

Como la ecuación es lineal respecto a los parámetros, puede estimarse en R mediante:

modelo_pol <- lm(
  Y ~ poly(
    X,
    4,
    raw = TRUE
  )
)

Finalmente, se reconstruye el modelo:

\[ \text{Área afectada} = a+ b(\text{Humedad})+ c(\text{Humedad})^2+ d(\text{Humedad})^3+ e(\text{Humedad})^4 \]

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

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

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

X_curva <- seq(
  x_min_calibracion,
  x_max_calibracion,
  length.out = 500
)

Y_curva <- a_est +
  b_est * X_curva +
  c_est * X_curva^2 +
  d_est * X_curva^3 +
  e_est * X_curva^4

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

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

ggplot() +
  geom_point(
    data = tabla_simplificada,
    aes(
      x = X,
      y = Y
    ),
    color = "#00A6D6",
    alpha = 1,
    size = 3.2
  ) +
  geom_line(
    data = curva_polinomica,
    aes(
      x = X,
      y = Y
    ),
    color = "red",
    linewidth = 1.5
  ) +
  labs(
    title = "Gráfica Nro. 3",
    subtitle = paste0(
      "Modelo polinómico de grado 4 entre el índice de humedad del suelo\n",
      "y el área afectada del deslizamiento"
    ),
    x = "Índice de humedad del suelo",
    y = "Área afectada del deslizamiento (m²)"
  ) +
  scale_x_continuous(
    limits = c(
      0,
      1
    ),
    breaks = seq(
      0,
      1,
      by = 0.1
    ),
    expand = expansion(
      mult = 0,
      add = 0
    )
  ) +
  scale_y_continuous(
    limits = c(
      0,
      y_max_modelo
    ),
    expand = expansion(
      mult = 0,
      add = 0
    )
  ) +
  coord_cartesian(
    xlim = c(
      0,
      1
    ),
    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(
  X,
  Y
) * 100

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

cat(
  "Coeficiente de correlación (r) =",
  round(
    r,
    2
  ),
  "%\n"
)

## Coeficiente de correlación (r) = 76.1 %

cat(
  "Coeficiente de determinación (R²) =",
  round(
    R2,
    2
  ),
  "%"
)

## Coeficiente de determinación (R²) = 89.75 %

El coeficiente de correlación es de 76.1 %, lo que representa una relación positiva fuerte entre el índice de humedad y el área afectada.

El coeficiente de determinación es de 89.75 %, lo que indica el porcentaje de variabilidad de los puntos representativos explicado por el modelo polinómico de grado 4.

10.- Restricciones

humedad_dominio <- seq(
  0,
  1,
  length.out = 10001
)

areas_dominio <- a_est +
  b_est * humedad_dominio +
  c_est * humedad_dominio^2 +
  d_est * humedad_dominio^3 +
  e_est * humedad_dominio^4

area_min_dominio <- min(
  areas_dominio,
  na.rm = TRUE
)

area_max_dominio <- max(
  areas_dominio,
  na.rm = TRUE
)

modelo_fuera_dominio <- any(
  areas_dominio < 0
)

respuesta_restriccion <- ifelse(
  modelo_fuera_dominio,
  "Sí",
  "No"
)

Dominios:

\[ D_X = \left\{ X\in\mathbb{R}: 0\leq X\leq1 \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 no genera áreas negativas dentro del dominio físico de la humedad. Sus estimaciones se encuentran aproximadamente entre 954.15 y 15836.51 m².

El modelo fue calibrado con índices de humedad representativos entre 0.2955 y 0.9971. Los valores externos a este intervalo, aunque permanezcan entre 0 y 1, corresponden a extrapolaciones y no deben interpretarse con la misma confianza.

11.- Estimación

# Pregunta de cantidad
humedad_objetivo <- 0.80

area_objetivo <- a_est +
  b_est * humedad_objetivo +
  c_est * humedad_objetivo^2 +
  d_est * humedad_objetivo^3 +
  e_est * humedad_objetivo^4

# Pregunta de porcentaje
incremento_pct <- 0.10

humedad_final <- humedad_objetivo *
  (
    1 +
      incremento_pct
  )

if (
  humedad_final > 1
) {
  stop(
    "El nuevo índice de humedad supera el dominio físico de 1."
  )
}

area_final <- a_est +
  b_est * humedad_final +
  c_est * humedad_final^2 +
  d_est * humedad_final^3 +
  e_est * humedad_final^4

porcentaje_cambio_area <- (
  (
    area_final -
      area_objetivo
  ) /
    area_objetivo
) * 100

Pregunta de cantidad

¿Cuál es el área afectada esperada cuando el índice de humedad del suelo es 0.80?

Resultado estimado: 1599.62 m²

Pregunta de porcentaje

¿En qué porcentaje cambia el área afectada esperada cuando el índice de humedad del suelo aumenta un 10 %, pasando de 0.80 a 0.88?

Resultado estimado: 19.4 %

12.- Conclusión

En conclusión:

Entre el índice de humedad del suelo y el área afectada del deslizamiento existe una relación no lineal representada mediante el modelo polinómico de grado 4: