PROYECTO DE REGRESIONES

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:

La latitud se considera la variable independiente \((X)\), porque representa la posición geográfica en la que ocurre el evento.

El tamaño del deslizamiento, representado por el área afectada, se considera la variable dependiente \((Y)\), porque es la respuesta que se desea analizar respecto a la ubicación latitudinal.

Dentro del modelo se establece:

• \(X=\) latitud.
• \(Y=\) área afectada del deslizamiento.

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

corregir_latitud <- function(variable) {
  variable <- limpiar_numerico(variable)

  indice <- is.finite(variable) &
    abs(variable) > 90

  while (
    any(indice)
  ) {
    variable[indice] <- variable[indice] / 10

    indice <- is.finite(variable) &
      abs(variable) > 90
  }

  variable
}

datos$latitude <- corregir_latitud(
  datos$latitude
)

datos$area_affected_m2 <- limpiar_numerico(
  datos$area_affected_m2
)

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

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

X_original <- datos_originales$latitude
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(
    "Latitud (grados)",
    "Área afectada del deslizamiento (m²)"
  ),
  caption = paste(
    "Tabla Nro. 1. Pares de valores originales de la latitud",
    "y el área afectada del deslizamiento"
  )
)

Tabla Nro. 1. Pares de valores originales de la latitud y el área afectada del deslizamiento

Latitud (grados)	Área afectada del deslizamiento (m²)
-89.8035	1621.5
-89.0797	7988.2
-88.9610	2143.6
-88.7500	676.8
-88.6150	1150.2
-88.5590	4474.5
-88.1890	6114.6
-88.1760	2154.1
-88.1370	5922.1
-88.1060	1318.5
-87.7346	1482.2
-87.5880	1466.2
-87.4010	2238.0
-86.9930	1091.3
-86.9740	718.3
-86.2648	1758.6
-86.2520	1580.7
-85.9550	2550.5
-85.5790	5069.9
-85.3480	2606.7

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 la latitud\n",
      "y el área afectada del deslizamiento"
    ),
    x = "Latitud (grados)",
    y = "Área afectada del deslizamiento (m²)"
  ) +
  scale_x_continuous(
    limits = c(
      -90,
      90
    ),
    breaks = seq(
      -90,
      90,
      by = 30
    ),
    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(
      -90,
      90
    ),
    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:

Inicialmente se eliminaron los valores faltantes y se ordenaron los registros por latitud. Después, los datos fueron divididos en 15 grupos con cantidades similares de observaciones. Para cada grupo se calculó la mediana de la latitud y la mediana del área afectada.

Este tratamiento reduce la dispersión y permite observar con mayor claridad los cambios de dirección de la relación polinómica sin modificar la tendencia general de los datos.

numero_grupos <- 15

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

cortes <- unique(
  cortes
)

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

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

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

X <- datos_prom$latitude
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(
    "Latitud representativa (grados)",
    "Área afectada representativa (m²)"
  ),
  caption = paste(
    "Tabla Nro. 2. Pares representativos de la latitud",
    "y el área afectada del deslizamiento"
  )
)

Tabla Nro. 2. Pares representativos de la latitud y el área afectada del deslizamiento

Latitud representativa (grados)	Área afectada representativa (m²)
-41.5621	1876.80
-0.9035	1814.60
12.1974	1842.55
17.4612	1861.10
24.5938	1896.80
27.5687	1896.90
30.2413	1860.75
33.3746	1883.40
36.0624	1862.65
38.6230	1843.20
40.9892	1785.30
44.1705	1698.10
45.8396	1733.30
48.7490	1811.50
64.0583	1780.60

5.2.- Gráfica simplificada de dispersión

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

