MODELO LINEAL

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

El índice de humedad del suelo se considera la variable independiente o causa (X), porque representa la humedad presente en el terreno y puede influir en su estabilidad.

El índice de susceptibilidad a deslizamientos se considera la variable dependiente o efecto (Y), porque representa la respuesta que se desea estimar a partir de la humedad del suelo.

Dentro del modelo se establece:

  • X = índice de humedad del suelo.
  • Y = índice de susceptibilidad a deslizamientos.

Las dos variables son cuantitativas continuas y se expresan mediante índices. 

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$susceptibility_index <- limpiar_numerico(
  datos$susceptibility_index
)

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

X_original <- datos_originales$soil_moisture_index
Y_original <- datos_originales$susceptibility_index

3.- Tabla pares de valores

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

cat(
  "Tamaño muestral de la tabla original =",
  nrow(tabla_original),
  "pares de datos.\n"
)
## Tamaño muestral de la tabla original = 11033 pares de datos.
knitr::kable(
  head(tabla_original, 20),
  digits = 4,
  col.names = c(
    "Índice de humedad del suelo (X)",
    "Índice de susceptibilidad (Y)"
  ),
  caption = "Tabla Nro. 1. Pares de valores originales"
)
Tabla Nro. 1. Pares de valores originales
Índice de humedad del suelo (X) Índice de susceptibilidad (Y)
0.6991 83.27
0.9316 93.48
0.8163 87.15
0.8226 74.52
0.5915 37.58
0.9407 89.74
0.9310 93.52
0.9337 80.28
0.8726 81.31
0.9837 89.49
0.7566 61.80
0.6429 80.24
0.3742 54.34
0.6196 58.96
0.8753 86.14
0.5863 49.00
0.5639 49.75
0.5146 69.49
0.4690 47.01
0.9472 92.30

4.- Gráfica de Dispersión

ggplot(
  tabla_original,
  aes(
    x = X,
    y = Y
  )
) +
  geom_point(
    color = "#00A6D6",
    alpha = 0.50,
    size = 1.5
  ) +
  geom_vline(
    xintercept = 0,
    color = "black",
    linewidth = 0.7
  ) +
  geom_hline(
    yintercept = 0,
    color = "black",
    linewidth = 0.7
  ) +
  labs(
    title = "Gráfica Nro. 1",
    subtitle = paste(
     "Diagrama de dispersión entre el índice de humedad del suelo\n",
    "y la susceptibilidad a deslizamientos"
    ),
    x = "Índice de humedad del suelo",
    y = "Índice de susceptibilidad a deslizamientos"
  ) +
  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, 100),
    breaks = seq(0, 100, by = 10),
    expand = expansion(mult = 0, add = 0)
  ) +
  coord_cartesian(
    xlim = c(0, 1),
    ylim = c(0, 100),
    expand = FALSE,
    clip = "on"
  ) +
  theme_bw(base_size = 14) +
  theme(
    plot.title = element_text(
      hjust = 0.5,
      face = "bold",
      color = "#6D213C"
    ),
    plot.subtitle = element_text(
      hjust = 0.5,
      face = "bold"
    ),
    axis.title = element_text(face = "bold"),
    panel.grid.minor = element_blank()
  )

5.- Tratamiento de datos

Para disminuir la influencia de valores extremos se aplica el rango intercuartílico IQR a las dos variables. Después, los valores de susceptibilidad se promedian para cada valor repetido del índice de humedad.

Q1_X <- quantile(X_original, 0.25, na.rm = TRUE)
Q3_X <- quantile(X_original, 0.75, na.rm = TRUE)
IQR_X <- Q3_X - Q1_X

lim_inf_X <- Q1_X - 1.5 * IQR_X
lim_sup_X <- Q3_X + 1.5 * IQR_X

