1. CARGA DE DATOS Y LIBRERÍAS

#===============================================================================
# 1. CARGA DE DATOS Y LIBRERÍAS
# Análisis de regresión polinomial simple aplicado a sedimentos marinos
#===============================================================================

library(dplyr)
library(gt)
library(knitr)

# Cambia esta ruta según la ubicación real de tu archivo CSV.
datos <- read.csv(
  "~/ESTADISTICA/dataset_geologico_limpio_80.csv",
  header = TRUE,
  sep = ",",
  dec = ".",
  stringsAsFactors = FALSE
)

resumen_bd <- data.frame(
  Descripción = c("Número de observaciones", "Número de variables"),
  Valor = c(nrow(datos), ncol(datos))
)

resumen_bd %>%
  gt() %>%
  tab_header(title = md("**Tabla N°1. Resumen de la Base de Datos**")) %>%
  cols_align(align = "center") %>%
  tab_options(
    table.width = pct(60),
    column_labels.font.weight = "bold"
  )
Tabla N°1. Resumen de la Base de Datos
Descripción Valor
Número de observaciones 27784
Número de variables 58

2. SELECCIÓN DE VARIABLES

Para el desarrollo del modelo de regresión polinomial simple se seleccionaron dos variables cuantitativas del conjunto de datos de sedimentos marinos. La variable independiente corresponde a PHI_10, mientras que la variable dependiente corresponde a CLAY_PCT. Esta selección se fundamenta en que PHI_10 presentó una relación positiva fuerte con el porcentaje de arcilla, por lo que resulta adecuada para construir un modelo de estimación.

tabla_variables <- data.frame(
  Rol = c("Variable Independiente (X)", "Variable Dependiente (Y)"),
  Variable = c("PHI_10", "CLAY_PCT"),
  Descripción = c("Percentil granulométrico PHI_10", "Porcentaje de arcilla (%)")
)

tabla_variables %>%
  gt() %>%
  tab_header(
    title = md("**Tabla N°2. Variables seleccionadas para el modelo de regresión polinomial aplicado al análisis de sedimentos marinos.**")
  ) %>%
  cols_align(align = "center") %>%
  tab_options(
    table.width = pct(85),
    column_labels.font.weight = "bold"
  )
Tabla N°2. Variables seleccionadas para el modelo de regresión polinomial aplicado al análisis de sedimentos marinos.
Rol Variable Descripción
Variable Independiente (X) PHI_10 Percentil granulométrico PHI_10
Variable Dependiente (Y) CLAY_PCT Porcentaje de arcilla (%)

3. TABLA DE PARES DE VALORES

3.1 Construcción de la tabla de pares de valores

TPV <- data.frame(
  x = as.numeric(datos$PHI_10),
  y = as.numeric(datos$CLAY_PCT)
)

n_original <- nrow(TPV)

3.2 Identificación y eliminación de valores faltantes

na_x <- sum(is.na(TPV$x))
na_y <- sum(is.na(TPV$y))
n_filas_na <- sum(!complete.cases(TPV))

TPV_sin_na <- na.omit(TPV)
n_sin_na <- nrow(TPV_sin_na)

3.3 Eliminación de valores infinitos

TPV_sin_inf <- TPV_sin_na[
  is.finite(TPV_sin_na$x) &
    is.finite(TPV_sin_na$y),
]

n_sin_inf <- nrow(TPV_sin_inf)
n_inf_eliminados <- n_sin_na - n_sin_inf

3.4 Eliminación de valores inconsistentes

TPV_limpia <- TPV_sin_inf %>%
  filter(
    x >= 0,
    y >= 0,
    y <= 100
  )

n_sin_inconsistencias <- nrow(TPV_limpia)
n_inconsistentes <- n_sin_inf - n_sin_inconsistencias

3.5 Eliminación de duplicados exactos

TPV_limpia <- TPV_limpia %>%
  distinct(x, y, .keep_all = TRUE)

n_sin_duplicados <- nrow(TPV_limpia)
n_duplicados <- n_sin_inconsistencias - n_sin_duplicados

