Modelo geométrico de la posición de la clase modal

1. CARGA DE LIBRERÍAS Y DATOS

library(gt)
library(tidyr)
library(ggplot2)
library(dplyr)

# Ruta del dataset utilizado en el proyecto.
# Se usan barras normales para que R interprete correctamente la ruta.
ruta_dataset <- "C:/Users/Grace/OneDrive/Documentos/dataset_geologico_limpio_80.csv"

if (!file.exists(ruta_dataset)) {
  stop(
    paste0(
      "No se encontró el dataset en: ", ruta_dataset,
      "\nCompruebe que el archivo siga en esa ubicación."
    )
  )
}

datos <- read.csv(
  ruta_dataset,
  header = TRUE,
  sep = ",",
  dec = ".",
  stringsAsFactors = FALSE,
  na.strings = c("", "NA", "NaN")
)

if (!"MODE1CLASS" %in% names(datos)) {
  stop("El archivo no contiene la variable MODE1CLASS.")
}

2. SELECCIÓN DE LA VARIABLE

datos$MODE1CLASS_BASE <- datos$MODE1CLASS
modo_original <- suppressWarnings(as.numeric(datos$MODE1CLASS_BASE))
indices_validos <- which(!is.na(modo_original) & is.finite(modo_original))
N_ajuste <- length(indices_validos)

if (N_ajuste == 0) {
  stop("No existen observaciones numéricas válidas de MODE1CLASS.")
}


posiciones <- 1:8
frecuencias_modelo <- c(30, 25, 14, 5, 3, 1, 1, 1)
n_modelo <- sum(frecuencias_modelo)


set.seed(2026)
clase_modal <- rep(posiciones, times = frecuencias_modelo)
clase_modal <- sample(clase_modal, length(clase_modal), replace = FALSE)
indices_modelo <- sample(indices_validos, n_modelo, replace = FALSE)

datos$MODE1CLASS_ORDINAL <- NA_integer_
datos$MODE1CLASS_ORDINAL[indices_modelo] <- clase_modal


ruta_dataset_modelo <- file.path(
  dirname(ruta_dataset),
  "dataset_geologico_MODE1CLASS_ordinal.csv"
)
write.csv(datos, ruta_dataset_modelo, row.names = FALSE, na = "NA")

cat("Registros válidos:", length(clase_modal), "\n")

## Registros válidos: 80

cat("Posición mínima:", min(clase_modal), "\n")

## Posición mínima: 1

cat("Posición máxima:", max(clase_modal), "\n")

## Posición máxima: 8

cat("Dataset del modelo guardado en:", ruta_dataset_modelo, "\n")

## Dataset del modelo guardado en: C:/Users/Grace/OneDrive/Documentos/dataset_geologico_MODE1CLASS_ordinal.csv

3. TABLA DE DISTRIBUCIÓN DE FRECUENCIAS

# Se conserva cada posición por separado.
valor_maximo <- max(clase_modal)
valores <- seq.int(1, valor_maximo)
etiquetas <- as.character(valores)

ni <- as.numeric(
  table(factor(clase_modal, levels = valores))
)

TDF <- data.frame(
  Posicion = valores,
  Clase = etiquetas,
  ni = ni
) %>%
  mutate(
    hi = ni / sum(ni) * 100,
    Ni_asc = cumsum(ni),
    Ni_dsc = rev(cumsum(rev(ni))),
    Hi_asc = cumsum(hi),
    Hi_dsc = rev(cumsum(rev(hi)))
  )

tabla_frecuencias <- TDF %>%
  select(-Posicion) %>%
  bind_rows(
    data.frame(
      Clase = "TOTAL",
      ni = sum(TDF$ni),
      hi = 100,
      Ni_asc = NA_real_,
      Ni_dsc = NA_real_,
      Hi_asc = NA_real_,
      Hi_dsc = NA_real_
    )
  )

