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 = SiO2 (wt%) (variable independiente) e y = MgO (wt%) (variable dependiente). Cuando un magma se enfría, primero cristalizan el olivino y el piroxeno, ricos en magnesio: el líquido que queda tiene cada vez más sílice y menos magnesio (diferenciación magmática).

x <- base[["SiO2 (wt%)"]]
y <- base[["MgO (wt%)"]]
completos <- !is.na(x) & !is.na(y) & x > 0 & y > 0    # filas con los dos datos
x <- x[completos]
y <- y[completos]
silicatada <- x >= 35 & x <= 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 (SiO2 (wt%), MgO (wt%))")
Tabla N° 1: primeros 10 pares de valores (SiO2 (wt%), MgO (wt%))
x y
50.84 7.88
48.78 8.77
64.40 1.65
49.20 5.32
49.28 7.59
50.64 7.57
48.48 6.34
45.54 44.02
50.81 6.75
51.41 4.20 
knitr::kable(data.frame(variable = c("SiO2 (wt%)", "MgO (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
SiO2 (wt%) 35 52.57 79.83
MgO (wt%) 0 11.17 54.76

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 = "SiO2 (wt%)", ylab = "MgO (wt%)")

El MgO cae rápido al principio y después cada vez más lento: una curva que se aplana, no una recta.

5. Conjetura del modelo

Una caída que se aplana es la forma de un modelo exponencial, y = a·e^(b·x), con b negativo: log y = log a + b·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.879 77.3
logaritmico -0.627 39.3
polinomico 0.726 52.7 
modelo <- lm(log(y) ~ x)
a <- exp(coef(modelo)[[1]])
b <- coef(modelo)[[2]]
c(a = a, b = b)
##             a             b 
## 10738.9722865    -0.1455165
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 exponencial", xlab = "SiO2 (wt%)", ylab = "MgO (wt%)")
curve(a * exp(b * x), from = min(x), to = max(x), add = TRUE, col = "#d95926", lwd = 3)
legend("topright", c("datos", "y = a·e^(b·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(x, log(y))
r2 <- r^2 * 100
knitr::kable(data.frame(modelo = "exponencial", 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
exponencial -0.879 77.27 fuerte 
residuo <- log(y) - log(a * exp(b * 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 = "SiO2 (wt%)", ylab = "Residuo en logaritmo (ln y real - ln y del modelo)")
abline(h = 0, col = "#d95926", lwd = 2)

7. Restricciones

El modelo se ajustó solo con rocas silicatadas (SiO2 entre 35 y 80 %): por debajo hay carbonatitas, sulfuros o análisis parciales, que no siguen esta tendencia. Da siempre y > 0, pero fuera de ese rango no vale: con SiO2 = 0 daría miles de % de MgO.

c(x_min_datos = min(x), x_max_datos = max(x), y_en_x_100 = a * exp(b * 100))
##  x_min_datos  x_max_datos   y_en_x_100 
## 35.000000000 79.830000000  0.005143531

8. Estimación de pronósticos

Si una roca tiene 60 % de SiO2 (una andesita), ¿cuánto MgO se espera?

x0 <- 60
y0 <- a * exp(b * x0)
y0
## [1] 1.734371

9. Conclusión

Entre SiO2 (wt%) e MgO (wt%) existe una relación que se representa con un modelo exponencial: y = 1.073897{4}·e(-0.1455·x), donde x es SiO2 (wt%) e y es MgO (wt%). La correlación es fuerte (r = -0.879): MgO (wt%) depende de SiO2 (wt%) en un 77.3 % y el 22.7 % restante se debe a otros factores que este modelo no considera. Si SiO2 (wt%) vale 60, se espera MgO (wt%) ≈ 1.73. A más sílice, menos magnesio: es la huella de la diferenciación magmática.