Proyecto estrictamente académico (Estadística y Probabilidades, FIGEMPA – UCE). Datos: EarthChem PetDB (CC BY 4.0), descargados el 13/10/2026. Web del proyecto: https://crisrumiguano-art.github.io/geoquimica-earthchem-ml/

1. Carga de datos

Se usa la base completa (159 760 análisis). Para cada par se toman solo las filas que tienen los dos datos: no se inventa ningún valor que falte.

knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE, fig.align = "center")
if (!file.exists("data/earthchem_limpio.csv.gz")) {
  dir.create("data", showWarnings = FALSE)
  download.file("https://github.com/crisrumiguano-art/geoquimica-earthchem-ml/raw/main/data/earthchem_limpio.csv.gz", "data/earthchem_limpio.csv.gz", mode = "wb")
}
base <- read.csv(gzfile("data/earthchem_limpio.csv.gz", encoding = "UTF-8"), check.names = FALSE)
dim(base)
## [1] 159760     75

2. Asignación de variables

x = MgO (wt%) (variable independiente) e y = SiO2 (wt%) (variable dependiente). Es la misma diferenciación magmática vista desde el otro lado (como en la plantilla del curso): las rocas con poco magnesio son las más evolucionadas y ricas en sílice.

x <- base[["MgO (wt%)"]]
y <- base[["SiO2 (wt%)"]]
completos <- !is.na(x) & !is.na(y) & x > 0 & y > 0    # filas con los dos datos
x <- x[completos]
y <- y[completos]
silicatada <- y >= 35 & y <= 80                    # rocas silicatadas (restricción del dominio)
x <- x[silicatada]
y <- y[silicatada]
length(x)                                              # número de pares
## [1] 53630

3. Tabla de pares de valores

knitr::kable(head(data.frame(x, y), 10), caption = "Tabla N° 1: primeros 10 pares de valores (MgO (wt%), SiO2 (wt%))")
Tabla N° 1: primeros 10 pares de valores (MgO (wt%), SiO2 (wt%))
x y
7.88 50.84
8.77 48.78
1.65 64.40
5.32 49.20
7.59 49.28
7.57 50.64
6.34 48.48
44.02 45.54
6.75 50.81
4.20 51.41 
knitr::kable(data.frame(variable = c("MgO (wt%)", "SiO2 (wt%)"), minimo = c(min(x), min(y)), media = c(mean(x), mean(y)),
                        maximo = c(max(x), max(y))), digits = 2, caption = "Tabla N° 2: resumen de los pares")
Tabla N° 2: resumen de los pares
variable minimo media maximo
MgO (wt%) 0 11.17 54.76
SiO2 (wt%) 35 52.57 79.83

4. Diagrama de dispersión

plot(x, y, pch = 16, cex = 0.4, col = rgb(0.22, 0.53, 0.9, 0.25), las = 1,
     main = "Gráfica N° 1: Diagrama de dispersión", xlab = "MgO (wt%)", ylab = "SiO2 (wt%)")

La sílice baja muy rápido cuando el MgO sube de 0 a 10 % y después casi se estabiliza: una curva que se aplana.

5. Conjetura del modelo

Una subida o bajada rápida que después se aplana es la forma del modelo logarítmico, y = a + b·ln x. Antes de elegir se probaron los cinco modelos con este par (prueba y error): el coeficiente de correlación de cada uno sale de la recta ajustada a las variables transformadas.

r_modelos <- c(lineal = cor(x, y),                     # y = a + b·x
               potencial = cor(log(x), log(y)),        # y = a·x^b      -> log y = log a + b·log x
               exponencial = cor(x, log(y)),           # y = a·e^(b·x)  -> log y = log a + b·x
               logaritmico = cor(log(x), y),           # y = a + b·ln x
               polinomico = sqrt(summary(lm(y ~ x + I(x^2) + I(x^3)))$r.squared))  # y = a + b·x + c·x² + d·x³
