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)
library(readxl)

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

df_rasgo <- data.frame(
  rasgo = toupper(trimws(datos$TIPO_DE_RASGO_GEOLOGICO))
)

df_rasgo$rasgo <- case_when(

  # =========================================================
  # RASGO MENOR
  # =========================================================

  df_rasgo$rasgo %in% c(
    "FRACTURAS",
    "DIACLASAS",
    "VENILLAS SECUNDARIAS",
    "DISEMINACIÓN MINERAL",
    "TEXTURA PORFÍDICA",
    "TEXTURA GRANULAR",
    "BANDEAMIENTO LOCAL",
    "MICROFRACTURAS",
    "MINERALES ACCESORIOS VISIBLES"
  ) ~ "Rasgo menor",

  # =========================================================
  # RASGO INTERMEDIO
  # =========================================================

  df_rasgo$rasgo %in% c(
    "FALLA LOCAL",
    "ZONA DE CIZALLA",
    "CONTACTO LITOLÓGICO",
    "INTRUSIÓN LOCAL",
    "DIQUE",
    "SILL",
    "VETA PRINCIPAL",
    "ZONA MINERALIZADA",
    "BRECHA HIDROTERMAL",
    "STOCKWORK"
  ) ~ "Rasgo intermedio",

  # =========================================================
  # RASGO MAYOR
  # =========================================================

  df_rasgo$rasgo %in% c(
    "PROVINCIA GEOLÓGICA",
    "CINTURÓN MINERALIZADO",
    "SISTEMA DE FALLAS REGIONALES",
    "ARCO MAGMÁTICO",
    "BATOLITO O CUERPO INTRUSIVO REGIONAL",
    "TERRENO TECTÓNICO",
    "CUENCA GEOLÓGICA",
    "ZONA ESTRUCTURAL PRINCIPAL"
  ) ~ "Rasgo mayor",

  # =========================================================
  # OTROS
  # =========================================================

  TRUE ~ "Otros"
)

# =========================================================
# ORDEN ORDINAL
# =========================================================

orden_rasgo <- c(
  "Rasgo menor",
  "Rasgo intermedio",
  "Rasgo mayor",
  "Otros"
)

df_rasgo$rasgo <- factor(
  df_rasgo$rasgo,
  levels = orden_rasgo,
  ordered = TRUE
)

# =========================================================
# FRECUENCIAS
# =========================================================

ni <- table(df_rasgo$rasgo)

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

P <- round(hi * 100, 2)

tabla_finalrasgo <- data.frame(

  Tipo_Rasgo_Geologico = names(ni),

  ni = as.numeric(ni),

  hi = as.numeric(hi),

  P = as.numeric(P)
)

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

fila_total <- data.frame(

  Tipo_Rasgo_Geologico = "TOTAL",

  ni = sum(tabla_finalrasgo$ni),

  hi = sum(tabla_finalrasgo$hi),

  P = sum(tabla_finalrasgo$P)
)

tabla_finalrasgo <- rbind(
  tabla_finalrasgo,
  fila_total
)

tabla_finalrasgo
##   Tipo_Rasgo_Geologico   ni     hi      P
## 1          Rasgo menor  987 0.3948  39.48
## 2     Rasgo intermedio 1130 0.4520  45.20
## 3          Rasgo mayor  383 0.1532  15.32
## 4                Otros    0 0.0000   0.00
## 5                TOTAL 2500 1.0000 100.00
# =========================================================
# TABLA GT
# =========================================================

