Análisis Multivariable de Sedimentos Marinos

1. LIBRERÍAS

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

La librería dplyr se emplea para filtrar, ordenar y transformar los datos sin alterar la lógica del análisis. La librería gt permite presentar los resultados en tablas formateadas para el informe HTML. Finalmente, knitr permite integrar código, resultados, tablas y gráficas durante la generación del documento mediante Knit HTML.

2. CARGA DE DATOS

# Importar la base de datos
# Cambia la ruta si tu archivo está en otra ubicación.
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

3. SELECCIÓN DE VARIABLES

En este análisis, PHI_10 se define como la variable explicativa (X) y CLAY_PCT como la variable respuesta (Y). Desde el punto de vista sedimentológico, PHI_10 es un percentil granulométrico expresado en la escala phi; por tanto, resume una característica de la distribución de tamaños de grano del sedimento. Los valores de esta escala se relacionan con la finura del material: una variación hacia tamaños más finos refleja condiciones de depósito en las que puede aumentar la presencia relativa de fracciones arcillosas.

El porcentaje de arcilla (CLAY_PCT) cuantifica directamente la proporción de la fracción más fina de la muestra. Por ello, existe una base física para emplear PHI_10 como variable capaz de explicar parte de la variación de CLAY_PCT: ambos indicadores describen, desde perspectivas complementarias, la granulometría del sedimento. La relación propuesta no implica que PHI_10 sea la única causa del contenido de arcilla; representa una relación explicativa empírica en la cual el comportamiento granulométrico resume procesos de transporte, selección y depositación que se asocian con la proporción de arcilla.

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 (%)

4. TABLA DE PARES DE VALORES

La tabla se construye exclusivamente con PHI_10 y CLAY_PCT, que son las dos variables definidas para el estudio.Se informa el tamaño muestral y se presentan únicamente las primeras 20 observaciones originales.

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

n_original <- nrow(TPV)

tabla_tpv_original <- head(TPV, 20)

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

tabla_tpv_original %>%
  gt() %>%
  cols_label(
    Nro = "N°",
    PHI_10 = "PHI_10",
    CLAY_PCT = "CLAY_PCT (%)"
  ) %>%
  tab_header(
    title = md("**Tabla N°3. Primeras 20 observaciones de los pares PHI_10 y CLAY_PCT.**"),
    subtitle = paste("Tamaño muestral original:", n_original, "observaciones")
  ) %>%
  fmt_number(columns = c(PHI_10, CLAY_PCT), decimals = 2) %>%
  cols_align(align = "center") %>%
  tab_options(
    table.width = pct(85),
    column_labels.font.weight = "bold"
  )

Tabla N°3. Primeras 20 observaciones de los pares PHI_10 y CLAY_PCT.
Tamaño muestral original: 27784 observaciones
N°	PHI_10	CLAY_PCT (%)
1	4.96	0.00
2	4.07	0.00
3	9.46	0.00
4	9.35	0.00
5	12.40	0.00
6	3.70	12.70
7	6.00	19.60
8	8.40	0.00
9	4.96	13.74
10	5.32	0.00
11	2.36	0.00
12	4.41	35.35
13	3.10	16.20
14	5.90	19.70
15	10.20	29.70
16	4.96	0.00
17	NA	0.00
18	6.45	0.00
19	0.74	0.00
20	10.62	0.00

5. GRÁFICA DE DISPERSIÓN

La siguiente gráfica representa todos los datos originales, sin ningún tratamiento previo ni muestreo para visualización. Por tanto, se construye antes de eliminar valores faltantes, valores físicamente imposibles, duplicados u observaciones atípicas.

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

plot(
  TPV$PHI_10,
  TPV$CLAY_PCT,
  pch = 16,
  cex = 0.7,
  col = rgb(0,0,1,0.35),
  main = "Gráfica N°1. Datos originales sin tratamiento",
  xlab = "PHI_10",
  ylab = "CLAY_PCT (%)"
)

grid()
box(which = "outer")

En los datos originales se observa una gran dispersión. La nube de puntos es compleja, pues incorpora registros físicamente inconsistentes, duplicados y observaciones extremas; los valores faltantes no se representan gráficamente, aunque aún pertenecen a la base original. En consecuencia, en esta etapa todavía no es posible seleccionar con seguridad un modelo matemático que represente la relación entre las variables.

6. CONJETURA

6.1 Tratamiento de datos

El tratamiento se aplica exclusivamente a los pares PHI_10 y CLAY_PCT. Primero se identifican y eliminan los valores faltantes, ya que una observación incompleta no permite evaluar la relación entre ambas variables. Después se excluyen los valores infinitos para asegurar que los cálculos numéricos sean definidos. A continuación, se eliminan los valores físicamente imposibles: PHI_10 debe ser no negativo , mientras que CLAY_PCT debe permanecer entre 0 % y 100 % por su naturaleza de porcentaje.

