MODELO LOGARÍTMICO

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 variable

Causa y efecto:

El incremento de la precipitación acumulada en 7 días puede aumentar la cantidad de terreno afectado por un deslizamiento.

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

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

# Variable independiente (X): precipitación acumulada en 7 días
X_original <- datos_modelo$precipitation_7d_mm

# Variable dependiente (Y): área afectada del deslizamiento
Y_original <- datos_modelo$area_affected_m2

3.- Tabla de pares de valores original

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

n_original <- nrow(
  tabla_original
)

cat(
  "Tamaño muestral del modelo =",
  n_original
)
## Tamaño muestral del modelo = 11033
knitr::kable(
  head(
    tabla_original,
    20
  ),
  digits = 2,
  col.names = c(
    "Precipitación acumulada en 7 días (mm)",
    "Área afectada (m²)"
  ),
  caption = paste(
    "Tabla Nro. 1. Pares de valores originales de la precipitación",
    "acumulada en 7 días y el área afectada del deslizamiento"
  )
)
Tabla Nro. 1. Pares de valores originales de la precipitación acumulada en 7 días y el área afectada del deslizamiento
Precipitación acumulada en 7 días (mm) Área afectada (m²)
60.21 1879.8
115.77 2183.0
47.14 784.5
52.52 1815.4
12.51 1399.4
84.97 2568.7
121.28 2725.3
83.39 2167.9
93.87 1983.9
105.27 1835.7
44.63 413.2
60.71 1609.5
9.67 1917.7
67.48 1624.4
55.84 2450.9
31.55 971.8
33.14 885.1
47.54 881.4
23.37 2208.4
122.71 4071.5

4.- Gráfica original de dispersión

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

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 precipitación acumulada en 7 días\n",
      "y el área afectada del deslizamiento"
    ),
    x = "Precipitación acumulada en 7 días (mm)",
    y = "Área afectada del deslizamiento (m²)"
  ) +
  scale_x_continuous(
    limits = c(
      0,
      x_max_original
    ),
    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,
      x_max_original
    ),
    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

Promediación de los datos:

Inicialmente se eliminan los valores faltantes. Después, como el área afectada es una variable continua, se agrupa en niveles de 100 m² y se calcula la precipitación promedio correspondiente a cada nivel.

datos_modelo$area_nivel_m2 <- round(
  datos_modelo$area_affected_m2 / 100
) * 100

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

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

# Variable independiente (X): precipitación promedio
X <- datos_prom$precipitation_7d_mm

# Variable dependiente (Y): nivel de área afectada
Y <- datos_prom$area_nivel_m2

5.1.- Tabla pares de valores simplificada

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

n_modelo <- nrow(
  tabla_simplificada
)

cat(
  "Tamaño muestral del modelo =",
  n_modelo
)
## Tamaño muestral del modelo = 164
knitr::kable(
  head(
    tabla_simplificada,
    20
  ),
  digits = 2,
  col.names = c(
    "Precipitación promedio en 7 días (mm)",
    "Nivel de área afectada (m²)"
  ),
  caption = paste(
    "Tabla Nro. 2. Pares de valores simplificados de la precipitación",
    "acumulada en 7 días y el área afectada del deslizamiento"
  )
)
Tabla Nro. 2. Pares de valores simplificados de la precipitación acumulada en 7 días y el área afectada del deslizamiento
Precipitación promedio en 7 días (mm) Nivel de área afectada (m²)
0.50 200
27.14 300
28.07 400
29.85 500
32.72 600
37.44 700
40.23 800
43.09 900
45.69 1000
47.41 1100
50.81 1200
56.27 1300
53.15 1400
58.53 1500
57.73 1600
62.62 1700
65.84 1800
66.03 1900
68.69 2000
71.36 2100

El tamaño muestral disminuye debido al tratamiento de los datos antes de construir el modelo.

5.2.- Gráfica simplificada de dispersión

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

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

