SIO2

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

## Loading required package: readr

if (!require("dplyr")) install.packages("dplyr")

## Loading required package: dplyr

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

if (!require("knitr")) install.packages("knitr")

## Loading required package: knitr

if (!require("moments")) install.packages("moments")

## Loading required package: moments

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

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

## Rows: 2500 Columns: 2

## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): EARTH_MA_2
## dbl (1): SIO2_WT_PE
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

cat("✓ Datos cargados exitosamente\n")

## ✓ Datos cargados exitosamente

cat("✓ Dimensiones:", nrow(datos), "observaciones\n")

## ✓ Dimensiones: 2500 observaciones

#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)

#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_WT_PE (Porcentaje de Sílice)")

Tabla 1: Parámetros de Sturges para la variable SIO2_WT_PE (Porcentaje de Sílice)

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

#Intervalos y frecuencias
Li <- seq(min(SIO2), max(SIO2) - A, by = A)
Ls <- seq(min(SIO2) + A, max(SIO2), by = A)
MC <- (Li + Ls) / 2

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

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)

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 de SIO2_WT_PE (Porcentaje de Sílice)")

Tabla 2: Distribución de frecuencias de SIO2_WT_PE (Porcentaje de Sílice)

Clase	Li	Ls	MC	ni	hi	Niasc	Nidsc	Hiasc	Hidsc
1	42.46	45.59	44.02	5	0.20	5	2499	0.20	99.96
2	45.59	48.72	47.15	21	0.84	26	2494	1.04	99.76
3	48.72	51.84	50.28	50	2.00	76	2473	3.04	98.92
4	51.84	54.97	53.41	190	7.60	266	2423	10.64	96.92
5	54.97	58.10	56.54	427	17.08	693	2233	27.72	89.32
6	58.10	61.23	59.67	389	15.56	1082	1806	43.28	72.24
7	61.23	64.36	62.79	259	10.36	1341	1417	53.64	56.68
8	64.36	67.49	65.92	218	8.72	1559	1158	62.36	46.32
9	67.49	70.61	69.05	302	12.08	1861	940	74.44	37.60
10	70.62	73.74	72.18	391	15.64	2252	638	90.08	25.52
11	73.74	76.87	75.31	207	8.28	2459	247	98.36	9.88
12	76.87	80.00	78.44	40	1.60	2499	40	99.96	1.60

#Histograma con media y mediana
media <- mean(SIO2)
mediana <- median(SIO2)

hist(
  SIO2,
  breaks = seq(min(SIO2), max(SIO2) + A, A),
  col = "steelblue",
  main = "Distribución del Porcentaje de Sílice (SIO2_WT_PE)",
  xlab = "Sílice (%)",
  ylab = "Frecuencia"
)

abline(v = media, col = "red", lwd = 2, lty = 2)
abline(v = mediana, col = "blue", lwd = 2, lty = 2)

legend("topright",
       legend = c(paste("Media =", round(media, 2)),
                  paste("Mediana =", round(mediana, 2))),
       col = c("red", "blue"),
       lty = 2,
       lwd = 2,
       bty = "n")

#Histograma global
hist(
  SIO2,
  breaks = seq(min(SIO2), max(SIO2) + A, A),
  col = "lightcoral",
  main = "Distribución global de SIO2_WT_PE",
  xlab = "Sílice (%)",
  ylab = "Frecuencia"
)
abline(v = media, col = "red", lwd = 2, lty = 2)
abline(v = mediana, col = "blue", lwd = 2, lty = 2)

legend("topright",
       legend = c(paste("Media =", round(media, 2)),
                  paste("Mediana =", round(mediana, 2))),
       col = c("red", "blue"),
       lty = 2,
       lwd = 2,
       bty = "n")

# Gráfica 3: Distribución porcentual
intervalos <- paste(round(Li, 2), round(Ls, 2), sep = " - ")
pos_media <- (media - min(Li)) / A
pos_mediana <- (mediana - min(Li)) / A

barplot(
  hi,
  names.arg = intervalos,
  col = "darkorange",
  ylim = c(0, max(hi) + 5),
  space = 0,
  cex.names = 0.7,
  ylab = "Frecuencia Relativa (%)",
  xlab = "Intervalos de Sílice (%)",
  main = "Gráfica 3: Distribución porcentual de SIO2_WT_PE por intervalos",
  las = 2,
  border = "white"
)

abline(v = pos_media, col = "red", lwd = 2, lty = 2)
abline(v = pos_mediana, col = "blue", lwd = 2, lty = 2)

legend("topright",
       legend = c(paste("Media =", round(media, 2)),
                  paste("Mediana =", round(mediana, 2))),
       col = c("red", "blue"),
       lty = 2,
       lwd = 2,
       bty = "n")

# Calcular media y mediana 
media_rel_global <- mean(SIO2)
mediana_rel_global <- median(SIO2)

# Calcular posición exacta en el barplot
pos_media_bar_global <- (media_rel_global - min(Li)) / A
pos_mediana_bar_global <- (mediana_rel_global - min(Li)) / A

barplot(
  hi,   
  names.arg = intervalos,
  col = "mediumpurple",
  ylim = c(0, 100),
  cex.names = 0.7,
  space = 0, 
  ylab = "Frecuencia Relativa Acumulada (%)",
  xlab = "Intervalos de Sílice (%)",
  main = "Gráfica 4: Distribución acumulada porcentual de SIO2_WT_PE",
  las = 2,
  border = "white"
)

# Líneas verticales de media y mediana
abline(v = pos_media_bar_global, col = "red", lwd = 2, lty = 2)
abline(v = pos_mediana_bar_global, col = "blue", lwd = 2, lty = 2)

# Leyenda
legend("topright",
       legend = c(paste("Media =", round(media_rel_global, 2)),
                  paste("Mediana =", round(mediana_rel_global, 2))),
       col = c("red", "blue"),
       lty = 2,
       lwd = 2,
       bty = "n")

# Gráfica 5: Ojiva Absoluta
plot(Ls, Niasc, 
     type = "o", 
     pch = 16, 
     col = "blue",
     main = "Gráfica 5: Ojiva de Frecuencias Absolutas Acumuladas \nde SIO2_WT_PE",
     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 6: Ojiva Relativa
plot(Ls, Hiasc, 
     type = "o", 
     pch = 16, 
     col = "blue",
     main = "Gráfica 6: Ojiva de Frecuencias Relativas Acumuladas \nde SIO2_WT_PE",
     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)

# Calcular todas las variables necesarias para el boxplot
media_box <- mean(SIO2)
mediana_box <- median(SIO2)
Q1_box <- quantile(SIO2, 0.25)
Q3_box <- quantile(SIO2, 0.75)
IQR_box <- IQR(SIO2)

# Calcular 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 7: Distribución de SIO2_WT_PE \n con detección de valores atípicos",
        xlab = "Sílice (%)",
        ylab = "",
        border = "darkblue")

# Añadir estadísticos importantes
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)

# 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