Proyecto estrictamente académico. 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 y librerías

Se lee la base limpia (159 760 análisis de roca total × 75 variables). La primera vez se descarga del repositorio del proyecto (26 MB).

knitr::opts_chunk$set(
  echo = TRUE,          # muestra el código en el documento
  message = FALSE,      # sin mensajes de carga de paquetes
  warning = FALSE,      # sin advertencias
  fig.align = "center" # gráficos centrados
)
library(e1071)         # asimetría y curtosis (install.packages("e1071"))
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)              # filas y columnas
## [1] 159760     75

2. Datos

Se toma la columna de la variable y se quitan las filas vacías (vacío = no medido, nunca cero).

x <- base[["N Analytes"]]
x <- x[!is.na(x)]                                      # quita vacíos
n <- length(x)
tipico <- x                                            # aquí se agrupan todos los datos

3. Intervalos (regla de Sturges)

El número de intervalos sale de la regla de Sturges, k = 1 + 3,3·log10(n); la amplitud se redondea a un número cómodo para que los límites sean fáciles de leer.

k <- round(1 + 3.3 * log10(n))                         # número de intervalos
a <- (max(tipico) - min(tipico)) / k                   # amplitud calculada
a <- 5                                               # se redondea hacia arriba a un número cómodo
desde <- floor(min(tipico) / a) * a                    # primer límite: múltiplo de a
k <- floor((max(tipico) - desde) / a) + 1              # intervalos que hacen falta para llegar al máximo
Li <- desde + (0:(k - 1)) * a                          # límites inferiores
Ls <- Li + a                                           # límites superiores
MC <- (Li + Ls) / 2                                    # marca de clase

4. Frecuencias

Tabla de distribución de frecuencias: límites, marca de clase, frecuencia absoluta y relativa, y acumuladas ascendentes y descendentes.

dentro <- x[x >= Li[1] & x <= Ls[k]]
ni <- as.vector(table(cut(dentro, breaks = c(Li, Ls[k]), right = FALSE, include.lowest = TRUE)))  # cuenta por intervalo [Li, Ls[
abajo <- sum(x < Li[1])                                # clase abierta "menor que" (outliers bajos)
arriba <- sum(x > Ls[k])                               # clase abierta "mayor o igual" (outliers altos)
hi <- round(ni / n * 100, 2)                           # frecuencia relativa (%)
todos <- c(abajo, ni, arriba)                          # con las clases abiertas, para acumular
cerradas <- 2:(k + 1)                                  # posiciones de los intervalos normales
Ni_asc <- cumsum(todos)[cerradas]                      # acumulada ascendente
Ni_dsc <- rev(cumsum(rev(todos)))[cerradas]            # acumulada descendente
Hi_asc <- round(cumsum(todos) / n * 100, 2)[cerradas]
Hi_dsc <- round(rev(cumsum(rev(todos))) / n * 100, 2)[cerradas]
TDF <- data.frame(Li, Ls, MC, ni, hi, Ni_asc, Ni_dsc, Hi_asc, Hi_dsc)
knitr::kable(TDF, caption = "Tabla N° 1: distribución de frecuencias")
Tabla N° 1: distribución de frecuencias
Li Ls MC ni hi Ni_asc Ni_dsc Hi_asc Hi_dsc
0 5 2.5 62170 38.91 62170 159760 38.91 100.00
5 10 7.5 33020 20.67 95190 97590 59.58 61.09
10 15 12.5 23092 14.45 118282 64570 74.04 40.42
15 20 17.5 9480 5.93 127762 41478 79.97 25.96
20 25 22.5 10366 6.49 138128 31998 86.46 20.03
25 30 27.5 5781 3.62 143909 21632 90.08 13.54
30 35 32.5 5607 3.51 149516 15851 93.59 9.92
35 40 37.5 3822 2.39 153338 10244 95.98 6.41
40 45 42.5 3382 2.12 156720 6422 98.10 4.02
45 50 47.5 1547 0.97 158267 3040 99.07 1.90
50 55 52.5 751 0.47 159018 1493 99.54 0.93
55 60 57.5 409 0.26 159427 742 99.79 0.46
60 65 62.5 131 0.08 159558 333 99.87 0.21
65 70 67.5 116 0.07 159674 202 99.95 0.13
70 75 72.5 56 0.04 159730 86 99.98 0.05
75 80 77.5 30 0.02 159760 30 100.00 0.02 
c(abajo = abajo, arriba = arriba)                      # datos en las clases abiertas (no se borra ninguno)
##  abajo arriba 
##      0      0

