library(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
library(stringr)
library(gt)

# =========================================================
# IMPORTAR BASE
# =========================================================

datos <- read.csv("D:/dataset_CMIO_geologico.csv")

# =========================================================
# CREAR VARIABLE ORDINAL
# =========================================================

df_deposito <- data.frame(
  deposito = toupper(trimws(datos$Deposit_Type))
)

# =========================================================
# RECLASIFICACIÓN ORDINAL
# =========================================================

df_deposito$deposito <- case_when(

  toupper(df_deposito$deposito) %in% c("EPITHERMAL")
    ~ "Bajo",

  toupper(df_deposito$deposito) %in% c("PORPHYRY", "VMS")
    ~ "Medio",

  toupper(df_deposito$deposito) %in% c("SEDIMENT-HOSTED", "SKARN")
    ~ "Alto",

  TRUE ~ NA_character_
)

# =========================================================
# ORDEN LÓGICO
# =========================================================

orden_deposit <- c(
  "Bajo",
  "Medio",
  "Alto"
)

# =========================================================
# CONVERTIR A FACTOR ORDENADO
# =========================================================

df_deposito$deposito <- factor(
  df_deposito$deposito,
  levels = orden_deposit,
  ordered = TRUE
)

# =========================================================
# TABLA DE FRECUENCIAS
# =========================================================

ni <- table(df_deposito$deposito)

hi <- round(prop.table(ni), 4)

P <- round(hi * 100, 2)

tabla_finaldeposito <- data.frame(

  Nivel_Deposito = names(ni),

  ni = as.numeric(ni),

  hi = as.numeric(hi),

  P = as.numeric(P)
)

# =========================================================
# FILA TOTAL
# =========================================================

fila_total <- data.frame(

  Nivel_Deposito = "TOTAL",

  ni = sum(tabla_finaldeposito$ni),

  hi = sum(tabla_finaldeposito$hi),

  P = sum(tabla_finaldeposito$P)
)

tabla_finaldeposito <- rbind(
  tabla_finaldeposito,
  fila_total
)

tabla_finaldeposito
##   Nivel_Deposito   ni     hi      P
## 1           Bajo  497 0.1988  19.88
## 2          Medio 1007 0.4028  40.28
## 3           Alto  996 0.3984  39.84
## 4          TOTAL 2500 1.0000 100.00
tabla_deposito_gt <- tabla_finaldeposito %>%

  gt() %>%

  tab_header(

    title = md("**Tabla Nº1**"),

    subtitle = md("Distribución ordinal del tipo de depósito mineral")
  ) %>%

  tab_source_note(

    source_note = md("Autor: Grupo 2")
  ) %>%

  tab_style(

    style = cell_text(weight = "bold"),

    locations = cells_body(
      rows = Nivel_Deposito == "TOTAL"
    )
  )

tabla_deposito_gt
Tabla Nº1
Distribución ordinal del tipo de depósito mineral
Nivel_Deposito ni hi P
Bajo 497 0.1988 19.88
Medio 1007 0.4028 40.28
Alto 996 0.3984 39.84
TOTAL 2500 1.0000 100.00
Autor: Grupo 2 
P_global <- tabla_finaldeposito$P[
  tabla_finaldeposito$Nivel_Deposito != "TOTAL"
]

barplot(

  P_global,

  main = "Gráfica Nº1: Distribución ordinal del tipo de depósito",

  xlab = "Nivel ordinal",

  ylab = "Probabilidad (%)",

  col = "blue",

  names.arg = tabla_finaldeposito$Nivel_Deposito[
    tabla_finaldeposito$Nivel_Deposito != "TOTAL"
  ],

  ylim = c(0,100),

  las = 1
)

# =========================================================
# TAMAÑO MUESTRAL
# =========================================================

n <- sum(
  tabla_finaldeposito$ni[
    tabla_finaldeposito$Nivel_Deposito != "TOTAL"
  ]
)

n
## [1] 2500
# =========================================================
# FRECUENCIAS OBSERVADAS
# =========================================================

x <- tabla_finaldeposito$ni[
  tabla_finaldeposito$Nivel_Deposito != "TOTAL"
]

x
## [1]  497 1007  996
# =========================================================
# ESCALA ORDINAL
# =========================================================

X <- 1:length(x)

X
## [1] 1 2 3
# =========================================================
# MEDIA OBSERVADA
# =========================================================

media_observada <- sum(X * x) / n

media_observada
## [1] 2.1996
# =========================================================
# PARÁMETROS
# =========================================================

p <- media_observada / length(x)

p
## [1] 0.7332
q <- 1 - p

q
## [1] 0.2668
# =========================================================
# MODELO BINOMIAL
# =========================================================

P_binomial <- dbinom(
  X,
  size = length(x),
  prob = p
)

P_binomial
## [1] 0.1565725 0.4302808 0.3941553
# =========================================================
# FRECUENCIAS
# =========================================================

Fo <- (x / n) * 100

Fo
## [1] 19.88 40.28 39.84
Fe <- P_binomial * 100

Fe
## [1] 15.65725 43.02808 39.41553
barplot(

  rbind(Fo, Fe),

  beside = TRUE,

  col = c("skyblue", "blue"),

  names.arg = tabla_finaldeposito$Nivel_Deposito[
    tabla_finaldeposito$Nivel_Deposito != "TOTAL"
  ],

  main = "Gráfica Nº2: Comparación entre la realidad\ny el modelo binomial",

  ylab = "Probabilidad (%)",

  xlab = "Nivel ordinal",

  ylim = c(0,100)
)

legend(

  "topright",

  legend = c("Real", "Modelo"),

  fill = c("skyblue", "blue")
)

plot(

  Fo,
  Fe,

  xlim = c(0, max(Fo)),

  ylim = c(0, max(Fe)),

  main = "Gráfica Nº3: Correlación de frecuencias",

  xlab = "Frecuencia observada (%)",

  ylab = "Frecuencia esperada (%)",

  pch = 19,

  col = "darkblue"
)

abline(
  a = 0,
  b = 1,

  col = "red",

  lwd = 2
)

Correlacion <- cor(Fo, Fe) * 100

Correlacion
## [1] 99.47112
K <- length(x)

K
## [1] 3
gl <- K - 1

gl
## [1] 2
x2 <- sum((Fo - Fe)^2 / Fe)

x2
## [1] 1.318959
vc <- qchisq(0.99, gl)

vc
## [1] 9.21034
x2 < vc
## [1] TRUE
Variable <- c("Nivel ordinal del depósito")

tabla_resumen <- data.frame(

  Variable,

  round(Correlacion,2),

  round(x2,2),

  round(vc,2)
)

colnames(tabla_resumen) <- c(

  "Variable",

  "Test Pearson (%)",

  "Chi Cuadrado",

  "Umbral de aceptación"
)

library(knitr)

kable(

  tabla_resumen,

  format = "markdown",

  caption = "Tabla Nº2: Resumen de bondad del modelo"
  
)
Tabla Nº2: Resumen de bondad del modelo
Variable Test Pearson (%) Chi Cuadrado Umbral de aceptación
Nivel ordinal del depósito 99.47 1.32 9.21 
# =========================================================
# CÁLCULO DE PROBABILIDAD
# =========================================================

# Eliminar fila TOTAL

tabla_sin_total <- tabla_finaldeposito[
  tabla_finaldeposito$Nivel_Deposito != "TOTAL",
]

# =========================================================
# EXTRAER PROBABILIDAD DE LA CATEGORÍA
# =========================================================

prob_alto <- tabla_sin_total$P[
  tabla_sin_total$Nivel_Deposito == "Alto"
]

prob_alto
## [1] 39.84
# =========================================================
# GRÁFICO EXPLICATIVO
# =========================================================

plot(
  1,

  type = "n",

  axes = FALSE,

  xlab = "",

  ylab = ""
)

text(

  x = 1,
  y = 1,

  labels = paste(

    "Cálculo de probabilidad\n",
    "(Estimación general)\n\n",

    "¿Qué probabilidad existe de que un depósito\n",
    "mineral analizado en Estados Unidos pertenezca\n",
    "a un nivel de deposito geológico ALTO?\n\n",

    "Probabilidad = ",
    prob_alto,
    " (%)",

    sep = ""
  ),

  cex = 1.3,

  col = "black",

  font = 2
)