PROFUNDIDAD DE LA MUESTRA

#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_with_depth.csv"
datos <- read_csv(ruta)

## Rows: 2500 Columns: 3

## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): EARTH_MA_2
## dbl (1): TOP_DEPTH_M
## lgl (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 Profundidad
PROFUNDIDAD_raw <- datos$TOP_DEPTH_M
PROFUNDIDAD <- as.numeric(PROFUNDIDAD_raw)
PROFUNDIDAD <- PROFUNDIDAD[!is.na(PROFUNDIDAD)]
# Nota: En profundidad, el valor 0 puede ser válido (superficie), 
# pero siguiendo el esqueleto previo de limpieza:
PROFUNDIDAD <- PROFUNDIDAD[PROFUNDIDAD != 0]

n <- length(PROFUNDIDAD)

#Parámetros de Sturges
R <- max(PROFUNDIDAD) - min(PROFUNDIDAD)
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 TOP_DEPTH_M (Profundidad de la muestra)")

Tabla 1: Parámetros de Sturges para la variable TOP_DEPTH_M (Profundidad de la muestra)

Parámetro	Valor
Rango (R)	805.4700
Número de datos (n)	2499.0000
Número de intervalos (k)	12.0000
Amplitud de clase (A)	67.1225

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

ni <- numeric(k)
for (i in 1:k) {
  if (i == k) {
    ni[i] <- sum(PROFUNDIDAD >= Li[i] & PROFUNDIDAD <= Ls[i])
  } else {
    ni[i] <- sum(PROFUNDIDAD >= Li[i] & PROFUNDIDAD < 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 TOP_DEPTH_M (Profundidad en metros)")

Tabla 2: Distribución de frecuencias de TOP_DEPTH_M (Profundidad en metros)

Clase	Li	Ls	MC	ni	hi	Niasc	Nidsc	Hiasc	Hidsc
1	10.47	77.59	44.03	81	3.24	81	2499	3.24	100.00
2	77.59	144.72	111.15	222	8.88	303	2418	12.12	96.76
3	144.72	211.84	178.28	240	9.60	543	2196	21.72	87.88
4	211.84	278.96	245.40	341	13.65	884	1956	35.37	78.28
5	278.96	346.08	312.52	439	17.57	1323	1615	52.94	64.63
6	346.08	413.21	379.64	391	15.65	1714	1176	68.59	47.06
7	413.21	480.33	446.77	369	14.77	2083	785	83.36	31.41
8	480.33	547.45	513.89	222	8.88	2305	416	92.24	16.64
9	547.45	614.57	581.01	119	4.76	2424	194	97.00	7.76
10	614.57	681.69	648.13	61	2.44	2485	75	99.44	3.00
11	681.70	748.82	715.26	10	0.40	2495	14	99.84	0.56
12	748.82	815.94	782.38	4	0.16	2499	4	100.00	0.16

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

hist(
  PROFUNDIDAD,
  breaks = seq(min(PROFUNDIDAD), max(PROFUNDIDAD) + A, A),
  col = "steelblue",
  main = "Distribución de la Profundidad de la Muestra (TOP_DEPTH_M)",
  xlab = "Profundidad (m)",
  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(
  PROFUNDIDAD,
  breaks = seq(min(PROFUNDIDAD), max(PROFUNDIDAD) + A, A),
  col = "lightcoral",
  main = "Distribución global de TOP_DEPTH_M",
  xlab = "Profundidad (m)",
  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 Profundidad (m)",
  main = "Gráfica 3: Distribución porcentual de TOP_DEPTH_M 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(PROFUNDIDAD)
mediana_rel_global <- median(PROFUNDIDAD)

# 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 Profundidad (m)",
  main = "Gráfica 4: Distribución acumulada porcentual de TOP_DEPTH_M",
  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 TOP_DEPTH_M",
     xlab = "Profundidad (m)",
     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 TOP_DEPTH_M",
     xlab = "Profundidad (m)",
     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(PROFUNDIDAD)
mediana_box <- median(PROFUNDIDAD)
Q1_box <- quantile(PROFUNDIDAD, 0.25)
Q3_box <- quantile(PROFUNDIDAD, 0.75)
IQR_box <- IQR(PROFUNDIDAD)

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

boxplot(PROFUNDIDAD,
        horizontal = TRUE,
        col = "lightblue",
        main = "Gráfica 7: Distribución de TOP_DEPTH_M \n con detección de valores atípicos",
        xlab = "Profundidad (m)",
        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(PROFUNDIDAD)
maximo <- max(PROFUNDIDAD)
rango <- maximo - minimo
media <- mean(PROFUNDIDAD)
mediana <- median(PROFUNDIDAD)
moda <- as.numeric(names(sort(table(PROFUNDIDAD), 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("m", 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 TOP_DEPTH_M",
      align = "l")

Tabla 3: Indicadores de Tendencia Central para la variable TOP_DEPTH_M

Indicador	Valor	Unidad	Interpretación
Mínimo	10.4700	m	Valor mínimo observado
Media	335.3845	m	Promedio de todos los valores
Mediana	335.9700	m	Valor que divide la muestra en dos partes iguales
Moda	233.2000	m	Valor más frecuente en la muestra
Máximo	815.9400	m	Valor máximo observado
Rango	805.4700	m	Diferencia entre el máximo y el mínimo

# Calcular indicadores de dispersión
varianza <- var(PROFUNDIDAD)
desv_est <- sd(PROFUNDIDAD)
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("m²", "m", ""),
  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 TOP_DEPTH_M",
      align = "l")

Tabla 4: Indicadores de Dispersión para la variable TOP_DEPTH_M

Indicador	Valor	Unidad	Interpretación
Varianza	21428.8446	m²	Medida de dispersión al cuadrado
Desviación Estándar	146.3859	m	Dispersión promedio respecto a la media
Coeficiente de Variación	43.65%		Dispersión relativa: ALTA (CV ≥ 30%)

# Calcular indicadores de posición
cuartiles <- quantile(PROFUNDIDAD)
Q1 <- cuartiles[2]
Q2 <- cuartiles[3]
Q3 <- cuartiles[4]
IQR_val <- IQR(PROFUNDIDAD)

# Detección de outliers
lim_inf_outlier <- Q1 - 1.5 * IQR_val
lim_sup_outlier <- Q3 + 1.5 * IQR_val
outliers <- PROFUNDIDAD[PROFUNDIDAD < lim_inf_outlier | PROFUNDIDAD > 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("m", 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 TOP_DEPTH_M",
      align = "l")

Tabla 5: Indicadores de Posición y detección de outliers en TOP_DEPTH_M

Indicador	Valor	Unidad	Interpretación
Cuartil 1 (Q1)	231.895	m	25% de datos por debajo de este valor
Cuartil 2 (Q2 - Mediana)	335.97	m	50% de datos por debajo de este valor (coincide con mediana)
Cuartil 3 (Q3)	438.4	m	75% de datos por debajo de este valor
Rango Intercuartílico (IQR)	206.505	m	Rango del 50% central de datos (Q3 - Q1)
Límite Inferior Outliers	-77.8625	m	Límite inferior para detección de valores atípicos
Límite Superior Outliers	748.1575	m	Límite superior para detección de valores atípicos
Número de Outliers	4 (0.16%)	observaciones	Cantidad y porcentaje de valores atípicos detectados

# Calcular coeficiente de asimetría de Fisher
asimetria <- moments::skewness(PROFUNDIDAD)

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(PROFUNDIDAD) - 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 TOP_DEPTH_M",
      align = "l")

Tabla 6: Indicadores de forma de la distribución de TOP_DEPTH_M

Indicador	Valor	Fórmula
Coeficiente de Asimetría (Fisher)	0.0728	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)	-0.5337	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