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 = La (ppm) (variable independiente) e y = Ce (ppm) (variable dependiente). Son tierras raras vecinas en la tabla periódica, con radio iónico y carga casi iguales: entran y salen juntas de los minerales, así que su proporción casi no cambia (Ce ≈ 2·La). Unos pocos análisis tienen La cercano a 1 000 ppm con apenas 26 ppm de Ce, una razón imposible en rocas naturales (errores de carga o anomalías de Ce extremas): tienen una distancia de Cook enorme y tiran la recta hacia abajo. Se ajusta con las rocas de razón Ce/La entre 0,5 y 5.

x <- base[["La (ppm)"]]
y <- base[["Ce (ppm)"]]
completos <- !is.na(x) & !is.na(y) & x > 0 & y > 0    # filas con los dos datos
x <- x[completos]
y <- y[completos]
# Influencia: la distancia de Cook mide cuánto cambia la recta si se quita cada punto
todos <- lm(y ~ x)
cook <- cooks.distance(todos)
head(data.frame(La = x, Ce = y, cook = round(cook, 2))[order(-cook), ], 5)   # los 5 puntos más influyentes
##           La    Ce  cook
## 17115 1441.0 370.0 31.41
## 15     965.0  26.0  8.49
## 3320   964.9  68.9  7.99
## 2558  1006.0 223.0  7.62
## 12857  917.4  64.4  6.52
razon <- y / x                                         # Ce/La: en rocas naturales está entre 0,5 y 5
natural <- razon >= 0.5 & razon <= 5
c(pares = length(x), fuera_de_0.5_a_5 = sum(!natural), porcentaje = round(mean(!natural) * 100, 1))
##            pares fuera_de_0.5_a_5       porcentaje 
##          31690.0            769.0              2.4
c(pendiente_con_todos = coef(todos)[[2]], pendiente_sin_esos = coef(lm(y[natural] ~ x[natural]))[[2]])
## pendiente_con_todos  pendiente_sin_esos 
##            1.525529            1.674608
x <- x[natural]                                        # el modelo se ajusta sin ellos (no se borran de la base)
y <- y[natural]
length(x)                                              # número de pares
## [1] 30921

3. Tabla de pares de valores

knitr::kable(head(data.frame(x, y), 10), caption = "Tabla N° 1: primeros 10 pares de valores (La (ppm), Ce (ppm))")
Tabla N° 1: primeros 10 pares de valores (La (ppm), Ce (ppm))
x y
9.000 23.900
0.267 0.983
7.000 11.000
8.750 20.550
4.160 8.970
5.340 13.790
6.650 17.500
2.640 6.610
0.660 2.300
30.000 58.000 
knitr::kable(data.frame(variable = c("La (ppm)", "Ce (ppm)"), 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
La (ppm) 0 22.23 2337
Ce (ppm) 0 44.10 3975

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, xlim = c(0, quantile(x, 0.995)), ylim = c(0, quantile(y, 0.995)),
     main = "Gráfica N° 1: Diagrama de dispersión", xlab = "La (ppm)", ylab = "Ce (ppm)")

Los puntos forman una franja recta que sale del origen: a más La, más Ce.

5. Conjetura del modelo

La nube es recta: se conjetura un modelo lineal, y = 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.970 94.1
potencial 0.990 98.0
exponencial 0.493 24.3
logaritmico 0.555 30.8
polinomico 0.977 95.4 
modelo <- lm(y ~ x)
a <- coef(modelo)[[1]]                              # intercepto
b <- coef(modelo)[[2]]                              # pendiente
c(a = a, b = b)
##        a        b 
## 6.869752 1.674608
plot(x, y, pch = 16, cex = 0.4, col = rgb(0.22, 0.53, 0.9, 0.25), las = 1, xlim = c(0, quantile(x, 0.995)), ylim = c(0, quantile(y, 0.995)),
     main = "Gráfica N° 2: Regresión lineal", xlab = "La (ppm)", ylab = "Ce (ppm)")
curve(a + b * x, from = min(x), to = max(x), add = TRUE, col = "#d95926", lwd = 3)
legend("topleft", c("datos", "y = a + b·x"), col = c(rgb(0.22, 0.53, 0.9), "#d95926"),
       pch = c(16, NA), lwd = c(NA, 3), bty = "n")

Mínimos cuadrados a mano. La recta que hace mínima la suma de los residuos al cuadrado tiene pendiente b = Sxy / Sxx e intercepto a = ȳ − b·x̄. Da lo mismo que lm():

Sxy <- sum((x - mean(x)) * (y - mean(y)))              # suma de productos cruzados
Sxx <- sum((x - mean(x))^2)                            # suma de cuadrados de x
b_manual <- Sxy / Sxx
a_manual <- mean(y) - b_manual * mean(x)
covarianza <- Sxy / (length(x) - 1)
c(a_manual = a_manual, b_manual = b_manual, a_lm = a, b_lm = b, covarianza = covarianza)
##    a_manual    b_manual        a_lm        b_lm  covarianza 
##    6.869752    1.674608    6.869752    1.674608 4189.738827

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, y)
r2 <- r^2 * 100
knitr::kable(data.frame(modelo = "lineal", 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
lineal 0.9699 94.07 fuerte 
residuo <- y - (a + 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, xlim = c(0, quantile(x, 0.995)),
     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 = "La (ppm)", ylab = "Residuo (y real - y del modelo)")
abline(h = 0, col = "#d95926", lwd = 2)

7. Restricciones

Se ajusta con rocas de razón Ce/La entre 0,5 y 5 (el 97 % de los pares). Una concentración no puede ser negativa: y > 0 exige x > −a/b. El modelo solo vale dentro del rango medido de La; fuera de él es extrapolar.

c(x_minimo_valido = max(0, -a / b), x_min_datos = min(x), x_max_datos = max(x))
## x_minimo_valido     x_min_datos     x_max_datos 
##          0.0000          0.0001       2337.0000

8. Estimación de pronósticos

Si un basalto tiene 20 ppm de La, ¿cuánto Ce se espera?

x0 <- 20
y0 <- a + b * x0
y0
## [1] 40.36192

9. Conclusión

Entre La (ppm) e Ce (ppm) existe una relación que se representa con un modelo lineal: y = 6.87 + 1.675·x, donde x es La (ppm) e y es Ce (ppm). La correlación es fuerte (r = 0.97): Ce (ppm) depende de La (ppm) en un 94.1 % y el 5.9 % restante se debe a otros factores que este modelo no considera. Si La (ppm) vale 20, se espera Ce (ppm) ≈ 40.36. El Ce es casi el doble del La en las rocas, como se espera de dos tierras raras livianas.