REGRESIÓN MÚLTIPLE LINEAL

0.- Carga de librerías

library(readxl)
library(knitr)
library(scatterplot3d)

1.- Carga de datos

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

2.- Selección de tres variables (Causa y Efecto)

Causa y efecto:

La precipitación acumulada en 7 días y la pendiente del terreno se consideran variables independientes porque ambas pueden influir en la inestabilidad del terreno.

El índice de susceptibilidad a deslizamientos se considera la variable dependiente porque representa la respuesta que se desea estimar.

Dentro del modelo se establece:

  • \(Y=\) índice de susceptibilidad a deslizamientos.
  • \(X_1=\) precipitación acumulada en 7 días.
  • \(X_2=\) pendiente del terreno.
limpiar_numerico <- function(variable) {
  variable <- as.character(variable)
  variable <- gsub(",", ".", variable)
  as.numeric(variable)
}

y_inicial <- limpiar_numerico(
  datos$susceptibility_index
)

x1_inicial <- limpiar_numerico(
  datos$precipitation_7d_mm
)

x2_inicial <- limpiar_numerico(
  datos$terrain_slope_deg
)

n_inicial <- nrow(
  datos
)

cat(
  "Tamaño muestral inicial =",
  n_inicial
)
## Tamaño muestral inicial = 11033

3.- Tabla de Tripleta de Valores

Se eliminan únicamente los registros incompletos o que quedan fuera del dominio físico de las variables.

TPV <- data.frame(
  y = y_inicial,
  x1 = x1_inicial,
  x2 = x2_inicial
)

TPV <- TPV[
  is.finite(TPV$y) &
    is.finite(TPV$x1) &
    is.finite(TPV$x2) &
    TPV$y >= 0 &
    TPV$y <= 100 &
    TPV$x1 >= 0 &
    TPV$x2 >= 0 &
    TPV$x2 <= 90,
]

row.names(
  TPV
) <- NULL

n_modelo <- nrow(
  TPV
)

cat(
  "Tamaño muestral del modelo =",
  n_modelo
)
## Tamaño muestral del modelo = 11033

El tamaño muestral disminuye únicamente cuando existen datos faltantes o valores que no pertenecen a los dominios físicos establecidos.

tabla_tripletas <- data.frame(
  Nro = seq_len(
    n_modelo
  ),
  Y = TPV$y,
  X1 = TPV$x1,
  X2 = TPV$x2
)

knitr::kable(
  head(
    tabla_tripletas,
    20
  ),
  digits = 2,
  col.names = c(
    "N°",
    "Susceptibilidad Y",
    "Precipitación 7 días X1 (mm)",
    "Pendiente X2 (grados)"
  ),
  caption = paste(
    "Tabla Nro. 1. Tripletas de valores del índice de susceptibilidad,",
    "la precipitación acumulada y la pendiente del terreno"
  )
)
Tabla Nro. 1. Tripletas de valores del índice de susceptibilidad, la precipitación acumulada y la pendiente del terreno
N° Susceptibilidad Y Precipitación 7 días X1 (mm) Pendiente X2 (grados)
1 83.27 60.21 28.14
2 93.48 115.77 33.00
3 87.15 47.14 19.77
4 74.52 52.52 22.85
5 37.58 12.51 15.65
6 89.74 84.97 23.76
7 93.52 121.28 14.44
8 80.28 83.39 11.02
9 81.31 93.87 16.17
10 89.49 105.27 25.24
11 61.80 44.63 14.28
12 80.24 60.71 28.05
13 54.34 9.67 25.72
14 58.96 67.48 8.34
15 86.14 55.84 31.26
16 49.00 31.55 8.05
17 49.75 33.14 8.53
18 69.49 47.54 15.82
19 47.01 23.37 18.10
20 92.30 122.71 22.80

4.- Gráfica de Dispersión

x1 <- TPV$x1
x2 <- TPV$x2
y <- TPV$y

x1_max_grafica <- ceiling(
  max(
    x1,
    na.rm = TRUE
  ) / 50
) * 50

x2_max_grafica <- ceiling(
  max(
    x2,
    na.rm = TRUE
  ) / 10
) * 10

color_puntos <- rgb(
  0,
  166,
  214,
  105,
  maxColorValue = 255
)

par(
  mar = c(
    5,
    6,
    5,
    3
  ),
  col.main = "#6D213C"
)

scatterplot3d(
  x = x1,
  y = x2,
  z = y,
  angle = 225,
  pch = 16,
  color = color_puntos,
  cex.symbols = 0.58,
  main = paste0(
    "Gráfica Nro. 1: Diagrama de dispersión entre la susceptibilidad,\n",
    "la precipitación acumulada en 7 días y la pendiente del terreno"
  ),
  xlab = "Precipitación acumulada 7 días X1 (mm)",
  ylab = "Pendiente del terreno X2 (grados)",
  zlab = "Índice de susceptibilidad Y",
  xlim = c(
    0,
    x1_max_grafica
  ),
  ylim = c(
    0,
    x2_max_grafica
  ),
  zlim = c(
    0,
    100
  ),
  grid = TRUE,
  box = TRUE,
  y.margin.add = 0.8,
  las = 1
)

