ANÁLISIS ESTADÍSTICO

CARGA DE DATOS Y LIBRERÍAS

CARGA DE DATOS

#Limpiar entorno
rm(list = ls())
#Cargar librerías
if (!require("readr")) install.packages("readr")
if (!require("dplyr")) install.packages("dplyr")
if (!require("knitr")) install.packages("knitr")
if (!require("moments")) install.packages("moments")

library(readr)
library(dplyr)
library(knitr)
library(moments)

#Cargar datos
ruta <- "D:/SIO2_filled.csv"
datos <- read_csv(ruta)

cat("✓ Datos cargados exitosamente\n")
## ✓ Datos cargados exitosamente
cat("✓ Dimensiones:", nrow(datos), "observaciones\n")
## ✓ Dimensiones: 2500 observaciones

EXTRAER Y LIMPIAR LA VARIABLE

ASIGNACION DE VARIABLES

#Extraer y limpiar variable de Sílice
SIO2_raw <- datos$SIO2_WT_PE
SIO2 <- as.numeric(SIO2_raw)
SIO2 <- SIO2[!is.na(SIO2)]
SIO2 <- SIO2[SIO2 != 0]

n <- length(SIO2)

TABLA DE DISTRIBUCION DE PARAMETROS POR STURGES

# Parámetros de Sturges
R <- max(SIO2) - min(SIO2)
k <- floor(1 + 3.3 * log10(n))
A <- R / k

parametros_sturges <- data.frame(
  Parámetro = c("Rango (R)", 
                "Número de datos (n)", 
                "Número de intervalos (k)", 
                "Amplitud de clase (A)"),
  Valor = c(round(R, 4), n, k, round(A, 4))
)

kable(parametros_sturges, 
      caption = "Tabla 1: Parámetros de Sturges para la variable SIO2"
)
Tabla 1: Parámetros de Sturges para la variable SIO2
Parámetro Valor
Rango (R) 37.5400
Número de datos (n) 2500.0000
Número de intervalos (k) 12.0000
Amplitud de clase (A) 3.1283 
# =====================================================
# TABLA COMPLETA DE STURGES 
# =====================================================

# Asegurar que n no tenga NAs
SIO2_clean <- na.omit(SIO2)
n <- length(SIO2_clean)

# Generación precisa de límites
Li <- seq(min(SIO2_clean), max(SIO2_clean) - A, length.out = k)
Ls <- seq(min(SIO2_clean) + A, max(SIO2_clean), length.out = k)

# Corrección de precisión: forzar que el último límite sea el máximo exacto
Ls[k] <- max(SIO2_clean)

MC <- (Li + Ls) / 2

ni <- numeric(k)

for (i in 1:k) {
  if (i == k) {
    ni[i] <- sum(SIO2_clean >= Li[i] & SIO2_clean <= Ls[i])
  } else {
    ni[i] <- sum(SIO2_clean >= Li[i] & SIO2_clean < Ls[i])
  }
}

# Recalcular n total real contado para evitar descalce
n_real <- sum(ni)

hi <- round((ni / n_real) * 100, 2)

Niasc <- cumsum(ni)
Nidsc <- rev(cumsum(rev(ni)))

# Acumulados calculados y forzados al 100 en el límite correspondiente
Hiasc <- round((Niasc / n_real) * 100, 2)
Hidsc <- round((Nidsc / n_real) * 100, 2)

# Forzado explícito de extremos para corregir redondeos visuales
Hiasc[k] <- 100.00
Hidsc[1] <- 100.00

TDF <- data.frame(
  Clase = 1:k,
  Li = round(Li, 2),
  Ls = round(Ls, 2),
  MC = round(MC, 2),
  ni = ni,
  hi = hi,
  Niasc = Niasc,
  Nidsc = Nidsc,
  Hiasc = Hiasc,
  Hidsc = Hidsc
)

kable(
  TDF,
  caption = "Tabla 2: Distribución de frecuencias"
)
Tabla 2: Distribución de frecuencias
Clase Li Ls MC ni hi Niasc Nidsc Hiasc Hidsc
1 42.46 45.59 44.02 5 0.20 5 2499 0.20 100.00
2 45.59 48.72 47.15 21 0.84 26 2494 1.04 99.80
3 48.72 51.84 50.28 50 2.00 76 2473 3.04 98.96
4 51.84 54.97 53.41 190 7.60 266 2423 10.64 96.96
5 54.97 58.10 56.54 427 17.09 693 2233 27.73 89.36
6 58.10 61.23 59.67 389 15.57 1082 1806 43.30 72.27
7 61.23 64.36 62.79 259 10.36 1341 1417 53.66 56.70
8 64.36 67.49 65.92 218 8.72 1559 1158 62.38 46.34
9 67.49 70.61 69.05 302 12.08 1861 940 74.47 37.62
10 70.62 73.74 72.18 391 15.65 2252 638 90.12 25.53
11 73.74 76.87 75.31 207 8.28 2459 247 98.40 9.88
12 76.87 80.00 78.44 40 1.60 2499 40 100.00 1.60