Q1_Y <- quantile(Y_original, 0.25, na.rm = TRUE)
Q3_Y <- quantile(Y_original, 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_originales[
  datos_originales$soil_moisture_index >= lim_inf_X &
    datos_originales$soil_moisture_index <= lim_sup_X &
    datos_originales$susceptibility_index >= lim_inf_Y &
    datos_originales$susceptibility_index <= lim_sup_Y,
]

datos_prom <- aggregate(
  susceptibility_index ~ soil_moisture_index,
  data = datos_filtrados,
  FUN = mean
)

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

X <- datos_prom$soil_moisture_index
Y <- datos_prom$susceptibility_index

cat("Datos originales =", nrow(datos_originales), "\n")
## Datos originales = 11033
cat("Datos después del filtrado IQR =", nrow(datos_filtrados), "\n")
## Datos después del filtrado IQR = 10995
cat("Pares promediados utilizados =", nrow(datos_prom), "\n")
## Pares promediados utilizados = 5343

5.1.- Tabla pares de valores simplificada

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

cat(
  "Tamaño muestral del modelo =",
  nrow(tabla_simplificada),
  "pares de datos promediados.\n"
)
## Tamaño muestral del modelo = 5343 pares de datos promediados.
knitr::kable(
  head(tabla_simplificada, 20),
  digits = 4,
  col.names = c(
    "Índice de humedad del suelo (X)",
    "Índice de susceptibilidad (Y)"
  ),
  caption = "Tabla Nro. 2. Pares de valores tratados y promediados"
)
Tabla Nro. 2. Pares de valores tratados y promediados
Índice de humedad del suelo (X) Índice de susceptibilidad (Y)
0.1229 37.130
0.1236 41.620
0.1240 35.070
0.1254 33.700
0.1494 52.475
0.1499 38.140
0.1502 49.300
0.1510 43.240
0.1520 50.710
0.1524 34.480
0.1559 34.120
0.1576 33.240
0.1586 54.490
0.1634 51.250
0.1640 37.570
0.1661 44.220
0.1673 34.830
0.1692 32.700
0.1710 57.950
0.1717 37.820

5.2.- Gráfica de dispersión simplificada

ggplot(
  tabla_simplificada,
  aes(
    x = X,
    y = Y
  )
) +
  geom_point(
    color = "#00A6D6",
    alpha = 0.60,
    size = 1.6
  ) +
  geom_vline(
    xintercept = 0,
    color = "black",
    linewidth = 0.7
  ) +
  geom_hline(
    yintercept = 0,
    color = "black",
    linewidth = 0.7
  ) +
  labs(
    title = "Gráfica Nro. 2",
    subtitle = paste(
      "Diagrama de dispersión simplificado entre el índice de\n",
      "humedad del suelo y la susceptibilidad a deslizamientos"
    ),
    x = "Índice de humedad del suelo",
    y = "Índice de susceptibilidad a deslizamientos"
  ) +
  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, 100),
    breaks = seq(0, 100, by = 10),
    expand = expansion(mult = 0, add = 0)
  ) +
  coord_cartesian(
    xlim = c(0, 1),
    ylim = c(0, 100),
    expand = FALSE,
    clip = "on"
  ) +
  theme_bw(base_size = 14) +
  theme(
    plot.title = element_text(
      hjust = 0.5,
      face = "bold",
      color = "#6D213C"
    ),
    plot.subtitle = element_text(
      hjust = 0.5,
      face = "bold"
    ),
    axis.title = element_text(face = "bold"),
    panel.grid.minor = element_blank()
  )

6.- Conjetura

La distribución de los puntos presenta una tendencia ascendente aproximadamente recta. Por ello, se considera que el índice de susceptibilidad aumenta de forma lineal cuando aumenta el índice de humedad del suelo.

El modelo propuesto es:

\[ Y = a + bX \]

Aplicado al estudio:

\[ \text{Susceptibilidad} = a + b(\text{Índice de humedad del suelo}) \]

7.- Parámetros

modelo_lineal <- lm(
  Y ~ X
)

a <- unname(coef(modelo_lineal)[1])
b <- unname(coef(modelo_lineal)[2])

Modelo lineal obtenido:

\[ \widehat{Y} = 28.6775 + 63.7284X \]

El intercepto \(a\) representa el índice teórico de susceptibilidad cuando el índice de humedad es igual a cero. La pendiente \(b\) indica cuánto cambia, en promedio, la susceptibilidad cuando la humedad aumenta una unidad.

Justificación matemática del uso de la regresión lineal

El modelo propuesto tiene la forma:

\[ Y = a + bX \]

La función lm(Y ~ X) estima los parámetros mediante mínimos cuadrados ordinarios, minimizando:

\[ \sum_{i=1}^{n} \left( Y_i-\widehat{Y}_i \right)^2 \]

La pendiente y el intercepto se expresan como:

\[ b = \frac{ \sum_{i=1}^{n} (X_i-\overline{X}) (Y_i-\overline{Y}) }{ \sum_{i=1}^{n} (X_i-\overline{X})^2 } \]

\[ a = \overline{Y} - b\overline{X} \]

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

La representación comienza en el origen de los ejes. La recta se dibuja desde \(X=0\) hasta \(X=1\); cuando \(X=0\), su valor es el intercepto \(a\), por lo que no necesariamente atraviesa el punto \((0,0)\).

Y_estimado <- predict(modelo_lineal)

tabla_comparacion <- data.frame(
  X = X,
  Y_real = Y,
  Y_modelo = Y_estimado
)