Posteriormente se suprimen los duplicados exactos para que una misma pareja de valores no tenga una influencia repetida en el ajuste. Finalmente, se realiza una depuración conservadora de valores atípicos. Se ajusta primero un modelo polinomial base de grado 5 y se calculan sus residuos studentizados. Se eliminan únicamente las observaciones con valor absoluto del residuo studentizado mayor a 3. Este criterio identifica registros extremadamente alejados del comportamiento estimado, preservando la mayor parte de la variabilidad natural de los datos geológicos. No se modifica el algoritmo de depuración ni el criterio de corte.

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

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

TPV_sin_inf <- TPV_sin_na[
  is.finite(TPV_sin_na$PHI_10) &
    is.finite(TPV_sin_na$CLAY_PCT),
]

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

TPV_limpia <- TPV_sin_inf %>%
  filter(
    PHI_10 >= 0,
    CLAY_PCT >= 0,
    CLAY_PCT <= 100
  )

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

TPV_limpia <- TPV_limpia %>%
  distinct(PHI_10, CLAY_PCT, .keep_all = TRUE)

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

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

TPV_limpia$residuo_student <- rstudent(modelo_base)

TPV_final <- TPV_limpia %>%
  filter(abs(residuo_student) <= 3) %>%
  select(-residuo_student) %>%
  arrange(PHI_10)

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

tabla_depuracion <- data.frame(
  Etapa = c(
    " datos originales",
    "Filas con valores faltantes eliminadas",
    "Valores infinitos eliminados",
    "Valores físicamente imposibles eliminados",
    "Duplicados exactos eliminados",
    "Outliers extremos eliminados por residuo studentizado > 3",
    " datos finales utilizados 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°4. Resumen del tratamiento conservador aplicado a PHI_10 y CLAY_PCT.**"),
    subtitle = paste0("Registros antes: ", n_original, "; registros después: ", n_final,
                      "; porcentaje total eliminado: ", round(porcentaje_eliminado, 2), "%")
  ) %>%
  cols_align(align = "center") %>%
  tab_options(
    table.width = pct(90),
    column_labels.font.weight = "bold"
  )

Tabla N°4. Resumen del tratamiento conservador aplicado a PHI_10 y CLAY_PCT.
Registros antes: 27784; registros después: 18362; porcentaje total eliminado: 33.91%
Etapa del tratamiento de datos	Número de registros
datos originales	27784
Filas con valores faltantes eliminadas	1300
Valores infinitos eliminados	0
Valores físicamente imposibles eliminados	0
Duplicados exactos eliminados	7554
Outliers extremos eliminados por residuo studentizado > 3	568
datos finales utilizados en la regresión	18362

En total, el procedimiento reduce la base de 27784 a 18362 registros, lo que corresponde a un 33.91 % de registros eliminados. Específicamente, la depuración por residuos studentizados elimina 568 observaciones (3 % de los registros disponibles antes de esa etapa).

6.2 Nueva gráfica

x <- TPV_final$PHI_10
y <- TPV_final$CLAY_PCT

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

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°2. Dispersión tratada entre PHI_10 y CLAY_PCT",
  xlab = "PHI_10",
  ylab = "CLAY_PCT (%)"
)

grid()
box(which = "outer")

Tras el tratamiento, la tendencia general puede observarse con mayor claridad. La reducción de registros incompletos, físicamente imposibles, repetidos y extremadamente atípicos permite examinar la relación predominante sin eliminar la dispersión natural propia de las muestras de sedimentos marinos.

6.3 Nueva conjetura

La nueva nube de puntos es compatible con una regresión polinomial porque la respuesta esperada de CLAY_PCT no parece variar a una tasa constante a lo largo de PHI_10. Se aprecia una trayectoria curvilínea, con cambios en la pendiente que una recta no representaría adecuadamente. Por esta razón se plantea una familia de modelos polinomiales y se conserva el mecanismo original de comparación y selección de grado.

7. CÁLCULO DE PARÁMETROS

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