tabla_rasgo_gt <- tabla_finalrasgo %>%

  gt() %>%

  tab_header(

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

    subtitle = md(
      "Distribución ordinal del tipo de rasgo geológico"
    )
  ) %>%

  tab_source_note(

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

  tab_options(

    table.border.top.color = "black",

    table.border.bottom.color = "black",

    heading.border.bottom.color = "black",

    heading.border.bottom.width = px(2),

    column_labels.border.top.color = "black",

    column_labels.border.bottom.color = "black",

    column_labels.border.bottom.width = px(2),

    table_body.hlines.color = "gray",

    table_body.border.bottom.color = "black",

    row.striping.include_table_body = TRUE
  ) %>%

  tab_style(

    style = cell_text(weight = "bold"),

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

tabla_rasgo_gt
Tabla Nº1
Distribución ordinal del tipo de rasgo geológico
Tipo_Rasgo_Geologico ni hi P
Rasgo menor 987 0.3948 39.48
Rasgo intermedio 1130 0.4520 45.20
Rasgo mayor 383 0.1532 15.32
Otros 0 0.0000 0.00
TOTAL 2500 1.0000 100.00
Autor: Grupo 2 
# =========================================================
# GRÁFICA DESCRIPTIVA
# =========================================================

P_global <- tabla_finalrasgo$P[
  tabla_finalrasgo$Tipo_Rasgo_Geologico != "TOTAL"
]

barplot(

  P_global,

  main = "Gráfica Nº1: Distribución ordinal del tipo\nde rasgo geológico",

  xlab = "Tipo de rasgo geológico",

  ylab = "Probabilidad (%)",

  col = "blue",

  names.arg = tabla_finalrasgo$Tipo_Rasgo_Geologico[
    tabla_finalrasgo$Tipo_Rasgo_Geologico != "TOTAL"
  ],

  las = 2,

  ylim = c(0,100)
)

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

n <- sum(
  tabla_finalrasgo$ni[
    tabla_finalrasgo$Tipo_Rasgo_Geologico != "TOTAL"
  ]
)

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

x <- tabla_finalrasgo$ni[
  tabla_finalrasgo$Tipo_Rasgo_Geologico != "TOTAL"
]

x
## [1]  987 1130  383    0
# =========================================================
# ESCALA ORDINAL
# =========================================================

X <- 1:length(x)

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

media_observada <- sum(X * x) / n

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

p <- media_observada / length(x)

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

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

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

P_binomial
## [1] 0.30946537 0.36413537 0.19042816 0.03734485
# =========================================================
# FRECUENCIAS OBSERVADAS Y ESPERADAS
# =========================================================

Fo <- (x / n) * 100

Fo
## [1] 39.48 45.20 15.32  0.00
Fe <- P_binomial * 100

Fe
## [1] 30.946537 36.413537 19.042816  3.734485
barplot(

  rbind(Fo, Fe),

  beside = TRUE,

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

  names.arg = tabla_finalrasgo$Tipo_Rasgo_Geologico[
    tabla_finalrasgo$Tipo_Rasgo_Geologico != "TOTAL"
  ],

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

  ylab = "Probabilidad (%)",

  xlab = "Tipo de rasgo geológico",

  ylim = c(0,100),

  las = 2
)

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] 98.81137
K <- length(x)

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

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

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

vc
## [1] 11.34487
x2 < vc
## [1] TRUE
# =========================================================
# ELIMINAR TOTAL
# =========================================================

tabla_sin_total <- tabla_finalrasgo[
  tabla_finalrasgo$Tipo_Rasgo_Geologico != "TOTAL",
]

# =========================================================
# PROBABILIDAD DE RASGO MAYOR
# =========================================================

prob_mayor <- tabla_sin_total$P[
  tabla_sin_total$Tipo_Rasgo_Geologico == "Rasgo mayor"
]

prob_mayor
## [1] 15.32
# =========================================================
# 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 rasgo\n",
    "geológico analizado corresponda a un\n",
    "RASGO MAYOR?\n\n",

    "Probabilidad = ",
    prob_mayor,
    " (%)",

    sep = ""
  ),

  cex = 1.2,

  col = "black",

  font = 2
)

Variable <- c("Tipo de rasgo geológico")

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
Tipo de rasgo geológico 98.81 8.94 11.34