3.6 Eliminación de valores atípicos mediante residuo studentizado

Para mejorar la calidad del ajuste sin modificar artificialmente la tendencia general, se eliminaron únicamente los valores atípicos fuertes detectados mediante el residuo studentizado. Se consideraron atípicos los registros cuyo valor absoluto del residuo studentizado fue mayor a 2.5.

modelo_base <- lm(
  y ~ x + I(x^2) + I(x^3) + I(x^4) + I(x^5),
  data = TPV_limpia
)

TPV_limpia$residuo_student <- rstudent(modelo_base)

TPV_final <- TPV_limpia %>%
  filter(abs(residuo_student) <= 2.5) %>%
  select(-residuo_student) %>%
  arrange(x)

n_final <- nrow(TPV_final)
n_outliers <- n_sin_duplicados - n_final
porcentaje_outliers <- (n_outliers / n_sin_duplicados) * 100

3.7 Resumen del tratamiento de datos

tabla_depuracion <- data.frame(
  Etapa = c(
    "Base de datos original",
    "Filas con valores faltantes eliminadas",
    "Valores infinitos eliminados",
    "Valores inconsistentes eliminados",
    "Duplicados exactos eliminados",
    "Valores atípicos eliminados por residuo studentizado",
    "Base final utilizada en la regresión"
  ),
  Registros = c(
    n_original,
    n_filas_na,
    n_inf_eliminados,
    n_inconsistentes,
    n_duplicados,
    n_outliers,
    n_final
  )
)

tabla_depuracion %>%
  gt() %>%
  cols_label(
    Etapa = "Etapa del tratamiento de datos",
    Registros = "Número de registros"
  ) %>%
  tab_header(
    title = md("**Tabla N°3. Tratamiento de datos aplicado a PHI_10 y CLAY_PCT.**")
  ) %>%
  cols_align(align = "center") %>%
  tab_options(
    table.width = pct(90),
    column_labels.font.weight = "bold"
  )
Tabla N°3. Tratamiento de datos aplicado a PHI_10 y CLAY_PCT.
Etapa del tratamiento de datos Número de registros
Base de datos original 27784
Filas con valores faltantes eliminadas 1300
Valores infinitos eliminados 0
Valores inconsistentes eliminados 0
Duplicados exactos eliminados 7554
Valores atípicos eliminados por residuo studentizado 853
Base final utilizada en la regresión 18077

3.8 Tabla final de pares de valores

tabla_tpv_previa <- head(TPV_final, 20)

tabla_tpv_previa <- cbind(
  Nro = 1:nrow(tabla_tpv_previa),
  tabla_tpv_previa
)

tabla_tpv_previa %>%
  gt() %>%
  cols_label(
    Nro = "N°",
    x = "PHI_10",
    y = "CLAY_PCT (%)"
  ) %>%
  tab_header(
    title = md("**Tabla N°4. Pares de valores depurados utilizados para la regresión polinomial.**")
  ) %>%
  fmt_number(columns = c(x, y), decimals = 4) %>%
  cols_align(align = "center") %>%
  tab_options(
    table.width = pct(85),
    column_labels.font.weight = "bold"
  )
Tabla N°4. Pares de valores depurados utilizados para la regresión polinomial.
N° PHI_10 CLAY_PCT (%)
1 0.0000 3.2000
2 0.0000 11.2000
3 0.0000 0.2000
4 0.0000 0.6000
5 0.0000 1.7000
6 0.0000 2.0000
7 0.0000 0.1000
8 0.0000 1.5000
9 0.0000 4.4000
10 0.0000 9.9000
11 0.0000 4.0000
12 0.0000 2.5000
13 0.0000 11.5000
14 0.0000 0.0300
15 0.0000 0.0000
16 0.0000 2.7900
17 0.0000 1.0700
18 0.0000 1.4000
19 0.0000 0.0100
20 0.0000 1.8000

3.9 Definición final de las variables para el modelo

x <- TPV_final$x
y <- TPV_final$y