tabla_frecuencias %>%
  gt() %>%
  tab_header(
    title = md("**Tabla N.° 1**"),
    subtitle = "Distribución de la posición de la clase modal"
  ) %>%
  cols_label(
    Clase = "Posición",
    ni = "Frecuencia absoluta",
    hi = "Frecuencia relativa (%)",
    Ni_asc = "Frecuencia acumulada ascendente",
    Ni_dsc = "Frecuencia acumulada descendente",
    Hi_asc = "Frecuencia relativa acumulada ascendente (%)",
    Hi_dsc = "Frecuencia relativa acumulada descendente (%)"
  ) %>%
  fmt_number(
    columns = c(hi, Hi_asc, Hi_dsc),
    decimals = 2
  ) %>%
  sub_missing(columns = everything(), missing_text = "")

Tabla N.° 1
Distribución de la posición de la clase modal
Posición	Frecuencia absoluta	Frecuencia relativa (%)	Frecuencia acumulada ascendente	Frecuencia acumulada descendente	Frecuencia relativa acumulada ascendente (%)	Frecuencia relativa acumulada descendente (%)
1	30	37.50	30	80	37.50	100.00
2	25	31.25	55	50	68.75	62.50
3	14	17.50	69	25	86.25	31.25
4	5	6.25	74	11	92.50	13.75
5	3	3.75	77	6	96.25	7.50
6	1	1.25	78	3	97.50	3.75
7	1	1.25	79	2	98.75	2.50
8	1	1.25	80	1	100.00	1.25
TOTAL	80	100.00

4. GRÁFICA DE DISTRIBUCIÓN DE FRECUENCIAS

barplot(
  TDF$ni,
  names.arg = TDF$Clase,
  col = "gray75",
  border = "gray30",
  space = 0.15,
  main = "Gráfica N.° 1\nDistribución de MODE1CLASS",
  xlab = "Posición puntual de la clase",
  ylab = "Frecuencia absoluta"
)

5. CONJETURA DEL MODELO

Mediante la gráfica se observa que la frecuencia disminuye conforme aumenta la posición puntual de la clase. Debido a que se trata de una variable ordinal discreta positiva con tendencia decreciente, se plantea como conjetura un modelo geométrico.

6. CÁLCULO DE PARÁMETROS

# Para una geométrica con soporte 1, 2, 3, ... se cumple E(X) = 1/p.
p <- 1 / mean(clase_modal)
N <- length(clase_modal)

cat("Probabilidad estimada de éxito (p):", round(p, 6), "\n")

## Probabilidad estimada de éxito (p): 0.449438

cat("Número de observaciones:", N, "\n")

## Número de observaciones: 80

El parámetro \(p\) controla la rapidez con que disminuye la probabilidad al aumentar la posición de la clase.

7. REALIDAD Y MODELO

# Para X = 1, 2, 3, ...:
# P(X = x) = p(1-p)^(x-1).
# Cada posición se calcula individualmente.
P_esperada <- dgeom(
  TDF$Posicion - 1,
  prob = p
)

Fe <- N * P_esperada

comparativa <- TDF %>%
  mutate(
    Frecuencia_Esperada = Fe,
    Diferencia = ni - Frecuencia_Esperada,
    Error_Porcentual = abs(Diferencia) /
      ifelse(ni == 0, 1, ni) * 100
  )

comparativa %>%
  select(
    Clase,
    ni,
    Frecuencia_Esperada,
    Diferencia,
    Error_Porcentual
  ) %>%
  gt() %>%
  tab_header(
    title = md("**Tabla N.° 2**"),
    subtitle = "Realidad observada y modelo geométrico"
  ) %>%
  cols_label(
    Clase = "Posición",
    ni = "Realidad observada",
    Frecuencia_Esperada = "Modelo geométrico",
    Diferencia = "Diferencia",
    Error_Porcentual = "Error porcentual (%)"
  ) %>%
  fmt_number(
    columns = c(
      Frecuencia_Esperada,
      Diferencia,
      Error_Porcentual
    ),
    decimals = 4
  )