Debido a que los límites de clase obtenidos mediante la regla de Sturges presentan números decimales difíciles de interpretar visualmente se simplifico la tabla.

# =====================================================
# TABLA RESUMEN BASADA EN HISTOGRAMA 
# =====================================================

# Se definen los límites explícitos para asegurar exactamente k intervalos (8)
cortes_k <- seq(min(SIO2), max(SIO2), length.out = k + 1)

h <- hist(
  SIO2,
  breaks = cortes_k,
  plot = FALSE
)

Li <- head(h$breaks, -1)
Ls <- tail(h$breaks, -1)

ni <- h$counts
n <- sum(ni)

hi <- round((ni / n) * 100, 2)

Niasc <- cumsum(ni)
Nidsc <- rev(cumsum(rev(ni)))

Hiasc <- round(cumsum(hi), 2)
Hidsc <- round(rev(cumsum(rev(hi))), 2)

MC <- h$mids

TDF <- data.frame(
  Clase = seq_along(ni),
  Li = round(Li, 0),
  Ls = round(Ls, 0),
  MC = round(MC, 0),
  ni = ni,
  hi = hi,
  Niasc = Niasc,
  Nidsc = Nidsc,
  Hiasc = Hiasc,
  Hidsc = Hidsc
)

kable(
  TDF,
  caption = "Tabla 3: Distribución de frecuencias resumen con acumulados exactos"
)
Tabla 3: Distribución de frecuencias resumen con acumulados exactos
Clase Li Ls MC ni hi Niasc Nidsc Hiasc Hidsc
1 42 46 44 5 0.20 5 2500 0.20 100.00
2 46 49 47 21 0.84 26 2495 1.04 99.80
3 49 52 50 50 2.00 76 2474 3.04 98.96
4 52 55 53 190 7.60 266 2424 10.64 96.96
5 55 58 57 427 17.08 693 2234 27.72 89.36
6 58 61 60 390 15.60 1083 1807 43.32 72.28
7 61 64 63 259 10.36 1342 1417 53.68 56.68
8 64 67 66 218 8.72 1560 1158 62.40 46.32
9 67 71 69 302 12.08 1862 940 74.48 37.60
10 71 74 72 391 15.64 2253 638 90.12 25.52
11 74 77 75 207 8.28 2460 247 98.40 9.88
12 77 80 78 40 1.60 2500 40 100.00 1.60

GRAFICAS DE DISTRIBUCION DE CANTIDAD

# Histograma inicial
hist(
  SIO2,
  breaks = h$breaks,
  col = "gray",
  border = "black",
  main = "Gráfica 1: Distribución de cantidad de SIO2\nen depósitos minerales de Estados Unidos",
  xlab = "Sílice (%)",
  ylab = "Frecuencia"
)

hist(
  SIO2,
  breaks = k,
  col = "gray",
  main = "Gráfica 2: Distribución global de SIO2\nen depósitos minerales de Estados Unidos",
  xlab = "Sílice (%)",
  ylab = "Frecuencia"
)

# Gráfica 3: Distribución porcentual
h <- hist(SIO2,
          breaks = k,
          plot = FALSE)

porcentaje <- h$counts / sum(h$counts) * 100

par(mar = c(7, 4, 4, 2))

barplot(
  porcentaje,
  names.arg = paste(round(head(h$breaks, -1), 0),
                    round(tail(h$breaks, -1), 0),
                    sep = "-"),
  col = "darkorange",
  ylim = c(0, max(porcentaje) * 1.15),
  ylab = "Frecuencia Relativa (%)",
  main = "Gráfica 3: Distribución porcentual de SIO2",
  las = 2,
  space = 0,
  border = "white"
)

title(
  xlab = "Sílice (%)",
  line = 5.5
)

# Gráfica 4: Frecuencia Relativa
media_rel_global <- mean(SIO2, na.rm = TRUE)
mediana_rel_global <- median(SIO2, na.rm = TRUE)

h <- hist(SIO2,
          breaks = k,
          plot = FALSE)