for(grado in 2:6){
  formula_modelo <- as.formula(
    paste("CLAY_PCT ~", 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(datos_modelo$PHI_10, datos_modelo$CLAY_PCT),
      Porcentaje_Pearson = cor(datos_modelo$PHI_10, datos_modelo$CLAY_PCT) * 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.005) %>%
  arrange(Grado)

mejor_grado <- candidatos$Grado[1]
modelo_final <- modelos[[paste0("Grado_", mejor_grado)]]

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 = 2) %>%
  cols_align(align = "center") %>%
  tab_header(
    title = md("**Tabla N°5. Coeficientes estimados del modelo de regresión polinomial.**")
  ) %>%
  tab_options(
    table.width = pct(80),
    column_labels.font.weight = "bold"
  )

Tabla N°5. Coeficientes estimados del modelo de regresión polinomial.
Parámetro	Término	Valor
Coeficiente 0	Intercepto	1.10
Coeficiente 1	X^1	2.83
Coeficiente 2	X^2	−0.02

Los coeficientes expresan el aporte estimado de cada potencia de PHI_10 al porcentaje de arcilla, manteniendo los demás términos del polinomio. El intercepto representa la respuesta estimada cuando PHI_10 es cero dentro de la formulación matemática; los demás coeficientes describen la curvatura de la relación.

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

for(i in 2:length(coeficientes)){
  signo <- ifelse(coeficientes[i] >= 0, " + ", " - ")
  ecuacion <- paste0(
    ecuacion,
    signo,
    abs(round(coeficientes[i],2)),
    "·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")

8. REALIDAD Y MODELO

x_modelo <- seq(
  min(x),
  max(x),
  length.out = 500
)

datos_pred <- data.frame(
  PHI_10 = 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°3. 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"
)

La superposición permite contrastar la realidad observada con el patrón estimado. Los puntos conservan una dispersión alrededor de la curva, como es esperable en datos naturales, mientras que la línea roja resume la tendencia promedio que el modelo polinomial atribuye a la relación entre PHI_10 y CLAY_PCT.

9. TEST

pearson_final <- cor.test(
  datos_modelo$PHI_10,
  datos_modelo$CLAY_PCT,
  method = "pearson"
)

r <- as.numeric(pearson_final$estimate)
porcentaje_pearson <- r * 100
R2 <- summary(modelo_final)$r.squared
R2_ajustado <- summary(modelo_final)$adj.r.squared
AIC_final <- AIC(modelo_final)
BIC_final <- BIC(modelo_final)
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²)"
  ),
  Valor = c(
    round(r,4),
    round(porcentaje_pearson,2),
    round(R2,4)
  )
)

tabla_tests %>%
  gt() %>%
  cols_label(
    Indicador = "Indicador estadístico",
    Valor = "Valor"
  ) %>%
  cols_align(align = "center") %>%
  tab_header(
    title = md("**Tabla N°6. Indicadores estadísticos del modelo polinomial .**")
  ) %>%
  tab_options(
    table.width = pct(85),
    column_labels.font.weight = "bold"
  )

Tabla N°6. Indicadores estadísticos del modelo polinomial .
Indicador estadístico	Valor
Coeficiente de correlación de Pearson (r)	0.9195
Porcentaje de Pearson (%)	91.9500
Coeficiente de determinación (R²)	0.8562

10. RESTRICCIONES

El modelo es válido únicamente dentro del intervalo de PHI_10 observado después del tratamiento de datos. Fuera de ese intervalo, el uso de la ecuación constituye una extrapolación y puede producir estimaciones inestables o físicamente poco plausibles, especialmente porque los polinomios pueden cambiar su comportamiento de forma marcada en los extremos.

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,66,"El modelo debe utilizarse únicamente",cex = 1.15)
text(50,58,"dentro del rango observado de PHI_10:",cex = 1.15)

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

text(50,33,"Modelo empírico: no demuestra causalidad.",cex = 1.05)
text(50,25,"No se recomienda realizar extrapolaciones.",cex = 1.05)
text(50,16,"CLAY_PCT debe interpretarse entre 0 % y 100 %.",cex = 0.95,col = "gray40")

11. ESTIMACIÓN

La estimación se conserva para un valor de PHI_10 igual a 8 y se presenta mediante el mismo procedimiento de predicción del modelo final.

x_estimacion <- 8

nuevo_dato <- data.frame(
  PHI_10 = 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,2), " %"),
  cex = 2,
  font = 2,
  col = "#C0392B"
)

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

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.41 %. Esta estimación debe interpretarse como el valor esperado según la relación ajustada y es válida únicamente cuando el valor de entrada se encuentra dentro del intervalo observado en los datos tratados.

12. CONCLUSIÓN

Entre el percentil granulométrico PHI_10 y el porcentaje de arcilla CLAY_PCT se observa una relación tipo polinomial de grado 2, representada por el modelo CLAY_PCT = 1.1 + 2.83·PHI_10^1 - 0.02·PHI_10^2, donde x corresponde al percentil granulométrico PHI_10 e y corresponde al porcentaje de arcilla CLAY_PCT (%). El modelo solo es confiable dentro del rango observado de PHI_10.

Ejemplo: cuando PHI_10 es igual a 8, el modelo predice un porcentaje de arcilla aproximado de 22.41 %.