knitr::kable(
  head(tabla_comparacion, 20),
  digits = 4,
  col.names = c(
    "Humedad del suelo",
    "Susceptibilidad real",
    "Susceptibilidad del modelo"
  ),
  caption = "Comparación de la realidad con el modelo"
)
Comparación de la realidad con el modelo
Humedad del suelo Susceptibilidad real Susceptibilidad del modelo
0.1229 37.130 36.5097
0.1236 41.620 36.5543
0.1240 35.070 36.5798
0.1254 33.700 36.6690
0.1494 52.475 38.1985
0.1499 38.140 38.2304
0.1502 49.300 38.2495
0.1510 43.240 38.3005
0.1520 50.710 38.3642
0.1524 34.480 38.3897
0.1559 34.120 38.6127
0.1576 33.240 38.7211
0.1586 54.490 38.7848
0.1634 51.250 39.0907
0.1640 37.570 39.1289
0.1661 44.220 39.2628
0.1673 34.830 39.3392
0.1692 32.700 39.4603
0.1710 57.950 39.5750
0.1717 37.820 39.6196 
ggplot(
  tabla_simplificada,
  aes(
    x = X,
    y = Y
  )
) +
  geom_point(
    color = "#00A6D6",
    alpha = 0.60,
    size = 1.6
  ) +
  geom_vline(
    xintercept = 0,
    color = "black",
    linewidth = 0.7
  ) +
  geom_hline(
    yintercept = 0,
    color = "black",
    linewidth = 0.7
  ) +
  geom_segment(
    aes(
      x = 0,
      y = a,
      xend = 1,
      yend = a + b
    ),
    color = "red",
    linewidth = 1.5,
    inherit.aes = FALSE
  ) +
  labs(
    title = "Gráfica Nro. 3",
    subtitle = paste(
      "Comparación entre el índice de humedad del suelo\n",
      "y la susceptibilidad a deslizamientos"
    ),
    x = "Índice de humedad del suelo",
    y = "Índice de susceptibilidad a deslizamientos"
  ) +
  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, 100),
    breaks = seq(0, 100, by = 10),
    expand = expansion(mult = 0, add = 0)
  ) +
  coord_cartesian(
    xlim = c(0, 1),
    ylim = c(0, 100),
    expand = FALSE,
    clip = "on"
  ) +
  theme_bw(base_size = 14) +
  theme(
    plot.title = element_text(
      hjust = 0.5,
      face = "bold",
      color = "#6D213C"
    ),
    plot.subtitle = element_text(
      hjust = 0.5,
      face = "bold"
    ),
    axis.title = element_text(face = "bold"),
    panel.grid.minor = element_blank()
  )

9.- Test de Bondad

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

R2 <- summary(modelo_lineal)$r.squared

cat(
  "Coeficiente de correlación de Pearson =",
  round(r * 100, 2),
  "%\n"
)
## Coeficiente de correlación de Pearson = 88.04 %
cat(
  "Coeficiente de determinación =",
  round(R2 * 100, 2),
  "%\n"
)
## Coeficiente de determinación = 77.5 %

El coeficiente de Pearson es de 88.04 %, por lo que existe una relación lineal positiva fuerte.

El coeficiente de determinación es de 77.5 %.

10.- Restricciones

Y_X0 <- a
Y_X1 <- a + b

Y_min <- min(Y_X0, Y_X1)
Y_max <- max(Y_X0, Y_X1)

restriccion <- ifelse(
  Y_min >= 0 &&
    Y_max <= 100,
  "No",
  "Sí"
)

Dominios físicos de las variables:

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

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

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

Respuesta: No.

Para \(X=0\):

\[ \widehat{Y} = 28.6775 \]

Para \(X=1\):

\[ \widehat{Y} = 92.4059 \]

El modelo genera valores entre 28.6775 y 92.4059, por lo que todos permanecen dentro del dominio físico de Y..

11.- Estimación

humedad_cantidad <- 0.80

susceptibilidad_cantidad <- a +
  b * humedad_cantidad

humedad_inicial <- 0.50
humedad_final <- 0.80

susceptibilidad_inicial <- a +
  b * humedad_inicial

susceptibilidad_final <- a +
  b * humedad_final

porcentaje_aumento <- (
  (
    susceptibilidad_final -
      susceptibilidad_inicial
  ) /
    susceptibilidad_inicial
) * 100

Pregunta de cantidad

¿Cuál es el índice de susceptibilidad esperado cuando el índice de humedad del suelo es 0.80?

Resultado: 79.66

Pregunta de porcentaje

¿En qué porcentaje aumenta la susceptibilidad esperada cuando el índice de humedad del suelo aumenta de 0.50 a 0.80?

Resultado: 31.58 %

12.- Conclusión

En conclusión:

Entre el índice de humedad del suelo y el índice de susceptibilidad a deslizamientos existe una relación lineal positiva. El modelo obtenido es Y = 28.6775 + 63.7284X .

El coeficiente de Pearson fue de 88.04 % y el coeficiente de determinación fue de 77.5 %. Los datos fueron filtrados mediante IQR y promediados para cada valor repetido de humedad.

La ecuación genera valores entre 28.6775 y 92.4059 dentro del dominio de X; por tanto, no produce valores fuera del dominio físico de la susceptibilidad.