hi <- h$counts / sum(h$counts) * 100

intervalos_graf <- paste(
  round(head(h$breaks, -1), 0),
  round(tail(h$breaks, -1), 0),
  sep = " - "
)

par(mar = c(7, 4, 4, 2))

barplot(
  hi,
  names.arg = intervalos_graf,
  col = "mediumpurple",
  ylim = c(0, 100),
  cex.names = 0.7,
  space = 0,
  ylab = "Frecuencia Relativa (%)",
  xlab = "",
  main = "Gráfica 4: Distribución de cantidad en porcentaje de SIO2\nen depósitos minerales de Estados Unidos",
  las = 2,
  border = "white"
)

title(
  xlab = "Intervalos de Sílice (%)",
  line = 6
)

breaks_sio2 <- seq(
  min(SIO2, na.rm = TRUE),
  max(SIO2, na.rm = TRUE) + A,
  by = A
)

histograma <- hist(
  SIO2,
  breaks = breaks_sio2,
  col = rgb(0.7, 0.7, 0.7, 0.7),
  border = "black",
  main = "Gráfica Nº5: Histograma con polígono de frecuencias\nde SIO2 en depósitos minerales de Estados Unidos",
  xlab = "Sílice (%)",
  ylab = "Cantidad"
)

lines(
  histograma$mids,
  histograma$counts,
  type = "o",
  pch = 16,
  lwd = 2
)

# Gráfica 6: Ojiva Absoluta
plot(Ls, Niasc, 
     type = "o", 
     pch = 16, 
     col = "blue",
     main = "Gráfica 6: Ojiva de Frecuencias Absolutas Acumuladas \nde SIO2 de la muestra",
     xlab = "Sílice (%)",
     ylab = "Frecuencia Acumulada Absoluta (Ni)")

lines(Ls, Nidsc, type = "o", pch = 16, col = "red")

legend("topleft",
       c("Ni Ascendente", "Ni Descendente"),
       col = c("blue", "red"),
       pch = 16)

# Gráfica 7: Ojiva Relativa
plot(Ls, Hiasc, 
     type = "o", 
     pch = 16, 
     col = "blue",
     main = "Gráfica 7: Ojiva de Frecuencias Relativas Acumuladas \nde SIO2 de la muestra",
     xlab = "Sílice (%)",
     ylab = "Frecuencia Acumulada Relativa (%)")

lines(Ls, Hidsc, type = "o", pch = 16, col = "red")

legend("bottomright", 
       c("Hi Ascendente", "Hi Descendente"),
       col = c("blue", "red"), 
       pch = 16)

# Gráfica 8 y 9: Diagrama de Caja y Superposición
# Estadísticos de SIO2
media_box <- mean(SIO2, na.rm = TRUE)
mediana_box <- median(SIO2, na.rm = TRUE)
Q1_box <- quantile(SIO2, 0.25, na.rm = TRUE)
Q3_box <- quantile(SIO2, 0.75, na.rm = TRUE)
IQR_box <- IQR(SIO2, na.rm = TRUE)

# Outliers 
lim_inf_outlier_box <- Q1_box - 1.5 * IQR_box
lim_sup_outlier_box <- Q3_box + 1.5 * IQR_box
outliers_box <- SIO2[SIO2 < lim_inf_outlier_box | SIO2 > lim_sup_outlier_box]
n_outliers_box <- length(outliers_box)

boxplot(SIO2,
        horizontal = TRUE,
        col = "lightblue",
        main = "Gráfica 8: Distribución de SIO2 de la muestra \n con detección de valores atípicos",
        xlab = "Sílice (%)",
        ylab = "",
        border = "darkblue")

# Añadir media
points(media_box, 1, pch = 23, bg = "red", cex = 1.2)

legend("topright",
       legend = c(paste("Media:", round(media_box, 2)),
                 paste("Mediana:", round(mediana_box, 2)),
                 paste("Q1:", round(Q1_box, 2)),
                 paste("Q3:", round(Q3_box, 2)),
                 paste("Outliers:", n_outliers_box)),
       bty = "n",
       cex = 0.8)

summary(SIO2)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   42.46   57.66   63.06   63.87   70.71   80.00
# Histograma con Boxplot superpuesto
hist(
  SIO2,
  breaks = h$breaks,
  col = "lightgray",
  border = "white",
  main = "Gráfica Nº9: Histograma con diagrama de caja superpuesto\nde SIO2 de depósitos minerales de Estados Unidos",
  xlab = "Sílice (%)",
  ylab = "Cantidad"
)