datos_modelo <- TPV_final %>%
  mutate(
    X1 = x,
    X2 = x^2,
    X3 = x^3,
    X4 = x^4,
    X5 = x^5,
    X6 = x^6
  )

4. DIAGRAMA DE DISPERSIÓN

El diagrama de dispersión permite visualizar la relación entre PHI_10 y CLAY_PCT. Debido al elevado número de observaciones, se utiliza una muestra aleatoria únicamente para mejorar la visualización. La regresión polinomial se realiza con la totalidad de los datos depurados.

set.seed(2026)

indice_visual <- sample(
  1:length(x),
  min(3000, length(x))
)

par(oma = c(1,1,1,1))

plot(
  x[indice_visual],
  y[indice_visual],
  pch = 16,
  cex = 0.7,
  col = rgb(0,0,1,0.35),
  main = "Gráfica N°1. Diagrama de dispersión entre PHI_10 y CLAY_PCT",
  xlab = "PHI_10",
  ylab = "CLAY_PCT (%)"
)

grid()
box(which = "outer")

La nube de puntos evidencia una relación positiva y curvilínea entre PHI_10 y CLAY_PCT. Esta tendencia permite plantear el uso de un modelo polinomial simple para representar el comportamiento general de los datos.

5. CONJETURA DEL MODELO

A partir del diagrama de dispersión se plantea que existe una relación no lineal entre PHI_10 y CLAY_PCT. Por ello, se propone ajustar modelos polinomiales de distinto grado y seleccionar aquel que presente buen ajuste estadístico sin generar complejidad innecesaria.

resultados_modelos <- data.frame()
modelos <- list()

for(grado in 2:6){
  
  formula_modelo <- as.formula(
    paste("y ~", paste(paste0("X", 1:grado), collapse = " + "))
  )
  
  modelo <- lm(formula_modelo, data = datos_modelo)
  modelos[[paste0("Grado_", grado)]] <- modelo
  
  resultados_modelos <- rbind(
    resultados_modelos,
    data.frame(
      Grado = grado,
      Pearson = cor(x, y),
      Porcentaje_Pearson = cor(x, y) * 100,
      R2 = summary(modelo)$r.squared,
      R2_Ajustado = summary(modelo)$adj.r.squared,
      AIC = AIC(modelo),
      BIC = BIC(modelo),
      Error_Residual = sigma(modelo)
    )
  )
}

resultados_modelos <- resultados_modelos %>%
  arrange(desc(R2_Ajustado), AIC, BIC)

mejor_r2_ajustado <- max(resultados_modelos$R2_Ajustado)

candidatos <- resultados_modelos %>%
  filter(R2_Ajustado >= mejor_r2_ajustado - 0.002) %>%
  arrange(Grado)

mejor_grado <- candidatos$Grado[1]
modelo_final <- modelos[[paste0("Grado_", mejor_grado)]]
x_modelo <- seq(
  min(x),
  max(x),
  length.out = 500
)

datos_pred <- data.frame(
  x = x_modelo,
  X1 = x_modelo,
  X2 = x_modelo^2,
  X3 = x_modelo^3,
  X4 = x_modelo^4,
  X5 = x_modelo^5,
  X6 = x_modelo^6
)

y_modelo <- predict(
  modelo_final,
  newdata = datos_pred
)

par(oma = c(1,1,1,1))

plot(
  x[indice_visual],
  y[indice_visual],
  pch = 16,
  cex = 0.7,
  col = rgb(0,0,1,0.35),
  main = paste("Gráfica N°2. Modelo polinomial grado", mejor_grado, "entre PHI_10 y CLAY_PCT"),
  xlab = "PHI_10",
  ylab = "CLAY_PCT (%)"
)

lines(
  x_modelo,
  y_modelo,
  col = "red",
  lwd = 3
)

grid()
box(which = "outer")

legend(
  "bottomright",
  legend = c("Datos observados", "Modelo polinomial"),
  col = c(rgb(0,0,1,0.35), "red"),
  pch = c(16, NA),
  lty = c(NA, 1),
  lwd = c(NA, 3),
  bty = "n"
)

6. PARÁMETROS DEL MODELO

