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 = CaO (wt%) (variable dependiente). El calcio sube con el magnesio en los basaltos (plagioclasa y clinopiroxeno) y vuelve a bajar en las peridotitas, donde domina el olivino, que casi no lleva calcio.

x <- base[["MgO (wt%)"]]
y <- base[["CaO (wt%)"]]
completos <- !is.na(x) & !is.na(y) & x > 0 & y > 0    # filas con los dos datos
x <- x[completos]
y <- y[completos]
length(x)                                              # número de pares
## [1] 54414

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%), CaO (wt%))")
Tabla N° 1: primeros 10 pares de valores (MgO (wt%), CaO (wt%))
x y
7.88 11.65
8.77 11.07
1.65 4.16
5.32 11.87
7.59 11.18
7.57 12.30
6.34 9.53
44.02 0.54
6.75 11.63
4.20 12.21 
knitr::kable(data.frame(variable = c("MgO (wt%)", "CaO (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.15 54.76
CaO (wt%) 0 7.68 97.29

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 = "CaO (wt%)")

El CaO sube hasta un máximo cerca de 14 % de MgO y después baja: una curva con joroba. Una recta no puede seguirla (r lineal = −0,28).

5. Conjetura del modelo

Una curva que sube y baja es la forma de un modelo polinómico, aquí de grado 3: y = a + b·x + c·x² + d·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.283 8.0
potencial 0.302 9.1
exponencial -0.329 10.8
logaritmico 0.275 7.6
polinomico 0.799 63.8 
modelo <- lm(y ~ x + I(x^2) + I(x^3))
a <- coef(modelo)[[1]]
b <- coef(modelo)[[2]]
c <- coef(modelo)[[3]]
d <- coef(modelo)[[4]]
c(a = a, b = b, c = c, d = d)
##            a            b            c            d 
##  2.042879364  1.770996854 -0.085452332  0.001018958
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 polinómico", xlab = "MgO (wt%)", ylab = "CaO (wt%)")
curve(a + b * x + c * x^2 + d * x^3, from = min(x), to = max(x), add = TRUE, col = "#d95926", lwd = 3)
legend("topright", c("datos", "y = a + b·x + c·x² + d·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 <- sqrt(summary(modelo)$r.squared)
r2 <- r^2 * 100
knitr::kable(data.frame(modelo = "polinómico", 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
polinómico 0.7987 63.79 fuerte 
residuo <- y - (a + b * x + c * x^2 + d * x^3)                                   # 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

Un polinomio se dispara fuera del rango de los datos (riesgo de sobreajuste si se sube el grado). Además el CaO debe ser positivo: se busca dónde la curva baja de cero.

xs <- seq(min(x), max(x), length.out = 2000)
ys <- a + b * xs + c * xs^2 + d * xs^3
c(x_min_datos = min(x), x_max_datos = max(x),
  primer_x_con_y_negativo = if (any(ys < 0)) xs[which(ys < 0)[1]] else NA,
  x_del_maximo = xs[which.max(ys)])
##             x_min_datos             x_max_datos primer_x_con_y_negativo 
##                0.000349               54.760000                      NA 
##            x_del_maximo 
##               13.751897

8. Estimación de pronósticos

Si un basalto tiene 8 % de MgO, ¿cuánto CaO se espera?

x0 <- 8
y0 <- a + b * x0 + c * x0^2 + d * x0^3
y0
## [1] 11.26361

9. Conclusión

Entre MgO (wt%) e CaO (wt%) existe una relación que se representa con un modelo polinómico: y = 2.043 + 1.771·x + (-0.0855)·x² + 0.00102·x³, donde x es MgO (wt%) e y es CaO (wt%). La correlación es fuerte (r = 0.799): CaO (wt%) depende de MgO (wt%) en un 63.8 % y el 36.2 % restante se debe a otros factores que este modelo no considera. Si MgO (wt%) vale 8, se espera CaO (wt%) ≈ 11.26. Una sola curva describe dos tipos de roca distintos: basaltos a la izquierda del máximo y peridotitas a la derecha.