1 Configuración y Carga de Datos

library(dplyr)
library(gt)
library(e1071)

col_principal <- "#0E6655"
col_barras    <- "#16A085"
col_acento    <- "#E67E22"
col_grid      <- "#D7DBDD"

setwd("C:/Users/ASUS/Desktop/Estadistica/new_york_exel")

archivo_csv <- "Oil__Gas____Other_Regulated_Wells__Beginning_1860.csv"
# Si tu archivo se llama distinto (p. ej. con sufijo "__1_"), cambia la línea de arriba

# Detección de separador (evita depender de un delimitador fijo)
separadores <- c(",", ";", "\t", "|")
mejor_sep <- NULL
mejor_ncol <- 1

for (s in separadores) {
  n_campos <- tryCatch(utils::count.fields(archivo_csv, sep = s)[1],
                        error = function(e) 1)
  if (!is.na(n_campos) && n_campos > mejor_ncol) {
    mejor_ncol <- n_campos
    mejor_sep  <- s
  }
}
if (is.null(mejor_sep)) mejor_sep <- ","
cat("Separador detectado:", ifelse(mejor_sep == "\t", "TAB", mejor_sep),
    "| Columnas detectadas:", mejor_ncol, "\n")
## Separador detectado: , | Columnas detectadas: 52
Datos_Brutos <- read.csv(archivo_csv, header = TRUE, sep = mejor_sep,
                          check.names = TRUE, stringsAsFactors = FALSE)

# Localización de la columna de fecha de emisión del permiso
col_permiso <- names(Datos_Brutos)[
  grepl("permit", names(Datos_Brutos), ignore.case = TRUE) &
  grepl("issued", names(Datos_Brutos), ignore.case = TRUE) &
  grepl("date",   names(Datos_Brutos), ignore.case = TRUE)
]

if (length(col_permiso) == 0) {
  columnas_txt <- paste(names(Datos_Brutos), collapse = " | ")
  stop("ERROR: No se encontró ninguna columna de fecha de emisión del permiso ",
       "(se buscó el patrón 'Permit'+'Issued'+'Date').\n",
       "Separador usado: '", mejor_sep, "'  |  N° de columnas leídas: ", ncol(Datos_Brutos), "\n",
       "COLUMNAS ENCONTRADAS EN EL ARCHIVO:\n", columnas_txt)
}

nombre_col_permiso <- col_permiso[1]
cat("Columna de fecha de emisión del permiso identificada:", nombre_col_permiso, "\n")
## Columna de fecha de emisión del permiso identificada: Permit.Issued.Date

2 Extracción y Limpieza de la Variable

# La fecha viene en formato MM/DD/AAAA; se extrae el año con expresión regular
Datos <- Datos_Brutos %>%
  mutate(Anio_Permiso = suppressWarnings(as.integer(
    sub(".*/([0-9]{4}).*", "\\1", .data[[nombre_col_permiso]])
  )))

Datos <- Datos %>%
  filter(!is.na(Anio_Permiso) & Anio_Permiso >= 1900 & Anio_Permiso <= 2026) %>%
  mutate(
    Decada_Permiso = floor(Anio_Permiso / 10) * 10,
    Periodo = paste0(Decada_Permiso, " - ", Decada_Permiso + 9)
  )

Variable_Exacta <- Datos$Anio_Permiso
if (length(Variable_Exacta) == 0) stop("ERROR: No hay datos válidos.")

n_sin_dato <- nrow(Datos_Brutos) - length(Variable_Exacta)

cat("N (pozos con fecha de emisión de permiso 1900-2026):", length(Variable_Exacta),
    "| Décadas:", length(unique(Datos$Decada_Permiso)), "\n")
## N (pozos con fecha de emisión de permiso 1900-2026): 18286 | Décadas: 8
cat("Registros sin fecha de emisión de permiso o fuera de rango:", n_sin_dato,
    "(", round(100 * n_sin_dato / nrow(Datos_Brutos), 1), "%)\n")