x_min_simplificado <- floor(
  min(
    X,
    na.rm = TRUE
  ) / 10
) * 10

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

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 la latitud\n",
      "y el área afectada del deslizamiento"
    ),
    x = "Latitud (grados)",
    y = "Área afectada del deslizamiento (m²)"
  ) +
  scale_x_continuous(
    limits = c(
      x_min_simplificado,
      x_max_simplificado
    ),
    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(
      x_min_simplificado,
      x_max_simplificado
    ),
    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 cambios de dirección y una tendencia no lineal. Por esta razón, se propone un modelo polinómico de grado 4 para representar la relación entre la latitud y el área afectada del deslizamiento.

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{Latitud})+ c(\text{Latitud})^2+ d(\text{Latitud})^3+ e(\text{Latitud})^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] 1800.63

b_est

## [1] 8.383504

c_est

## [1] -0.1161255

d_est

## [1] -0.005412836

e_est

## [1] 0.00007954896

Ecuación polinómica obtenida:

\[ \widehat{Y} = 1800.6302 + 8.384e+00X - 1.161e-01X^2 - 5.413e-03X^3 + 7.955e-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 \]

Para aplicar regresión lineal, se define cada término como una variable explicativa:

\[ 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, el modelo queda expresado como:

\[ 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 polinómico aplicado al estudio:

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

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

X_curva <- seq(
  min(
    X,
    na.rm = TRUE
  ),
  max(
    X,
    na.rm = TRUE
  ),
  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 la latitud\n",
      "y el área afectada del deslizamiento"
    ),
    x = "Latitud (grados)",
    y = "Área afectada del deslizamiento (m²)"
  ) +
  scale_x_continuous(
    limits = c(
      x_min_simplificado,
      x_max_simplificado
    ),
    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(
      x_min_simplificado,
      x_max_simplificado
    ),
    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) = -45.58 %

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

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

El coeficiente de correlación es de -45.58 %, lo que representa una relación negativa moderada entre la latitud y el área afectada.

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

10.- Restricciones

latitudes_dominio <- seq(
  -90,
  90,
  length.out = 10001
)

areas_dominio <- a_est +
  b_est * latitudes_dominio +
  c_est * latitudes_dominio^2 +
  d_est * latitudes_dominio^3 +
  e_est * latitudes_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"
)

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

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

Dominios:

\[ D_X = \left\{ X\in\mathbb{R}: -90\leq X\leq90 \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 polinómico no genera áreas negativas dentro del dominio físico de la latitud. Las estimaciones se encuentran aproximadamente entre 1633.74 y 9270.66 m².

El modelo fue calibrado con valores representativos de latitud entre -41.5621° y 64.0583°. Los valores fuera de este intervalo, aunque pertenezcan al dominio físico, corresponden a extrapolaciones.

11.- Estimación

# Pregunta de cantidad
latitud_objetivo <- median(
  X,
  na.rm = TRUE
)

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

# Pregunta de porcentaje
incremento_pct <- 0.10

nueva_latitud <- latitud_objetivo *
  (
    1 +
      incremento_pct
  )

if (
  nueva_latitud > 90 ||
    nueva_latitud < -90
) {
  stop(
    "La nueva latitud queda fuera del dominio físico."
  )
}

area_actual <- area_objetivo

area_nueva <- a_est +
  b_est * nueva_latitud +
  c_est * nueva_latitud^2 +
  d_est * nueva_latitud^3 +
  e_est * nueva_latitud^4

porcentaje_cambio_area <- (
  (
    area_nueva -
      area_actual
  ) /
    area_actual
) * 100

Pregunta de cantidad

¿Cuál es el área afectada esperada cuando la latitud es 33.3746°?

Resultado estimado: 1848.55 m²

Pregunta de porcentaje

¿En qué porcentaje cambia el área afectada esperada al incrementar la latitud seleccionada un 10 %?

Resultado estimado: -1.08 %

12.- Conclusión

En conclusión:

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

\[ \widehat{Y} = 1800.6302 + 8.384e+00X - 1.161e-01X^2 - 5.413e-03X^3 + 7.955e-05X^4 \]

El coeficiente de correlación fue de -45.58 % y el coeficiente de determinación fue de 64.17 %. El tratamiento mediante 15 grupos representativos permitió observar con mayor claridad los cambios de dirección de la relación. El modelo no genera áreas negativas dentro del dominio físico de la latitud.