5.- Conjetura

Se plantea que la precipitación acumulada en 7 días \((X_1)\) y la pendiente del terreno \((X_2)\) están asociadas con un aumento del índice de susceptibilidad \((Y)\). Por ello, la nube tridimensional se ajustará mediante un plano lineal.

6.- Parámetros

Modelo matemático:

\[ Y=a+bX_1+cX_2 \]

Ajuste del modelo:

regresion_multiple <- lm(
  y ~ x1 + x2,
  data = TPV
)

a <- unname(
  coef(
    regresion_multiple
  )[1]
)

b <- unname(
  coef(
    regresion_multiple
  )[2]
)

c <- unname(
  coef(
    regresion_multiple
  )[3]
)

a
## [1] 42.87213
b
## [1] 0.2759242
c
## [1] 0.5934665

Ecuación múltiple obtenida:

\[ \widehat{Y} = 42.87213 + 0.27592X_1 + 0.59347X_2 \]

El intercepto \(a\) representa la susceptibilidad teórica cuando la precipitación y la pendiente son iguales a cero. El coeficiente \(b\) representa el cambio promedio de la susceptibilidad por cada milímetro adicional de precipitación, manteniendo constante la pendiente. El coeficiente \(c\) representa el cambio promedio de la susceptibilidad por cada grado adicional de pendiente, manteniendo constante la precipitación.

7.- Gráfica de comparación de la realidad con el modelo

Vista en perspectiva

par(
  mar = c(
    5,
    6,
    5,
    3
  ),
  col.main = "#6D213C"
)

grafica_perspectiva <- scatterplot3d(
  x = x1,
  y = x2,
  z = y,
  angle = 225,
  pch = 16,
  color = color_puntos,
  cex.symbols = 0.58,
  main = paste0(
    "Gráfica Nro. 2: Modelo múltiple lineal — vista en perspectiva\n",
    "Susceptibilidad, precipitación acumulada y pendiente"
  ),
  xlab = "Precipitación acumulada 7 días X1 (mm)",
  ylab = "Pendiente del terreno X2 (grados)",
  zlab = "Índice de susceptibilidad Y",
  xlim = c(
    0,
    x1_max_grafica
  ),
  ylim = c(
    0,
    x2_max_grafica
  ),
  zlim = c(
    0,
    z_max_modelo
  ),
  grid = TRUE,
  box = TRUE,
  y.margin.add = 0.8,
  las = 1
)

grafica_perspectiva$plane3d(
  regresion_multiple,
  col = "#C0392B",
  lty = 2,
  lwd = 1.5,
  draw_polygon = FALSE
)

La vista en perspectiva permite observar simultáneamente la nube de puntos y la inclinación general del plano de regresión.

Vista frontal

par(
  mar = c(
    5,
    6,
    5,
    3
  ),
  col.main = "#6D213C"
)

grafica_frontal <- scatterplot3d(
  x = x1,
  y = x2,
  z = y,
  angle = 45,
  pch = 16,
  color = color_puntos,
  cex.symbols = 0.58,
  main = paste0(
    "Gráfica Nro. 3: Modelo múltiple lineal — vista frontal\n",
    "Susceptibilidad, precipitación acumulada y pendiente"
  ),
  xlab = "Precipitación acumulada 7 días X1 (mm)",
  ylab = "Pendiente del terreno X2 (grados)",
  zlab = "Índice de susceptibilidad Y",
  xlim = c(
    0,
    x1_max_grafica
  ),
  ylim = c(
    0,
    x2_max_grafica
  ),
  zlim = c(
    0,
    z_max_modelo
  ),
  grid = TRUE,
  box = TRUE,
  y.margin.add = 0.8,
  las = 1
)

grafica_frontal$plane3d(
  regresion_multiple,
  col = "#C0392B",
  lty = 2,
  lwd = 1.5,
  draw_polygon = FALSE
)

La vista frontal facilita la observación de la relación entre la precipitación, la susceptibilidad y la posición del plano.

Vista lateral

par(
  mar = c(
    5,
    6,
    5,
    3
  ),
  col.main = "#6D213C"
)

grafica_lateral <- scatterplot3d(
  x = x1,
  y = x2,
  z = y,
  angle = 135,
  pch = 16,
  color = color_puntos,
  cex.symbols = 0.58,
  main = paste0(
    "Gráfica Nro. 4: Modelo múltiple lineal — vista lateral\n",
    "Susceptibilidad, precipitación acumulada y pendiente"
  ),
  xlab = "Precipitación acumulada 7 días X1 (mm)",
  ylab = "Pendiente del terreno X2 (grados)",
  zlab = "Índice de susceptibilidad Y",
  xlim = c(
    0,
    x1_max_grafica
  ),
  ylim = c(
    0,
    x2_max_grafica
  ),
  zlim = c(
    0,
    z_max_modelo
  ),
  grid = TRUE,
  box = TRUE,
  y.margin.add = 0.8,
  las = 1
)