## Registros sin fecha de emisión de permiso o fuera de rango: 29104 ( 61.4 %)

3 Identificación de la Variable

cat("Variable      : Fecha de Emisión del Permiso (año de emisión del permiso de perforación)\n")
## Variable      : Fecha de Emisión del Permiso (año de emisión del permiso de perforación)
cat("Columna base  :", nombre_col_permiso, "\n")
## Columna base  : Permit.Issued.Date
cat("Tipo          : Cuantitativa continua (escala de intervalo; agrupada en clases de 10 años para las gráficas)\n")
## Tipo          : Cuantitativa continua (escala de intervalo; agrupada en clases de 10 años para las gráficas)
cat("Rango         :", min(Variable_Exacta), "-", max(Variable_Exacta), "\n")
## Rango         : 1953 - 2026
cat("N             :", length(Variable_Exacta), "\n")
## N             : 18286

4 Tabla de Distribución de Frecuencias

TDF_Raw <- Datos %>%
  group_by(Decada_Permiso, Periodo) %>%
  summarise(ni = n(), .groups = "drop") %>%
  arrange(Decada_Permiso)

ni <- TDF_Raw$ni
N  <- sum(ni)
hi <- (ni / N) * 100
Ni_asc  <- cumsum(ni)
Ni_desc <- rev(cumsum(rev(ni)))
Hi_asc  <- cumsum(hi)
Hi_desc <- rev(cumsum(rev(hi)))

TDF_Decadas <- data.frame(
  Periodo = TDF_Raw$Periodo,
  ni = ni,
  hi = round(hi, 2),
  Ni_asc = Ni_asc,
  Ni_desc = Ni_desc,
  Hi_asc = round(Hi_asc, 2),
  Hi_desc = round(Hi_desc, 2)
)

totales <- c("TOTAL", sum(ni), round(sum(hi), 2), "-", "-", "-", "-")
TDF_Char <- TDF_Decadas %>% mutate(across(everything(), as.character))
TDF_Show <- rbind(TDF_Char, totales)

modal_row <- which.max(TDF_Decadas$ni)

TDF_Show %>%
  gt() %>%
  tab_header(
    title = md("**DISTRIBUCIÓN POR DÉCADAS**"),
    subtitle = md(paste0("Variable: **Fecha de Emisión del Permiso (",
                          nombre_col_permiso, ")** · Nueva York"))
  ) %>%
  tab_spanner(label = "Frecuencias acumuladas",
              columns = c(Ni_asc, Ni_desc, Hi_asc, Hi_desc)) %>%
  cols_label(
    Periodo = "Periodo (Década)",
    ni = "Cant. Pozos (ni)",
    hi = "Porcentaje (hi%)",
    Ni_asc = "Ni (Asc)", Ni_desc = "Ni (Desc)",
    Hi_asc = "Hi (Asc)", Hi_desc = "Hi (Desc)"
  ) %>%
  cols_align(align = "center", columns = everything()) %>%
  tab_style(
    style = list(cell_fill(color = col_principal), cell_text(color = "white", weight = "bold")),
    locations = cells_title()
  ) %>%
  tab_style(
    style = list(cell_fill(color = "#148F77"), cell_text(color = "white", weight = "bold")),
    locations = cells_column_labels()
  ) %>%
  tab_style(
    style = list(cell_fill(color = "#148F77"), cell_text(color = "white", weight = "bold")),
    locations = cells_column_spanners()
  ) %>%
  tab_style(
    style = list(cell_fill(color = "#FDEBD0"), cell_text(weight = "bold")),
    locations = cells_body(rows = modal_row)
  ) %>%
  tab_style(
    style = list(cell_fill(color = "#D0ECE7"), cell_text(weight = "bold")),
    locations = cells_body(rows = Periodo == "TOTAL")
  ) %>%
  opt_row_striping() %>%
  opt_table_font(font = google_font("Roboto")) %>%
  tab_options(
    table.font.size = px(13),
    heading.align = "left",
    heading.title.font.size = px(17),
    data_row.padding = px(7),
    table.border.top.color = col_principal,
    table.border.bottom.color = col_principal,
    column_labels.border.bottom.color = col_principal
  ) %>%
  tab_source_note(md("*Fuente: NYS DEC — Oil, Gas & Other Regulated Wells. Elaboración: Dallyanna Lozano (Grupo 1).*"))