boxplot(
  SIO2,
  horizontal = TRUE,
  at = par("usr")[4] * 0.70,
  add = TRUE,
  xaxt = "n",
  yaxt = "n",
  boxwex = par("usr")[4] * 0.30,
  col = rgb(0.2, 0.2, 1, 0.30),
  border = "blue"
)

INDICADORES ESTADISTICOS Y OUTLIERS

# Calcular indicadores de tendencia central
minimo <- min(SIO2)
maximo <- max(SIO2)
rango <- maximo - minimo
media <- mean(SIO2)
mediana <- median(SIO2)
moda <- as.numeric(names(sort(table(SIO2), decreasing = TRUE)[1]))

# Crear tabla
tendencia_central <- data.frame(
  Indicador = c("Mínimo", "Media", "Mediana", "Moda", "Máximo", "Rango"),
  Valor = c(
    round(minimo, 4),
    round(media, 4),
    round(mediana, 4),
    round(moda, 4),
    round(maximo, 4),
    round(rango, 4)
  ),
  Unidad = rep("%", 6),
  Interpretación = c(
    "Valor mínimo observado",
    "Promedio de todos los valores",
    "Valor que divide la muestra en dos partes iguales",
    "Valor más frecuente en la muestra",
    "Valor máximo observado",
    "Diferencia entre el máximo y el mínimo"
  )
)

kable(tendencia_central,
      caption = "Tabla 3: Indicadores de Tendencia Central para la variable SIO2_WT_PE",
      align = "l")
Tabla 3: Indicadores de Tendencia Central para la variable SIO2_WT_PE
Indicador Valor Unidad Interpretación
Mínimo 42.4600 % Valor mínimo observado
Media 63.8746 % Promedio de todos los valores
Mediana 63.0650 % Valor que divide la muestra en dos partes iguales
Moda 71.0700 % Valor más frecuente en la muestra
Máximo 80.0000 % Valor máximo observado
Rango 37.5400 % Diferencia entre el máximo y el mínimo 
# Calcular indicadores de dispersión
varianza <- var(SIO2)
desv_est <- sd(SIO2)
CV <- (desv_est / media) * 100

# Interpretación del CV
if(CV < 15) {
  interpretacion_CV <- "BAJA (CV < 15%)"
} else if(CV < 30) {
  interpretacion_CV <- "MODERADA (15% ≤ CV < 30%)"
} else {
  interpretacion_CV <- "ALTA (CV ≥ 30%)"
}

# Crear tabla
dispersion <- data.frame(
  Indicador = c("Varianza", "Desviación Estándar", "Coeficiente de Variación"),
  Valor = c(
    round(varianza, 4),
    round(desv_est, 4),
    paste0(round(CV, 2), "%")
  ),
  Unidad = c("%²", "%", ""),
  Interpretación = c(
    "Medida de dispersión al cuadrado",
    "Dispersión promedio respecto a la media",
    paste("Dispersión relativa:", interpretacion_CV)
  )
)

kable(dispersion,
      caption = "Tabla 4: Indicadores de Dispersión para la variable SIO2_WT_PE",
      align = "l")
Tabla 4: Indicadores de Dispersión para la variable SIO2_WT_PE
Indicador Valor Unidad Interpretación
Varianza 54.6243 %² Medida de dispersión al cuadrado
Desviación Estándar 7.3908 % Dispersión promedio respecto a la media
Coeficiente de Variación 11.57% Dispersión relativa: BAJA (CV < 15%) 
# Calcular indicadores de posición
cuartiles <- quantile(SIO2)
Q1 <- cuartiles[2]
Q2 <- cuartiles[3]
Q3 <- cuartiles[4]
IQR_val <- IQR(SIO2)

# Detección de outliers
lim_inf_outlier <- Q1 - 1.5 * IQR_val
lim_sup_outlier <- Q3 + 1.5 * IQR_val
outliers <- SIO2[SIO2 < lim_inf_outlier | SIO2 > lim_sup_outlier]
n_outliers <- length(outliers)
porc_outliers <- round((n_outliers / n) * 100, 2)

