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

y = Nd (ppm), explicada por x1 = La (ppm) y x2 = Sm (ppm): el neodimio está entre el lantano y el samario en la serie de las tierras raras, así que debería poder estimarse con sus dos vecinos. Además se agrega una variable cualitativa, el ambiente tectónico, como variable ficticia (dummy).

datos <- data.frame(Nd = base[["Nd (ppm)"]], La = base[["La (ppm)"]], Sm = base[["Sm (ppm)"]],
                    grupo = base[["Tectonic Group"]])
datos <- na.omit(datos)                                # filas con todo completo
datos <- datos[datos$Nd > 0 & datos$La > 0 & datos$Sm > 0 & datos$grupo != "mixto", ]
nrow(datos)
## [1] 26765

3. Tabla de valores

knitr::kable(head(datos, 10), caption = "Tabla N° 1: primeras 10 filas")
Tabla N° 1: primeras 10 filas
Nd La Sm grupo
2 0.470 0.180 0.190 dorsal y corteza oceanica
7 14.300 9.000 3.980 margen convergente
9 1.544 0.267 0.789 dorsal y corteza oceanica
19 6.750 1.480 2.340 margen convergente
21 12.770 8.750 3.810 intraplaca y rift
27 6.260 4.160 0.890 dorsal y corteza oceanica
35 10.190 5.340 3.100 margen convergente
38 11.400 6.650 3.410 intraplaca y rift
50 4.890 2.640 1.590 intraplaca y rift
62 15.820 10.790 3.930 margen convergente

4. Diagrama de dispersión y matriz de correlación

set.seed(2026)
dibujo <- datos[sample(nrow(datos), 3000), c("Nd", "La", "Sm")]   # 3000 puntos para que se vea
dibujo <- dibujo[dibujo$Nd < quantile(datos$Nd, 0.995), ]          # sin el 0,5 % más extremo
pairs(dibujo, pch = 16, cex = 0.4, col = rgb(0.22, 0.53, 0.9, 0.3), main = "Gráfica N° 1: Dispersión por pares")

R <- cor(datos[, c("Nd", "La", "Sm")])
knitr::kable(round(R, 3), caption = "Tabla N° 2: matriz de correlación")
Tabla N° 2: matriz de correlación
Nd La Sm
Nd 1.000 0.797 0.765
La 0.797 1.000 0.707
Sm 0.765 0.707 1.000

Las dos explicativas están muy correlacionadas entre sí (La–Sm r = 0.71): es multicolinealidad. No impide predecir, pero hace que los coeficientes por separado sean menos estables.

r2_La_Sm <- summary(lm(La ~ Sm, data = datos))$r.squared
VIF <- 1 / (1 - r2_La_Sm)                              # factor de inflación de la varianza
VIF                                                    # más de 10 = multicolinealidad fuerte
## [1] 1.999439

5. Conjetura del modelo

Modelo lineal múltiple: y = a + b1·x1 + b2·x2. En forma matricial, β = (XᵀX)⁻¹·Xᵀy; lm() hace ese cálculo.

modelo <- lm(Nd ~ La + Sm, data = datos)
coef(modelo)
## (Intercept)          La          Sm 
##   1.7099733   0.3680367   2.4465774
# El mismo cálculo con matrices:
X <- cbind(1, datos$La, datos$Sm)
beta <- solve(t(X) %*% X) %*% t(X) %*% datos$Nd
round(t(beta), 4)
##      [,1]  [,2]   [,3]
## [1,] 1.71 0.368 2.4466
predicho <- predict(modelo)
plot(predicho, datos$Nd, pch = 16, cex = 0.4, col = rgb(0.22, 0.53, 0.9, 0.25), las = 1,
     xlim = c(0, quantile(datos$Nd, 0.995)), ylim = c(0, quantile(datos$Nd, 0.995)),   # zoom al 99,5 %
     main = "Gráfica N° 2: Nd real vs Nd del modelo", xlab = "Nd del modelo (ppm)", ylab = "Nd real (ppm)")
abline(a = 0, b = 1, col = "#d95926", lwd = 3)

Con variable ficticia. El ambiente tectónico entra como dos columnas 0/1 (la tercera categoría es la referencia).

modelo_dummy <- lm(Nd ~ La + Sm + grupo, data = datos)
round(coef(modelo_dummy), 4)
##                    (Intercept)                             La 
##                         0.8153                         0.3664 
##                             Sm grupodorsal y corteza oceanica 
##                         2.4358                        -0.0883 
##         grupointraplaca y rift        grupomargen convergente 
##                         1.4272                         1.6127
c(R2_sin_dummy = summary(modelo)$r.squared, R2_con_dummy = summary(modelo_dummy)$r.squared)
## R2_sin_dummy R2_con_dummy 
##    0.7166281    0.7169916

6. Test de aprobación

R2 <- summary(modelo)$r.squared * 100
r <- sqrt(R2 / 100)                                    # correlación múltiple
knitr::kable(data.frame(modelo = "múltiple", R = round(r, 4), R2_porcentaje = round(R2, 2)),
             caption = "Tabla N° 3: test de aprobación")
Tabla N° 3: test de aprobación
modelo R R2_porcentaje
múltiple 0.8465 71.66 
plot(predicho, residuals(modelo), pch = 16, cex = 0.4, col = rgb(0.22, 0.53, 0.9, 0.25), las = 1,
     xlim = c(0, quantile(datos$Nd, 0.995)), ylim = quantile(residuals(modelo), c(0.0025, 0.9975)),
     main = "Gráfica N° 3: Residuos", xlab = "Nd del modelo (ppm)", ylab = "Residuo (ppm)")
abline(h = 0, col = "#d95926", lwd = 2)

7. Restricciones

Las concentraciones son positivas: el modelo vale dentro del rango medido de La y Sm, y para combinaciones que existen en la naturaleza (La y Sm suben juntos; un La muy alto con Sm muy bajo no se da).

sapply(datos[, c("La", "Sm")], range)
##             La       Sm
## [1,]    0.0001 8.01e-05
## [2,] 2337.0000 2.90e+02

8. Estimación de pronósticos

Si una roca tiene 20 ppm de La y 5 ppm de Sm, ¿cuánto Nd se espera?

nuevo <- data.frame(La = 20, Sm = 5)
y0 <- predict(modelo, nuevo)
y0
##        1 
## 21.30359

9. Conclusión

El Nd se representa con un modelo lineal múltiple: y = 1.71 + 0.368·x1 + 2.447·x2, donde x1 es La y x2 es Sm. El modelo explica el 71.7 % de la variación del Nd. Para La = 20 ppm y Sm = 5 ppm se espera Nd ≈ 21.3 ppm. Agregar el ambiente tectónico como variable ficticia casi no mejora el ajuste: las tierras raras se mueven juntas en cualquier ambiente.