knitr::kable(data.frame(modelo = names(r_modelos), r = round(r_modelos, 3), r2_porcentaje = round(r_modelos^2 * 100, 1)),
             row.names = FALSE, caption = "Tabla N° 3: correlación de cada modelo con este par")
Tabla N° 3: correlación de cada modelo con este par
modelo r r2_porcentaje
lineal -0.587 34.4
potencial -0.873 76.2
exponencial -0.627 39.3
logaritmico -0.879 77.3
polinomico 0.853 72.8 
modelo <- lm(y ~ log(x))
a <- coef(modelo)[[1]]
b <- coef(modelo)[[2]]
c(a = a, b = b)
##         a         b 
## 61.235874 -5.309858
plot(x, y, pch = 16, cex = 0.4, col = rgb(0.22, 0.53, 0.9, 0.25), las = 1,
     main = "Gráfica N° 2: Regresión logarítmico", xlab = "MgO (wt%)", ylab = "SiO2 (wt%)")
curve(a + b * log(x), from = min(x), to = max(x), add = TRUE, col = "#d95926", lwd = 3)
legend("topright", c("datos", "y = a + b·ln x"), col = c(rgb(0.22, 0.53, 0.9), "#d95926"),
       pch = c(16, NA), lwd = c(NA, 3), bty = "n")

6. Test de aprobación

El coeficiente de correlación de Pearson (r) mide qué tan fuerte es la relación: |r| cerca de 1 es fuerte. El coeficiente de determinación (r² en %) dice qué parte de la variación de y explica el modelo.

r <- cor(log(x), y)
r2 <- r^2 * 100
knitr::kable(data.frame(modelo = "logarítmico", r = round(r, 4), r2_porcentaje = round(r2, 2),
                        fuerza = ifelse(abs(r) >= 0.7, "fuerte", ifelse(abs(r) >= 0.4, "moderada", "débil"))),
             caption = "Tabla N° 4: test de aprobación")
Tabla N° 4: test de aprobación
modelo r r2_porcentaje fuerza
logarítmico -0.879 77.27 fuerte 
residuo <- y - (a + b * log(x))                                   # lo que el modelo no explica
plot(x, residuo, pch = 16, cex = 0.4, col = rgb(0.22, 0.53, 0.9, 0.25), las = 1,
     ylim = quantile(residuo, c(0.0025, 0.9975)),           # sin el 0,5 % más extremo, para que se vea
     main = "Gráfica N° 3: Residuos del modelo", xlab = "MgO (wt%)", ylab = "Residuo (y real - y del modelo)")
abline(h = 0, col = "#d95926", lwd = 2)

7. Restricciones

Se usan solo rocas silicatadas (SiO2 entre 35 y 80 %). El logaritmo exige x > 0. Además y debe ser positivo: con b negativo, y > 0 solo si x < e^(−a/b). Fuera de ese valor el modelo daría sílice negativa, imposible.

c(x_debe_ser_mayor_que = 0, x_maximo_valido = exp(-a / b), x_max_datos = max(x))
## x_debe_ser_mayor_que      x_maximo_valido          x_max_datos 
##                 0.00            101975.59                54.76

8. Estimación de pronósticos

Si una roca tiene 3 % de MgO, ¿cuánta sílice se espera?

x0 <- 3
y0 <- a + b * log(x0)
y0
## [1] 55.4024

9. Conclusión

Entre MgO (wt%) e SiO2 (wt%) existe una relación que se representa con un modelo logarítmico: y = 61.236 + (-5.31)·ln x, donde x es MgO (wt%) e y es SiO2 (wt%). La correlación es fuerte (r = -0.879): SiO2 (wt%) depende de MgO (wt%) en un 77.3 % y el 22.7 % restante se debe a otros factores que este modelo no considera. Si MgO (wt%) vale 3, se espera SiO2 (wt%) ≈ 55.4. El MgO sirve para estimar qué tan evolucionada está una roca.