# Crear tabla
posicion <- data.frame(
  Indicador = c("Cuartil 1 (Q1)", "Cuartil 2 (Q2 - Mediana)", "Cuartil 3 (Q3)", 
                "Rango Intercuartílico (IQR)", "Límite Inferior Outliers", 
                "Límite Superior Outliers", "Número de Outliers"),
  Valor = c(
    round(Q1, 4),
    round(Q2, 4),
    round(Q3, 4),
    round(IQR_val, 4),
    round(lim_inf_outlier, 4),
    round(lim_sup_outlier, 4),
    paste0(n_outliers, " (", porc_outliers, "%)")
  ),
  Unidad = c(rep("%", 6), "observaciones"),
  Interpretación = c(
    "25% de datos por debajo de este valor",
    "50% de datos por debajo de este valor (coincide con mediana)",
    "75% de datos por debajo de este valor",
    "Rango del 50% central de datos (Q3 - Q1)",
    "Límite inferior para detección de valores atípicos",
    "Límite superior para detección de valores atípicos",
    "Cantidad y porcentaje de valores atípicos detectados"
  )
)

kable(posicion,
      caption = "Tabla 5: Indicadores de Posición y detección de outliers en SIO2_WT_PE",
      align = "l")
Tabla 5: Indicadores de Posición y detección de outliers en SIO2_WT_PE
Indicador Valor Unidad Interpretación
Cuartil 1 (Q1) 57.6575 % 25% de datos por debajo de este valor
Cuartil 2 (Q2 - Mediana) 63.065 % 50% de datos por debajo de este valor (coincide con mediana)
Cuartil 3 (Q3) 70.71 % 75% de datos por debajo de este valor
Rango Intercuartílico (IQR) 13.0525 % Rango del 50% central de datos (Q3 - Q1)
Límite Inferior Outliers 38.0788 % Límite inferior para detección de valores atípicos
Límite Superior Outliers 90.2887 % Límite superior para detección de valores atípicos
Número de Outliers 0 (0%) observaciones Cantidad y porcentaje de valores atípicos detectados 
# Calcular coeficiente de asimetría de Fisher
asimetria <- moments::skewness(SIO2)

if(abs(asimetria) < 0.5) {
  interpretacion_asimetria <- "Distribución simétrica"
} else if(asimetria > 0) {
  interpretacion_asimetria <- "Asimetría positiva (sesgo a la derecha)"
} else {
  interpretacion_asimetria <- "Asimetría negativa (sesgo a la izquierda)"
}

# Curtosis
curtosis <- moments::kurtosis(SIO2) - 3

if(abs(curtosis) < 0.5) {
  interpretacion_curtosis <- "Distribución mesocúrtica (normal)"
} else if(curtosis > 0) {
  interpretacion_curtosis <- "Distribución leptocúrtica (picuda)"
} else {
  interpretacion_curtosis <- "Distribución platicúrtica (aplanada)"
}

# Crear tabla
forma <- data.frame(
  Indicador = c("Coeficiente de Asimetría (Fisher)", "Interpretación Asimetría",
                "Coeficiente de Curtosis (Exceso)", "Interpretación Curtosis"),
  Valor = c(
    round(asimetria, 4),
    interpretacion_asimetria,
    round(curtosis, 4),
    interpretacion_curtosis
  ),
  Fórmula = c(
    "g₁ = E[(X-μ)³]/σ³",
    "|g₁| < 0.5: Simétrica; g₁ > 0: Positiva; g₁ < 0: Negativa",
    "g₂ = E[(X-μ)⁴]/σ⁴ - 3",
    "|g₂| < 0.5: Mesocúrtica; g₂ > 0: Leptocúrtica; g₂ < 0: Platicúrtica"
  )
)

kable(forma,
      caption = "Tabla 6: Indicadores de forma de la distribución de SIO2_WT_PE",
      align = "l")
Tabla 6: Indicadores de forma de la distribución de SIO2_WT_PE
Indicador Valor Fórmula
Coeficiente de Asimetría (Fisher) 0.0101 g₁ = E[(X-μ)³]/σ³
Interpretación Asimetría Distribución simétrica |g₁| < 0.5: Simétrica; g₁ > 0: Positiva; g₁ < 0: Negativa
Coeficiente de Curtosis (Exceso) -1.0922 g₂ = E[(X-μ)⁴]/σ⁴ - 3
Interpretación Curtosis Distribución platicúrtica (aplanada) |g₂| < 0.5: Mesocúrtica; g₂ > 0: Leptocúrtica; g₂ < 0: Platicúrtica

Conclusión

La variable SIO₂ presenta valores que fluctúan entre 42.46% y 80%, con valores en torno a la mediana de 63.07%. La desviación estándar de 7.39% y el coeficiente de variación de 11.57% indican un conjunto homogéneo. La distribución es simétrica y no presenta valores atípicos, lo que resulta favorable para la caracterización geoquímica, debido a que existe estabilidad en la concentración de sílice dentro de los depósitos analizados.