grafica_lateral$plane3d(
  regresion_multiple,
  col = "#C0392B",
  lty = 2,
  lwd = 1.5,
  draw_polygon = FALSE
)

La vista lateral permite apreciar con mayor claridad el efecto de la pendiente del terreno sobre el plano de regresión.

8.- Test de Bondad

8.1.- Coeficiente de correlación múltiple

R_multiple_porcentaje <- r_multiple * 100

cat(
  "Coeficiente de correlación múltiple (R) =",
  round(
    R_multiple_porcentaje,
    2
  ),
  "%"
)
## Coeficiente de correlación múltiple (R) = 82.6 %

El coeficiente de correlación múltiple es de 82.6 %, lo que indica una asociación fuerte entre la susceptibilidad y el conjunto formado por la precipitación y la pendiente.

8.2.- Coeficiente de determinación

R2_porcentaje <- r2 * 100

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

El coeficiente de determinación es de 68.22 %. Este valor representa el porcentaje de variabilidad del índice de susceptibilidad explicado conjuntamente por la precipitación acumulada y la pendiente del terreno.

9.- Restricciones

min_x1 <- min(
  x1,
  na.rm = TRUE
)

max_x1 <- max(
  x1,
  na.rm = TRUE
)

min_x2 <- min(
  x2,
  na.rm = TRUE
)

max_x2 <- max(
  x2,
  na.rm = TRUE
)

predicciones_calibracion <- predict(
  regresion_multiple
)

prediccion_min_calibracion <- min(
  predicciones_calibracion,
  na.rm = TRUE
)

prediccion_max_calibracion <- max(
  predicciones_calibracion,
  na.rm = TRUE
)

resultados_fuera_dominio <- sum(
  predicciones_calibracion < 0 |
    predicciones_calibracion > 100
)

limite_x1_pendiente_0 <- (
  100 -
    a
) / b

limite_x1_pendiente_90 <- (
  100 -
    a -
    c * 90
) / b

Dominios:

\[ D_{X_1} = \left\{ X_1\in\mathbb{R}: X_1\geq0 \right\} \]

\[ D_{X_2} = \left\{ X_2\in\mathbb{R}: 0\leq X_2\leq90 \right\} \]

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

Pregunta: ¿Existe alguna combinación de \(X_1\) y \(X_2\) que, reemplazada en el plano, genere un valor fuera del dominio de \(Y\)?

Respuesta: Sí.

El modelo solo produce resultados físicamente válidos cuando:

\[ 0 \leq a+bX_1+cX_2 \leq 100 \]

Como los coeficientes de precipitación y pendiente son positivos, algunas combinaciones elevadas generan susceptibilidades mayores que 100. El límite superior de precipitación depende de la pendiente:

\[ X_1 \leq \frac{ 100-a-cX_2 }{ b } \]

Cuando la pendiente es \(0^\circ\), el límite de precipitación es aproximadamente 207.04 mm. Cuando la pendiente es \(90^\circ\), el límite es aproximadamente 13.47 mm.

Dentro de los datos utilizados, el plano genera estimaciones entre 43.95 y 187.48. Se identificaron 612 combinaciones cuya estimación supera el dominio de la susceptibilidad.

El modelo fue calibrado con precipitaciones entre 0.5 y 507.07 mm, y pendientes entre 1.5° y 54.99°. Las combinaciones externas a estos intervalos corresponden a extrapolaciones.

10.- Estimación

precipitacion_objetivo <- 100
pendiente_objetivo <- 30

susceptibilidad_estimada <- a +
  b * precipitacion_objetivo +
  c * pendiente_objetivo

¿Qué índice de susceptibilidad se espera cuando la precipitación acumulada en 7 días es de 100 mm y la pendiente del terreno es de 30°?

Resultado: 88.27

11.- Conclusión

En conclusión:

Entre el índice de susceptibilidad a deslizamientos, la precipitación acumulada en 7 días y la pendiente del terreno existe una relación lineal multivariable definida mediante el plano:

\[ \widehat{Y} = 42.87213 + 0.27592X_1 + 0.59347X_2 \]

El coeficiente de correlación múltiple fue de 82.6 % y el coeficiente de determinación fue de 68.22 %. Esto indica el porcentaje de la variabilidad de la susceptibilidad explicado conjuntamente por la precipitación y la pendiente.

El modelo presenta restricciones porque algunas combinaciones elevadas de precipitación y pendiente generan susceptibilidades superiores a 100. Por ello debe respetarse la condición \(0\leq a+bX_1+cX_2\leq100\) y evitarse la extrapolación fuera de los intervalos de calibración.