5. Gráficos: histograma + polígono, ojiva y caja

Histograma con polígono de frecuencias, ojivas ascendente y descendente (se cruzan en la mediana) y diagrama de caja.

etiq <- paste0("[", Li, "-", Ls, "[")                  # texto de cada intervalo
nota <- c(if (abajo > 0) paste0("< ", Li[1], ": ", abajo, " (", round(abajo / n * 100, 2), " %)"),
          if (arriba > 0) paste0("≥ ", Ls[k], ": ", arriba, " (", round(arriba / n * 100, 2), " %)"))
layout(matrix(c(1, 1, 2, 3), 2, 2, byrow = TRUE), heights = c(1.15, 1))   # arriba 1 gráfico, abajo 2
par(mar = c(8, 5, 4, 1))
h <- barplot(hi, names.arg = etiq, space = 0, las = 2, col = "#3987e5", border = "white",   # histograma
             ylab = "Porcentaje (%)", cex.names = .75, ylim = c(0, max(hi) * 1.15),
             main = "Gráfica No. 4: Histograma y polígono de frecuencias de Número de analitos (analitos)")
lines(c(h[1] - 1, h, h[length(h)] + 1), c(0, hi, 0), col = "#d95926", lwd = 2, type = "b", pch = 16, cex = .6)  # polígono
mtext(if (length(nota) > 0) paste0("Clases abiertas (fuera del gráfico): ", paste(nota, collapse = " · ")) else "",
      side = 1, line = 6.5, cex = .7, col = "gray30")  # nota bajo el histograma
par(mar = c(5, 5, 4, 1))
plot(Ls, Hi_asc, type = "b", pch = 16, col = "#3987e5", ylim = c(0, 100),            # ojiva ascendente
     xlab = "Número de analitos (analitos)", ylab = "Porcentaje acumulado (%)",
     main = "Gráfica No. 5: Ojiva ascendente y descendente")
lines(Li, Hi_dsc, type = "b", pch = 17, col = "#d95926")                              # ojiva descendente
abline(h = 50, lty = 3, col = "gray50")                # línea del 50 %: donde se cruzan está la mediana
legend("right", c("ascendente (≤)", "descendente (≥)"), col = c("#3987e5", "#d95926"), pch = c(16, 17), bty = "n", cex = .85)
boxplot(x, horizontal = TRUE, col = "#86b6ef", xlab = "Número de analitos (analitos)",       # diagrama de caja
        main = "Gráfica No. 6: Diagrama de caja (todos los datos)")

6. Indicadores (con los datos sin agrupar)

Medidas de tendencia central, dispersión y forma con los datos sin agrupar (asimetría y curtosis con el paquete e1071).

ind <- c(n = n, media = mean(x), mediana = median(x), desviacion = sd(x), varianza = var(x),
  CV = sd(x) / mean(x) * 100, asimetria = skewness(x), curtosis = kurtosis(x))
as.numeric(names(which.max(table(x))))                # moda: el valor que más se repite
## [1] 1
knitr::kable(data.frame(indicador = names(ind), valor = round(ind, 4)), row.names = FALSE,
             caption = "Tabla N° 2: indicadores estadísticos")
Tabla N° 2: indicadores estadísticos
indicador valor
n 159760.0000
media 11.2085
mediana 9.0000
desviacion 11.7317
varianza 137.6325
CV 104.6674
asimetria 1.5847
curtosis 2.3482

7. Conclusión

Se analizaron 159 760 casos con dato. La mediana es 9 analitos; la asimetría es 1,58 (cola hacia los valores altos) y el coeficiente de variación es 104,67 %. La moda es 1. Muchas muestras tienen pocos analitos: no todas tienen un análisis completo. Por eso, al comparar dos elementos se usan solo las filas que tienen ambos datos, sin inventar los que faltan.