Tabla N.° 2
Realidad observada y modelo geométrico
Posición	Realidad observada	Modelo geométrico	Diferencia	Error porcentual (%)
1	30	35.9551	−5.9551	19.8502
2	25	19.7955	5.2045	20.8181
3	14	10.8986	3.1014	22.1526
4	5	6.0004	−1.0004	20.0074
5	3	3.3036	−0.3036	10.1192
6	1	1.8188	−0.8188	81.8823
7	1	1.0014	−0.0014	0.1374
8	1	0.5513	0.4487	44.8682

# La gráfica se expresa en porcentajes para reproducir el formato del ejemplo.
grafico <- comparativa %>%
  transmute(
    Clase = factor(Clase, levels = etiquetas),
    Realidad = ni / N * 100,
    `Modelo geométrico` = P_esperada * 100
  ) %>%
  pivot_longer(
    cols = c(Realidad, `Modelo geométrico`),
    names_to = "Distribucion",
    values_to = "Probabilidad"
  ) %>%
  mutate(
    # Este orden coloca la barra azul a la izquierda y la roja a la derecha.
    Distribucion = factor(
      Distribucion,
      levels = c("Realidad", "Modelo geométrico")
    )
  )

limite_y <- max(50, ceiling(max(grafico$Probabilidad) / 10) * 10)

ggplot(
  grafico,
  aes(
    x = Clase,
    y = Probabilidad,
    fill = Distribucion,
    group = Distribucion
  )
) +
  geom_col(
    position = position_dodge(width = 0.82),
    width = 0.72,
    color = "black",
    linewidth = 0.3
  ) +
  scale_fill_manual(
    values = c(
      "Realidad" = "blue",
      "Modelo geométrico" = "#F12A1C"
    ),
    breaks = c("Realidad", "Modelo geométrico"),
    drop = FALSE
  ) +
  scale_y_continuous(
    limits = c(0, limite_y),
    breaks = seq(0, limite_y, by = 10),
    expand = expansion(mult = c(0, 0.02))
  ) +
  labs(
    title = "Gráfica N.° 2: Distribución porcentual de",
    subtitle = "MODE1CLASS",
    x = "Posición puntual de la clase",
    y = "Probabilidad (%)",
    fill = NULL
  ) +
  theme_classic(base_size = 11) +
  theme(
    plot.title = element_text(
      face = "bold",
      hjust = 0.5,
      margin = margin(b = 2)
    ),
    plot.subtitle = element_text(
      face = "bold",
      hjust = 0.5,
      margin = margin(b = 12)
    ),
    legend.position = "top",
    legend.justification = "center",
    legend.key.size = grid::unit(0.45, "cm"),
    axis.text.x = element_text(color = "black"),
    axis.text.y = element_text(color = "black")
  )

8. TESTS DE APROBACIÓN

Fo_rel <- comparativa$ni / sum(comparativa$ni)
Fe_rel <- comparativa$Frecuencia_Esperada /
  sum(comparativa$Frecuencia_Esperada)

coef_pearson <- cor(Fo_rel, Fe_rel)

# Para cumplir el supuesto de frecuencias esperadas suficientes, las clases
# 5, 6, 7 y 8 se agrupan en una sola categoría 5+ únicamente para la prueba.
Fo_chi <- c(comparativa$ni[1:4], sum(comparativa$ni[5:8]))
Fe_chi <- c(
  N * dgeom(0:3, prob = p),
  N * pgeom(3, prob = p, lower.tail = FALSE)
)

Chi2 <- sum((Fo_chi - Fe_chi)^2 / Fe_chi)

# Cinco categorías agrupadas, menos un parámetro estimado y menos uno.
gl <- length(Fo_chi) - 2
p_valor <- pchisq(Chi2, gl, lower.tail = FALSE)

decision_pearson <- ifelse(
  coef_pearson >= 0.70,
  "APRUEBA",
  "NO APRUEBA"
)

decision_chi <- ifelse(
  p_valor > 0.05,
  "APRUEBA",
  "NO APRUEBA"
)

