Análisis estadístico de la variable cuantitativa continua: SiO2 (wt%)

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[["SiO2 (wt%)"]]
tipo <- base[["Sample Type"]]                        # tipo de roca de cada fila
tiene <- !is.na(x)                                     # filas con dato
x <- x[tiene]
tipo <- tipo[tiene]
n <- length(x)

3. Límites de outliers (regla intercuartílica)

Con el rango intercuartílico (IR = Q3 − Q1) se fijan los límites Q1 − 1,5·IR y Q3 + 1,5·IR. Los datos que quedan fuera no se eliminan: se cuentan en clases abiertas.

Q1 <- quantile(x, 0.25, names = FALSE)
Q3 <- quantile(x, 0.75, names = FALSE)
IR <- Q3 - Q1                                          # rango intercuartílico
lim_inf <- Q1 - 1.5 * IR                               # menos que esto = outlier bajo
lim_sup <- Q3 + 1.5 * IR                               # más que esto = outlier alto
tipico <- x[x >= lim_inf & x <= lim_sup]               # los intervalos se arman con estos datos
# Los outliers NO se borran: se cuentan en las clases abiertas (son otras rocas, no errores).

4. 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 <- 2                                               # 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

5. 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
36	38	37	478	0.84	1366	56143	2.40	98.44
38	40	39	1110	1.95	2476	55665	4.34	97.60
40	42	41	1471	2.58	3947	54555	6.92	95.66
42	44	43	3035	5.32	6982	53084	12.24	93.08
44	46	45	4680	8.21	11662	50049	20.45	87.76
46	48	47	5034	8.83	16696	45369	29.28	79.55
48	50	49	11797	20.69	28493	40335	49.96	70.72
50	52	51	10632	18.64	39125	28538	68.60	50.04
52	54	53	3196	5.60	42321	17906	74.21	31.40
54	56	55	1641	2.88	43962	14710	77.08	25.79
56	58	57	1357	2.38	45319	13069	79.46	22.92
58	60	59	1146	2.01	46465	11712	81.47	20.54
60	62	61	1255	2.20	47720	10566	83.67	18.53
62	64	63	1157	2.03	48877	9311	85.70	16.33
64	66	65	932	1.63	49809	8154	87.34	14.30

c(abajo = abajo, arriba = arriba)                      # datos en las clases abiertas (no se borra ninguno)

##  abajo arriba 
##    888   7222

6. 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. 1: Histograma y polígono de frecuencias de SiO2 (wt%)")
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 = "SiO2 (wt%)", ylab = "Porcentaje acumulado (%)",
     main = "Gráfica No. 2: 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 = "SiO2 (wt%)",       # diagrama de caja
        main = "Gráfica No. 3: Diagrama de caja (todos los datos)")

7. Error de agrupación: media y desviación con la TDF vs con los datos

Se comparan la media y la desviación calculadas con la TDF (cada dato = la marca de clase de su intervalo) con las calculadas con los datos; menos de 5 % de error indica que la agrupación representa bien los datos.

xc <- x[x >= Li[1] & x <= Ls[k]]                       # los datos de los intervalos normales
media_agrupada <- sum(MC * ni) / sum(ni)               # cada dato = la marca de clase de su intervalo
desv_agrupada <- sqrt(sum((MC - media_agrupada)^2 * ni) / sum(ni))
c(error_media = abs(media_agrupada - mean(xc)) / mean(xc) * 100,   # en %; menos de 5 % = bien
  error_desviacion = abs(desv_agrupada - sd(xc)) / sd(xc) * 100)

##      error_media error_desviacion 
##       0.03220892       0.31294309

8. 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))

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	57031.0000
media	52.2481
mediana	50.0000
desviacion	10.2309
varianza	104.6706
CV	19.5813
asimetria	0.3559
curtosis	3.1224

9. Corrección de la agrupación

Con Sturges sobre todo el rango, los valores extremos estiran los intervalos y casi todos los datos caen en uno o dos. Por eso la TDF de arriba usa el rango típico (sin outliers) y clases abiertas.

9.1. TDF “por defecto” (como la plantilla): Sturges sobre TODO el rango, en la escala original

Para justificar la corrección, se arma la TDF como en la plantilla: Sturges sobre todo el rango, sin tratar los extremos.

x0 <- base[["SiO2 (wt%)"]]
x0 <- x0[!is.na(x0)]
k0 <- round(1 + 3.3 * log10(length(x0)))
a0 <- (max(x0) - min(x0)) / k0                         # amplitud sin redondear
Li0 <- min(x0) + (0:(k0 - 1)) * a0
Ls0 <- Li0 + a0
ni0 <- as.vector(table(cut(x0[x0 <= Ls0[k0]], breaks = c(Li0, Ls0[k0]), right = FALSE, include.lowest = TRUE)))
arriba0 <- sum(x0 > Ls0[k0])                           # por redondeo, el máximo puede quedar afuera
hi0 <- round(c(ni0, arriba0[arriba0 > 0]) / length(x0) * 100, 2)
num <- function(v) trimws(formatC(round(v, 1), format = "fg", digits = 6))   # número a texto: 1 decimal, 6 cifras
etiq0 <- c(paste0("[", num(Li0), "-", num(Ls0), "["), if (arriba0 > 0) paste0("≥ ", num(Ls0[k0])))

9.2. Gráfico: por defecto (izquierda) vs corregida (derecha)

A la izquierda, la agrupación por defecto (casi todo cae en uno o dos intervalos); a la derecha, la corregida.

etiq <- paste0("[", Li, "-", Ls, "[")
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), " %)"))
par(mfrow = c(1, 2), mar = c(8, 5, 4, 1))              # 2 histogramas lado a lado
barplot(hi0, names.arg = etiq0, space = 0, las = 2, cex.names = .65, col = "gray60", border = "white",
        ylab = "Porcentaje (%)", main = paste0("Por defecto: Sturges sobre el rango completo\nSiO2, k = ", length(hi0)))
barplot(hi, names.arg = etiq, space = 0, las = 2, cex.names = .65, col = "#3987e5", border = "white",
        ylab = "Porcentaje (%)", main = paste0("Corregida: rango típico (IQR) + clases abiertas\nSiO2 (wt%), k = ", k, " cerradas"))
mtext(if (length(nota) > 0) paste0("Clases abiertas (fuera del gráfico): ", paste(nota, collapse = " · ")) else "",
      side = 1, line = 6.5, cex = .65, col = "gray30")

10. Estrato: rocas volcánicas máficas

El mismo análisis para un solo tipo de roca (basaltos y afines). Al separar por estrato, la distribución deja de mezclar rocas distintas.

10.1. Datos