Análisis estadístico de la variable cuantitativa continua: MgO (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[["MgO (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
0	2	1	10312	18.08	10312	57039	18.08	100.00
2	4	3	4566	8.01	14878	46727	26.08	81.92
4	6	5	7139	12.52	22017	42161	38.60	73.92
6	8	7	13705	24.03	35722	35022	62.63	61.40
8	10	9	7075	12.40	42797	21317	75.03	37.37
10	12	11	2012	3.53	44809	14242	78.56	24.97
12	14	13	1082	1.90	45891	12230	80.46	21.44
14	16	15	758	1.33	46649	11148	81.78	19.54
16	18	17	491	0.86	47140	10390	82.65	18.22
18	20	19	437	0.77	47577	9899	83.41	17.35

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

##  abajo arriba 
##      0   9462

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. 7: Histograma y polígono de frecuencias de MgO (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 = "MgO (wt%)", ylab = "Porcentaje acumulado (%)",
     main = "Gráfica No. 8: 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 = "MgO (wt%)",       # diagrama de caja
        main = "Gráfica No. 9: 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 
##         1.903710         3.049593

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	57039.0000
media	11.3087
mediana	7.0500
desviacion	13.0648
varianza	170.6891
CV	115.5286
asimetria	1.6670
curtosis	1.4113

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[["MgO (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\nMgO, 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\nMgO (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