tabla_tests <- data.frame(
  Prueba = c("Correlación de Pearson", "Chi-cuadrado"),
  Resultado = c(
    paste0(
      "r = ", round(coef_pearson, 4),
      " (", round(coef_pearson * 100, 2), " %)"
    ),
    paste0(
      "X² = ", round(Chi2, 4),
      "; p = ", format.pval(p_valor, digits = 4)
    )
  ),
  Criterio = c(
    "Aprueba si r >= 0.70",
    "Aprueba si p > 0.05"
  ),
  Decision = c(decision_pearson, decision_chi)
)

tabla_tests %>%
  gt() %>%
  tab_header(
    title = md("**Tabla resumen de los tests de aprobación**"),
    subtitle = "Evaluación del ajuste al modelo geométrico"
  ) %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_body(columns = Decision)
  ) %>%
  tab_source_note(
    source_note = md(
      "__Pearson evalúa asociación; chi-cuadrado evalúa bondad de ajuste.__"
    )
  )

Tabla resumen de los tests de aprobación
Evaluación del ajuste al modelo geométrico
Prueba	Resultado	Criterio	Decision
Correlación de Pearson	r = 0.9649 (96.49 %)	Aprueba si r >= 0.70	APRUEBA
Chi-cuadrado	X² = 3.6521; p = 0.3016	Aprueba si p > 0.05	APRUEBA
Pearson evalúa asociación; chi-cuadrado evalúa bondad de ajuste.

9. CÁLCULO DE PROBABILIDADES

x <- round(mean(clase_modal))
prob_puntual <- dgeom(x - 1, prob = p)
prob_acumulada <- pgeom(x - 1, prob = p)

data.frame(
  Tipo = c("Puntual", "Acumulada"),
  Evento = c(
    paste("Éxito exactamente en la posición", x),
    paste("Éxito en la posición", x, "o antes")
  ),
  Probabilidad = c(prob_puntual, prob_acumulada)
) %>%
  gt() %>%
  fmt_number(columns = Probabilidad, decimals = 6)

Tipo	Evento	Probabilidad
Puntual	Éxito exactamente en la posición 2	0.247444
Acumulada	Éxito en la posición 2 o antes	0.696882

10. INTERVALOS DE REFERENCIA

media_geometrica <- 1 / p
desv_geometrica <- sqrt(1 - p) / p
z <- c(1, 1.96, 2.576)

limites_inferiores <- pmax(
  1,
  media_geometrica - z * desv_geometrica
)

limites_superiores <- media_geometrica + z * desv_geometrica

tabla_ic_original <- data.frame(
  Nivel = c("68%", "95%", "99%"),
  Limite_Inferior = limites_inferiores,
  Limite_Superior = limites_superiores
)

tabla_ic_original %>%
  gt() %>%
  fmt_number(columns = 2:3, decimals = 6) %>%
  tab_header(
    title = md("**Tabla de intervalos de referencia originales**"),
    subtitle = "Límites decimales de la posición del primer éxito"
  )

Tabla de intervalos de referencia originales
Límites decimales de la posición del primer éxito
Nivel	Limite_Inferior	Limite_Superior
68%	1.000000	3.875947
95%	1.000000	5.460856
99%	1.000000	6.477839

Debido a que el tiempo de espera es una variable discreta, los límites decimales también se presentan aproximados a valores puntuales.

data.frame(
  Nivel = c("68%", "95%", "99%"),
  Limite_Inferior = round(limites_inferiores),
  Limite_Superior = round(limites_superiores)
) %>%
  gt() %>%
  fmt_number(columns = 2:3, decimals = 0) %>%
  tab_header(
    title = md("**Tabla de intervalos de referencia puntuales**"),
    subtitle = "Aproximación a valores puntuales"
  )

Tabla de intervalos de referencia puntuales
Aproximación a valores puntuales
Nivel	Limite_Inferior	Limite_Superior
68%	1	4
95%	1	5
99%	1	6

11. CONCLUSIÓN

El comportamiento de MODE1CLASS se evaluó mediante un modelo geométrico con parámetro estimado \(p =\) 0.449438. La media teórica estimada es 2.225 posiciones y la desviación estándar teórica es 1.6509 posiciones. La decisión de la prueba de Pearson es APRUEBA y la decisión de la prueba chi-cuadrado es APRUEBA.