DISTRIBUCIÓN POR DÉCADAS
Variable: Fecha de Emisión del Permiso (Permit.Issued.Date) · Nueva York
Periodo (Década) Cant. Pozos (ni) Porcentaje (hi%)
Frecuencias acumuladas
Ni (Asc) Ni (Desc) Hi (Asc) Hi (Desc)
1950 - 1959 7 0.04 7 18286 0.04 100
1960 - 1969 566 3.1 573 18279 3.13 99.96
1970 - 1979 4057 22.19 4630 17713 25.32 96.87
1980 - 1989 6213 33.98 10843 13656 59.3 74.68
1990 - 1999 1577 8.62 12420 7443 67.92 40.7
2000 - 2009 3750 20.51 16170 5866 88.43 32.08
2010 - 2019 1691 9.25 17861 2116 97.68 11.57
2020 - 2029 425 2.32 18286 425 100 2.32
TOTAL 18286 100 - - - -
Fuente: NYS DEC — Oil, Gas & Other Regulated Wells. Elaboración: Dallyanna Lozano (Grupo 1).

5 Representación Gráfica

5.1 Gráfica N°1 — Histograma (Frecuencia Absoluta)

vals_x <- TDF_Decadas$Periodo
vals_y <- TDF_Decadas$ni
ylim_max <- max(vals_y) * 1.15

par(mar = c(10, 5, 4, 2))
bp1 <- barplot(vals_y,
        main = "Gráfica N°1: Histograma de la fecha de emisión del permiso por década (NY)",
        cex.main = 0.9, ylab = "Cantidad de pozos (ni)",
        col = col_barras, border = "white", axes = FALSE, space = 0,
        ylim = c(0, ylim_max), axisnames = FALSE)
axis(2, col = col_principal, col.axis = col_principal)
axis(1, at = bp1, labels = vals_x, col = col_principal, col.axis = col_principal, las = 2, cex.axis = 0.8)
text(x = bp1, y = vals_y, label = vals_y, pos = 3, cex = 0.7, col = col_principal)
title(xlab = "Década de la emisión", line = 8)
grid(nx = NA, ny = NULL, col = col_grid, lty = "dotted")
box(bty = "l", col = col_principal)

5.2 Gráfica N°2 — Histograma (Porcentaje)

vals_y_pct <- TDF_Decadas$hi

par(mar = c(10, 5, 4, 2))
bp2 <- barplot(vals_y_pct,
        main = "Gráfica N°2: Histograma porcentual de la fecha de emisión del permiso (NY)",
        cex.main = 0.9, ylab = "% del total (hi)",
        col = col_barras, border = "white", axes = FALSE, space = 0,
        ylim = c(0, max(vals_y_pct) * 1.2), axisnames = FALSE)
axis(2, col = col_principal, col.axis = col_principal)
axis(1, at = bp2, labels = vals_x, col = col_principal, col.axis = col_principal, las = 2, cex.axis = 0.8)
text(x = bp2, y = vals_y_pct, label = paste0(round(vals_y_pct, 1), "%"), pos = 3, cex = 0.7, col = col_principal)
title(xlab = "Década de la emisión", line = 8)
grid(nx = NA, ny = NULL, col = col_grid, lty = "dotted")
box(bty = "l", col = col_principal)

5.3 Gráfica N°3 — Histograma de Clases Finas (Año Exacto)