ggplot(
  tabla_simplificada,
  aes(
    x = X,
    y = Y
  )
) +
  geom_point(
    color = "#00A6D6",
    alpha = 0.70,
    size = 2
  ) +
  labs(
    title = "Gráfica Nro. 2",
    subtitle = paste0(
      "Diagrama simplificado de dispersión entre la precipitación acumulada en 7 días\n",
      "y el área afectada del deslizamiento"
    ),
    x = "Precipitación acumulada en 7 días (mm)",
    y = "Área afectada del deslizamiento (m²)"
  ) +
  scale_x_continuous(
    limits = c(
      0,
      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(
      0,
      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 muestra una curva ascendente, lo que sugiere un modelo logarítmico. El área afectada aumenta a medida que se incrementa la precipitación acumulada en 7 días, indicando una relación no lineal.

Modelo logarítmico general:

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

Modelo logarítmico aplicado al estudio:

\[ \text{Área afectada} = a+ b\ln( \text{Precipitación}_{7d} ) \]

7.- Parámetros

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

a_est <- coef(
  modelo_log
)[1]

b_est <- coef(
  modelo_log
)[2]

a_est
## (Intercept) 
##   -28241.19
b_est
##   log(X) 
## 7754.239

Ecuación logarítmica:

\[ \widehat{Y} = -28241.1891 + 7754.2395 \ln(X) \]

El intercepto \(a\) representa el valor teórico del área afectada y el coeficiente \(b\) determina el cambio asociado al logaritmo natural de la precipitación.

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

Aunque el modelo presenta una relación logarítmica, se utiliza la función lm() porque puede expresarse como un modelo lineal respecto a sus parámetros.

Modelo logarítmico general:

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

Aplicando la transformación:

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

Se obtiene el modelo lineal:

\[ Y= \beta_0+ \beta_1X_1 \]

Esta ecuación puede estimarse mediante:

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

Finalmente, los coeficientes permiten reconstruir el modelo logarítmico:

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

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

x_inicio_curva <- max(
  0.01,
  min(
    X,
    na.rm = TRUE
  ) * 0.20
)

X_curva <- seq(
  x_inicio_curva,
  x_max_simplificado,
  length.out = 500
)

Y_curva <- a_est +
  b_est * log(
    X_curva
  )

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

ggplot() +
  geom_point(
    data = tabla_simplificada,
    aes(
      x = X,
      y = Y
    ),
    color = "#00A6D6",
    alpha = 0.70,
    size = 2
  ) +
  geom_line(
    data = curva_logaritmica,
    aes(
      x = X,
      y = Y
    ),
    color = "red",
    linewidth = 1.5
  ) +
  labs(
    title = "Gráfica Nro. 3",
    subtitle = paste0(
      "Modelo logarítmico entre la precipitación acumulada en 7 días\n",
      "y el área afectada del deslizamiento"
    ),
    x = "Precipitación acumulada en 7 días (mm)",
    y = "Área afectada del deslizamiento (m²)"
  ) +
  scale_x_continuous(
    limits = c(
      0,
      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(
      0,
      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(),
    legend.position = "none"
  )

9.- Test de bondad

r <- cor(
  X,
  Y
) * 100

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

El coeficiente de correlación es de 83.68 %, lo que evidencia una relación positiva fuerte entre la precipitación acumulada y el área afectada.

10.- Restricciones

x_limite_fisico <- exp(
  -a_est / b_est
)

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

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

Dominios:

\[ D_X= \left\{ X\in\mathbb{R}: X>0 \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: Sí.

El modelo presenta una restricción matemática porque el logaritmo natural solo está definido para \(X>0\). También presenta una restricción física: cuando la precipitación es menor que 38.17 mm, el modelo genera áreas afectadas negativas, lo cual no es coherente con la naturaleza de la variable.

Por tanto, el modelo es físicamente válido para:

\[ X \geq 38.17 \text{ mm} \]

y debe utilizarse preferentemente dentro del intervalo de calibración:

\[ 0.5 \leq X \leq 507.07 \]

11.- Estimación

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

if (
  precipitacion_cantidad <= 0
) {
  stop(
    "La precipitación debe ser mayor que cero."
  )
}

area_cantidad <- a_est +
  b_est * log(
    precipitacion_cantidad
  )

# Pregunta de porcentaje
incremento_precipitacion <- 0.10

precipitacion_final <- precipitacion_cantidad *
  (
    1 +
      incremento_precipitacion
  )

area_inicial <- a_est +
  b_est * log(
    precipitacion_cantidad
  )

area_final <- a_est +
  b_est * log(
    precipitacion_final
  )

incremento_area_porcentaje <- (
  (
    area_final -
      area_inicial
  ) /
    area_inicial
) * 100

Pregunta de cantidad

¿Cuál es el área afectada esperada cuando la precipitación acumulada en 7 días es de 157.64 mm?

Resultado: 10997.85 m²

Pregunta de porcentaje

¿En qué porcentaje aumenta el área afectada esperada cuando la precipitación acumulada en 7 días aumenta un 10 %?

Resultado: 6.72 %

12.- Conclusión

Entre la precipitación acumulada en 7 días y el área afectada del deslizamiento existe una relación de tipo logarítmico, representada por el modelo \(\widehat{Y}=-28241.1891 + 7754.2395\ln(X)\). El modelo presenta una restricción matemática porque solo está definido para valores de precipitación mayores que cero y una restricción física porque únicamente genera áreas afectadas no negativas para valores de precipitación mayores o iguales a 38.17 mm.