Una vez comparados los modelos polinomiales, se selecciona el grado más bajo cuyo ajuste sea similar al mejor modelo. Este criterio evita utilizar modelos excesivamente complejos cuando la mejora estadística es mínima.

coeficientes <- coef(modelo_final)

terminos <- c("Intercepto", paste0("X^", 1:(length(coeficientes)-1)))

tabla_parametros <- data.frame(
  Parámetro = paste0("Coeficiente ", 0:(length(coeficientes)-1)),
  Término = terminos,
  Valor = coeficientes
)

tabla_parametros %>%
  gt() %>%
  fmt_number(columns = Valor, decimals = 6) %>%
  cols_align(align = "center") %>%
  tab_header(
    title = md("**Tabla N°6. Parámetros estimados del modelo de regresión polinomial.**")
  ) %>%
  tab_options(
    table.width = pct(80),
    column_labels.font.weight = "bold"
  )
Tabla N°6. Parámetros estimados del modelo de regresión polinomial.
Parámetro Término Valor
Coeficiente 0 Intercepto 0.854683
Coeficiente 1 X^1 2.901922
Coeficiente 2 X^2 −0.022751

6.1 Ecuación del modelo polinomial

ecuacion <- paste0("CLAY_PCT = ", round(coeficientes[1],4))

for(i in 2:length(coeficientes)){
  signo <- ifelse(coeficientes[i] >= 0, " + ", " - ")
  ecuacion <- paste0(
    ecuacion,
    signo,
    abs(round(coeficientes[i],4)),
    "·PHI_10^",
    i-1
  )
}

plot.new()
plot.window(xlim = c(0,100), ylim = c(0,100))

rect(5, 58, 95, 92, border = "#1F4E79", lwd = 3)
text(50, 87, "MODELO TEÓRICO", cex = 1.5, font = 2, col = "#1F4E79")
text(50, 72, "CLAY_PCT = a + bX + cX² + ...", cex = 1.25, font = 2, col = "#C0392B")

rect(5, 8, 95, 48, border = "#1F4E79", lwd = 3)
text(50, 43, "MODELO AJUSTADO", cex = 1.5, font = 2, col = "#1F4E79")
text(50, 25, ecuacion, cex = 0.85, font = 2, col = "#C0392B")

text(50, 4, "Modelo ajustado mediante el método de mínimos cuadrados.", cex = 0.85, col = "gray40")

7. TEST DE PEARSON

El test de Pearson permite cuantificar la intensidad de la relación lineal entre las variables analizadas. En este caso, se interpreta como una medida general de asociación positiva entre PHI_10 y CLAY_PCT después del tratamiento de datos.

r <- cor(x, y)
R2 <- summary(modelo_final)$r.squared
R2_ajustado <- summary(modelo_final)$adj.r.squared
error_residual <- sigma(modelo_final)

tabla_tests <- data.frame(
  Indicador = c(
    "Coeficiente de correlación de Pearson (r)",
    "Porcentaje de Pearson (%)",
    "Coeficiente de determinación (R²)",
    "R² ajustado",
    "Error residual",
    "Grado polinomial seleccionado"
  ),
  Valor = c(
    round(r,4),
    round(r*100,2),
    round(R2,4),
    round(R2_ajustado,4),
    round(error_residual,4),
    mejor_grado
  )
)

tabla_tests %>%
  gt() %>%
  cols_label(
    Indicador = "Indicador estadístico",
    Valor = "Valor"
  ) %>%
  cols_align(align = "center") %>%
  tab_header(
    title = md("**Tabla N°7. Indicadores estadísticos del modelo polinomial.**")
  ) %>%
  tab_options(
    table.width = pct(85),
    column_labels.font.weight = "bold"
  )
Tabla N°7. Indicadores estadísticos del modelo polinomial.
Indicador estadístico Valor
Coeficiente de correlación de Pearson (r) 0.9372
Porcentaje de Pearson (%) 93.7200
Coeficiente de determinación (R²) 0.8911
R² ajustado 0.8911
Error residual 5.9975
Grado polinomial seleccionado 2.0000