par(mar = c(6, 5, 4, 2))
h3 <- hist(Variable_Exacta, plot = FALSE)

hp3 <- barplot(h3$counts,
        main = "Gráfica N°3: Histograma de clases finas de la fecha de emisión del permiso (NY)",
        cex.main = 0.9, ylab = "Cantidad de pozos (ni)",
        col = col_acento, border = "white", axes = FALSE, space = 0,
        ylim = c(0, max(h3$counts) * 1.15))
axis(2, col = col_principal, col.axis = col_principal)
etiquetas_clase <- paste0(h3$breaks[-length(h3$breaks)], "-", h3$breaks[-1])
axis(1, at = hp3, labels = etiquetas_clase, col = col_principal, col.axis = col_principal,
     las = 2, cex.axis = 0.7)
text(x = hp3, y = h3$counts, label = h3$counts, pos = 3, cex = 0.7, col = col_principal)
title(xlab = "Clase (años)", line = 4.5)
grid(nx = NA, ny = NULL, col = col_grid, lty = "dotted")
box(bty = "l", col = col_principal)

5.4 Diagrama de Caja (Año Exacto de la Fecha de Emisión del Permiso)

Muestra dispersión y valores atípicos del año exacto (previo a agrupar en décadas).

par(mar = c(6, 5, 4, 2))
boxplot(Variable_Exacta, horizontal = TRUE, col = col_barras,
        main = "Diagrama de caja del año de la fecha de emisión del permiso (NY)",
        cex.main = 0.9, xlab = "", outline = TRUE, outpch = 19, outcol = col_acento,
        boxwex = 0.5, frame.plot = FALSE, xaxt = "n")
eje_x <- pretty(Variable_Exacta, n = 8)
axis(1, at = eje_x, labels = format(eje_x, scientific = FALSE), cex.axis = 0.8,
     col = col_principal, col.axis = col_principal)
title(xlab = "Año exacto de la fecha de emisión del permiso", line = 4)

5.5 Ojivas (Frecuencias Acumuladas por Década)

Muestra la frecuencia acumulada ascendente y descendente por década.

par(mar = c(10, 5, 4, 8), xpd = TRUE)
x_vals_num <- sort(unique(Datos$Decada_Permiso))
x_labels   <- TDF_Decadas$Periodo
y_asc  <- TDF_Decadas$Ni_asc
y_desc <- TDF_Decadas$Ni_desc

ylim_ojivas <- c(0, max(y_asc, y_desc) * 1.1)

plot(x_vals_num, y_asc, type = "o", col = col_principal, lwd = 2, pch = 19,
     main = "Ojivas ascendente y descendente de la fecha de emisión del permiso (NY)",
     cex.main = 0.9, ylab = "Frecuencia acumulada", xlab = "",
     axes = FALSE, frame.plot = FALSE, ylim = ylim_ojivas)
axis(1, at = x_vals_num, labels = x_labels, las = 2, cex.axis = 0.8,
     col = col_principal, col.axis = col_principal)
axis(2, col = col_principal, col.axis = col_principal)
title(xlab = "Década de la emisión", line = 8)
lines(x_vals_num, y_desc, type = "o", col = col_acento, lwd = 2, pch = 19)
legend("right", legend = c("Ascendente", "Descendente"),
       col = c(col_principal, col_acento), lty = 1, pch = 19, cex = 0.7, lwd = 2,
       inset = c(-0.15, 0), bty = "n")
grid(nx = NULL, ny = NULL, col = col_grid, lty = "dotted")

6 Tabla de Indicadores

media_val   <- mean(Variable_Exacta)
mediana_val <- median(Variable_Exacta)

t_moda     <- table(Variable_Exacta)
freq_max   <- max(t_moda)
modas_calc <- as.numeric(names(t_moda)[t_moda == freq_max])
moda_txt   <- paste(modas_calc, collapse = ", ")

rango_txt    <- paste0("[", min(Variable_Exacta), "; ", max(Variable_Exacta), "]")
varianza_val <- var(Variable_Exacta)
sd_val       <- sd(Variable_Exacta)
cv_val       <- (sd_val / abs(media_val)) * 100
asimetria_val <- skewness(Variable_Exacta, type = 2)
curtosis_val  <- kurtosis(Variable_Exacta, type = 2)

vals_atipicos <- boxplot.stats(Variable_Exacta)$out
num_atipicos  <- length(vals_atipicos)
status_atipicos <- if (num_atipicos > 0) paste0(num_atipicos, " outliers") else "0 (Sin atípicos)"

df_indicadores <- data.frame(
  "Variable"  = "Fecha de emisión del permiso (año exacto base)",
  "Rango"     = rango_txt,
  "Media"     = media_val,
  "Mediana"   = mediana_val,
  "Moda"      = moda_txt,
  "Varianza"  = varianza_val,
  "Desv_Std"  = sd_val,
  "CV_Porc"   = cv_val,
  "Asimetria" = asimetria_val,
  "Curtosis"  = curtosis_val,
  "Atipicos"  = status_atipicos
)

df_indicadores %>%
  gt() %>%
  tab_header(
    title = md("**TABLA DE INDICADORES**"),
    subtitle = "Estadísticos descriptivos de la variable Fecha de Emisión del Permiso"
  ) %>%
  tab_source_note(source_note = "Autor: EDUARDO") %>%
  fmt_number(columns = c(Media, Mediana, Varianza, Desv_Std, CV_Porc, Curtosis), decimals = 2) %>%
  fmt_number(columns = c(Asimetria), decimals = 4) %>%
  cols_label(
    Variable = "Variable", Rango = "Rango [Min; Max]",
    Media = "Media (X̄)", Mediana = "Mediana (Me)", Moda = "Moda (Mo)",
    Varianza = "Varianza (S²)", Desv_Std = "Desv. Est. (S)", CV_Porc = "C.V. (%)",
    Asimetria = "Asimetría (As)", Curtosis = "Curtosis (K)", Atipicos = "Outliers"
  ) %>%
  cols_align(align = "center", columns = everything()) %>%
  tab_style(
    style = list(cell_fill(color = col_principal), cell_text(color = "white", weight = "bold")),
    locations = cells_title()
  ) %>%
  tab_style(
    style = list(cell_fill(color = "#148F77"), cell_text(color = "white", weight = "bold")),
    locations = cells_column_labels()
  ) %>%
  tab_style(
    style = list(cell_fill(color = "#FDEBD0"), cell_text(weight = "bold")),
    locations = cells_body(columns = Moda)
  ) %>%
  opt_table_font(font = google_font("Roboto")) %>%
  tab_options(
    table.font.size = px(13),
    heading.align = "left",
    data_row.padding = px(9),
    table.border.top.color = col_principal,
    table.border.bottom.color = col_principal,
    column_labels.border.bottom.color = col_principal
  )
TABLA DE INDICADORES
Estadísticos descriptivos de la variable Fecha de Emisión del Permiso
Variable Rango [Min; Max] Media (X̄) Mediana (Me) Moda (Mo) Varianza (S²) Desv. Est. (S) C.V. (%) Asimetría (As) Curtosis (K) Outliers
Fecha de emisión del permiso (año exacto base) [1953; 2026] 1,990.33 1,985.00 1981 210.11 14.50 0.73 0.4850 −0.98 0 (Sin atípicos)
Autor: EDUARDO

7 Conclusión

El periodo analizado va de 1953 a 2026 (N = 18286 pozos con fecha de emisión de permiso registrada; 29104 pozos no tienen este dato). La década con más registros de emisión de permiso fue 1980 - 1989 (moda: 1981). CV = 0.73% (ritmo relativamente constante). Asimetría = 0.485 (un sesgo hacia los años antiguos (asimetría positiva)). Sin años atípicos.