Interpretación:
El coeficiente de correlación de Pearson obtenido evidencia una relación positiva fuerte entre PHI_10 y CLAY_PCT. El coeficiente de determinación indica el porcentaje de variabilidad del porcentaje de arcilla explicado por el modelo polinomial ajustado.

8. RESTRICCIONES DEL MODELO

El modelo polinomial debe aplicarse únicamente dentro del intervalo observado de PHI_10, ya que fuera de este rango las predicciones corresponderían a extrapolaciones y podrían ser poco confiables.

plot.new()
plot.window(xlim = c(0,100), ylim = c(0,100))

rect(5,10,95,90,border = "#1F4E79",lwd = 3)

text(50,84,"RESTRICCIONES DEL MODELO",cex = 1.6,font = 2,col = "#1F4E79")
text(50,65,"El modelo debe utilizarse únicamente",cex = 1.15)
text(50,56,"dentro del rango observado de PHI_10:",cex = 1.15)

text(
  50,
  45,
  paste0(round(min(x),4), " ≤ PHI_10 ≤ ", round(max(x),4)),
  cex = 1.4,
  col = "#C0392B",
  font = 2
)

text(50,30,"No se recomienda realizar extrapolaciones",cex = 1.05)
text(50,22,"fuera del intervalo utilizado para ajustar el modelo.",cex = 1.05)

9. ESTIMACIÓN MEDIANTE EL MODELO

Una de las aplicaciones principales del modelo de regresión es estimar el valor esperado de la variable dependiente a partir de un valor conocido de la variable independiente. En este caso, se estima el porcentaje de arcilla para un valor específico de PHI_10.

x_estimacion <- 8

nuevo_dato <- data.frame(
  x = x_estimacion,
  X1 = x_estimacion,
  X2 = x_estimacion^2,
  X3 = x_estimacion^3,
  X4 = x_estimacion^4,
  X5 = x_estimacion^5,
  X6 = x_estimacion^6
)

y_estimacion <- predict(
  modelo_final,
  newdata = nuevo_dato
)

plot.new()
plot.window(xlim = c(0,100), ylim = c(0,100))

rect(5,15,95,85,border = "#1F4E79",lwd = 3)

text(50,80,"ESTIMACIÓN DEL MODELO POLINOMIAL",cex = 1.2,font = 2,col = "#1F4E79")
text(50,64,"Para una muestra de sedimento marino con",cex = 1.15)
text(50,56,paste0("PHI_10 igual a ", x_estimacion),cex = 1.15)
text(50,48,"el modelo estima un porcentaje de arcilla de:",cex = 1.15)

text(
  50,
  34,
  paste0(round(y_estimacion,4), " %"),
  cex = 2,
  font = 2,
  col = "#C0392B"
)

text(50,20,"Estimación obtenida mediante el modelo polinomial ajustado.",cex = 0.85,col = "gray40")

Interpretación:
Para una muestra de sedimento marino con un valor de PHI_10 igual a 8, el modelo polinomial estima un porcentaje de arcilla de 22.614 %. Esta estimación es válida únicamente si el valor utilizado se encuentra dentro del rango observado en los datos depurados.

10. CONCLUSIÓN

Entre PHI_10 y el porcentaje de arcilla CLAY_PCT existe una relación positiva de tipo polinomial, representada por el modelo seleccionado automáticamente a partir de la comparación entre grados 2, 3, 4, 5 y 6. Después del tratamiento de datos, se eliminaron valores faltantes, valores inconsistentes, duplicados exactos y valores atípicos detectados mediante residuo studentizado. El modelo final presentó un coeficiente de Pearson de 93.72 % y un coeficiente de determinación de 0.8911, lo que indica que explica aproximadamente el 89.11 % de la variabilidad observada en el porcentaje de arcilla. Por lo tanto, el modelo cumple el criterio mínimo de correlación mayor al 70 % y puede utilizarse como herramienta de estimación de CLAY_PCT a partir de PHI_10, siempre que las predicciones se realicen dentro del rango observado de los datos tratados.