Enlace de acceso al informe HTML publicado: https://rpubs.com/DavidCapetillo/Actividad3_GRD_H25

1 Resultado 01: Una variable aleatoria discreta construida desde los datos

1.1 Correcciones aplicadas a partir de las Actividades 1 y 2

En cumplimiento de la retroalimentación recibida en las Evaluaciones de las Actividades 1 y 2, se implementaron y verificaron sobre el flujo reproducible de datos las siguientes correcciones metodológicas y de comunicación antes de estimar los modelos probabilísticos de esta actividad:

Eliminación explícita de registros duplicados exactos: En la Actividad 1 se habían detectado $330$ filas idénticas en la base bruta sin declarar una regla de decisión sobre ellas. En esta entrega se aplica el filtro !duplicated() sobre las $129$ columnas de la base original de $1.085.813$ registros, descartando las $330$ filas duplicadas ($329$ pertenecientes a Cirugía Mayor Ambulatoria y $1$ a hospitalización general sin código H25), dejando $1.085.483$ registros únicos.
Trazabilidad exacta de exclusiones por calidad de datos sobre $\Omega$: Tras excluir los $211.712$ registros de Cirugía Mayor Ambulatoria (CMA), quedan $873.771$ egresos de atención cerrada. Sobre este conjunto se excluyen los $53$ episodios con fechas faltantes, estancias negativas ($< 0$ días: $52$ casos) o edades fuera del rango biológico compatible ($< 0$ o $> 110$ años: $53$ casos, $52$ de los cuales coinciden con fechas inválidas). El espacio muestral depurado definitivo es $\Omega = 873.718$ episodios hospitalarios.
Revisión de edades menores de 40 años en la cohorte H25: Dentro de los $389$ episodios que registran el código CIE-10 H25 (Catarata senil), se auditó la presencia de $9$ episodios con edades entre $15$ y $39$ años ($6$ en rol diagnóstico secundario asociados a diabetes mellitus E11.9, desprendimiento de retina H33.0 o glaucoma H40.5, y $3$ en rol principal codificados como H25.9 en el Hospital Base San José de Osorno). Se documenta la decisión de mantener los $389$ episodios publicados en la base oficial de FONASA para preservar la trazabilidad exacta con los parámetros de la Actividad 2 ($n_E = 389$, con $16$ fallecimientos y $p = 16/389 = 4,1131\%$), dejando constancia de que estos $9$ registros reflejan imprecisiones de codificación administrativa donde se empleó H25.9 en lugar de H26 (otras cataratas).
Cruce oficial de nombres de establecimientos con la tabla maestra: Se eliminó la asignación manual de nombres de hospitales observada en la Actividad 1 y se realizó un cruce relacional (left_join) por COD_HOSPITAL contra la hoja Hospitales del archivo auxiliar oficial Tablas_maestras_bases_GRD.xlsx. Esto corrige la identificación del código 118100 como Hospital Clínico Regional Dr. Guillermo Grant Benavente (Concepción) ($83$ episodios) y del código 116105 como Hospital Dr. César Garavagno Burotto (Talca) ($22$ episodios).
Ajuste del marcador de registro $M$ y precisión terminológica: Dado que el marcador original de la Actividad 2 ($\text{edad} \ge 65 \cap \overline{U}$) presentó un valor predictivo positivo casi nulo ($\text{VPP}_{\text{orig}} = 110 / 57.303 = 0,1920\%$), se ajustó $M$ incorporando la especialidad o servicio clínico de Oftalmología en episodios de $\ge 65$ años (criterio administrativo no circular que no utiliza campos diagnósticos), elevando el valor predictivo positivo a $\text{VPP} = 99 / 962 = 10,2911\%$, sin perjuicio de reportar comparativamente ambos marcadores en los Resultados 04 y 06. Asimismo, todo el informe se refiere estrictamente a episodios de egreso hospitalario (nunca a personas, prevalencia, incidencia, riesgo individual ni métodos MCMC).

ruta_base <- "GRD_PUBLICO_2024.txt"
ruta_maestra <- "Tablas_maestras_bases_GRD.xlsx"

# 1. Lectura de la base GRD 2024 con ruta relativa
grd_base <- fread(ruta_base, sep = "|", encoding = "Latin-1")
n_01_bruta <- nrow(grd_base)

# 2. Eliminación de duplicados exactos en las 129 columnas
dup_idx <- duplicated(grd_base)
n_dup <- sum(dup_idx)

cols_diag <- paste0("DIAGNOSTICO", 1:35)
cols_keep <- c(
  "COD_HOSPITAL", "SEXO", "FECHA_NACIMIENTO", "SERVICIO_SALUD",
  "TIPO_INGRESO", "ESPECIALIDAD_MEDICA", "TIPO_ACTIVIDAD",
  "FECHA_INGRESO", "SERVICIOINGRESO", "FECHAALTA", "SERVICIOALTA",
  "TIPOALTA", "IR_29301_PESO", "IR_29301_SEVERIDAD", "IR_29301_MORTALIDAD",
  cols_diag
)

grd_sin_dup <- grd_base[!dup_idx, ..cols_keep]
rm(grd_base)
n_02_sin_dup <- nrow(grd_sin_dup)

# 3. Exclusión de Cirugía Mayor Ambulatoria (CMA)
es_cma_vec <- grd_sin_dup$TIPO_ACTIVIDAD != "HOSPITALIZACIÓN"
n_cma <- sum(es_cma_vec)
grd_egresos <- grd_sin_dup[!es_cma_vec]
rm(grd_sin_dup)
n_03_egresos <- nrow(grd_egresos)

# Conversión vectorizada de fechas mediante tabla de valores únicos
fast_date <- function(x) {
  ux <- unique(x)
  ud <- as.Date(ux, format = "%Y-%m-%d")
  ud[match(x, ux)]
}

grd_egresos[, fecha_ing := fast_date(FECHA_INGRESO)]
grd_egresos[, f_alta := fast_date(FECHAALTA)]
grd_egresos[, f_nac := fast_date(FECHA_NACIMIENTO)]
grd_egresos[, estancia := as.numeric(f_alta - fecha_ing)]
grd_egresos[, edad := floor(as.numeric(fecha_ing - f_nac) / 365.25)]

# 4. Filtro de consistencia lógica de estancia y edad
cond_calidad <- !is.na(grd_egresos$estancia) & grd_egresos$estancia >= 0 &
                !is.na(grd_egresos$edad) & grd_egresos$edad >= 0 & grd_egresos$edad <= 110
n_excl_calidad <- sum(!cond_calidad)
omega <- grd_egresos[which(cond_calidad)]
rm(grd_egresos)
N_omega <- nrow(omega)

# 5. Identificación de la cohorte E (H25) y eventos basales sobre Omega
es_h25_mat <- sapply(cols_diag, function(col) substr(omega[[col]], 1, 3) == "H25")
omega[, E_p := es_h25_mat[, 1]]
omega[, E := rowSums(es_h25_mat) > 0]
omega[, E_s := E & !E_p]
rm(es_h25_mat)

corte_L <- as.numeric(quantile(omega$estancia, 0.75, na.rm = TRUE))
fast_grep <- function(x, pat) {
  ux <- unique(x)
  ug <- grepl(pat, ux, ignore.case = TRUE)
  ug[match(x, ux)]
}
omega[, F_evt := fast_grep(TIPOALTA, "FALLE|MUERTE|3")]
omega[, U_evt := fast_grep(TIPO_INGRESO, "URGENCIA|1")]
omega[, L_evt := estancia > corte_L]

# Marcador original (Actividad 2) y Marcador ajustado (Actividad 3)
omega[, M_orig := (edad >= 65 & !U_evt)]
omega[, es_oft := fast_grep(ESPECIALIDAD_MEDICA, "OFTALM") | fast_grep(SERVICIOALTA, "OFTALM")]
omega[, M := (es_oft & edad >= 65)]

# Subconjunto de la cohorte E y conteo de diagnósticos no vacíos (Variable X)
cohorte_E <- omega[which(omega$E)]
no_vacio_mat <- sapply(cols_diag, function(col) {
  v <- trimws(cohorte_E[[col]])
  !is.na(v) & v != ""
})
cohorte_E[, X := rowSums(no_vacio_mat)]

# Cruce con tabla maestra oficial de establecimientos
hosp_maestra <- read_excel(ruta_maestra, sheet = "Hospitales")
colnames(hosp_maestra) <- c("COD_HOSPITAL", "NOMBRE_HOSPITAL")
hosp_maestra <- hosp_maestra %>% mutate(COD_HOSPITAL = as.numeric(COD_HOSPITAL))

cohorte_E <- merge(cohorte_E, as.data.table(hosp_maestra), by = "COD_HOSPITAL", all.x = TRUE)
cohorte_E[is.na(NOMBRE_HOSPITAL), NOMBRE_HOSPITAL := paste0("Hospital Cód. ", COD_HOSPITAL)]

n_E <- nrow(cohorte_E)
n_Ep <- sum(cohorte_E$E_p)
n_Es <- sum(cohorte_E$E_s)
P_E <- n_E / N_omega

tabla_correcciones <- tibble(
  Etapa = c(
    "1. Base bruta GRD 2024 cargada desde archivo fuente",
    "2. Eliminación de filas duplicadas exactas (!duplicated)",
    "3. Exclusión de Cirugía Mayor Ambulatoria (CMA)",
    "4. Exclusión por estancia (< 0 o NA) o edad incompatible (< 0 o > 110) [Espacio muestral Ω]",
    "5. Cohorte E: Episodios con Catarata Senil (CIE-10: H25 en cualquier diagnóstico)",
    "   5a. Subcohorte Ep: H25 registrado como diagnóstico principal (DIAGNOSTICO1)",
    "   5b. Subcohorte Es: H25 registrado solo como diagnóstico secundario (DIAGNOSTICO2..35)"
  ),
  Cambio_Registros = c(
    format(n_01_bruta, big.mark = "."),
    paste0("-", format(n_dup, big.mark = ".")),
    paste0("-", format(n_cma, big.mark = ".")),
    paste0("-", format(n_excl_calidad, big.mark = ".")),
    format(n_E, big.mark = "."),
    format(n_Ep, big.mark = "."),
    format(n_Es, big.mark = ".")
  ),
  Total_Acumulado = c(
    format(n_01_bruta, big.mark = "."),
    format(n_02_sin_dup, big.mark = "."),
    format(n_03_egresos, big.mark = "."),
    format(N_omega, big.mark = "."),
    paste0(n_E, " (", sprintf("%.4f%%", P_E * 100), " de Ω)"),
    paste0(n_Ep, " (", sprintf("%.2f%%", (n_Ep / n_E) * 100), " de E)"),
    paste0(n_Es, " (", sprintf("%.2f%%", (n_Es / n_E) * 100), " de E)")
  )
)

kable(
  tabla_correcciones,
  col.names = c("Etapa del flujo de depuración y correcciones aplicadas", "Variación", "Tamaño resultante y proporción"),
  align = c("l", "r", "r"),
  caption = "Tabla 1: Trazabilidad del espacio muestral Ω y tamaño final de la cohorte E tras aplicar las correcciones de las Actividades 1 y 2."
) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)

Tabla 1: Trazabilidad del espacio muestral Ω y tamaño final de la cohorte E tras aplicar las correcciones de las Actividades 1 y 2.
Etapa del flujo de depuración y correcciones aplicadas	Variación	Tamaño resultante y proporción
Base bruta GRD 2024 cargada desde archivo fuente	1.085.813	1.085.813
Eliminación de filas duplicadas exactas (!duplicated)	-330	1.085.483
Exclusión de Cirugía Mayor Ambulatoria (CMA)	-211.712	873.771
Exclusión por estancia (< 0 o NA) o edad incompatible (< 0 o > 110) [Espacio muestral Ω]	-53	873.718
Cohorte E: Episodios con Catarata Senil (CIE-10: H25 en cualquier diagnóstico)	389	389 (0.0445% de Ω)
5a. Subcohorte Ep: H25 registrado como diagnóstico principal (DIAGNOSTICO1)	108	108 (27.76% de E)
5b. Subcohorte Es: H25 registrado solo como diagnóstico secundario (DIAGNOSTICO2..35)	281	281 (72.24% de E)

Interpretación de la Tabla 1: La depuración secuencial sobre los $1.085.813$ registros iniciales descarta $330$ filas duplicadas exactas, $211.712$ atenciones ambulatorias (CMA) y $53$ registros con fechas o edades inconsistentes, consolidando un espacio muestral de $N = 873.718$ episodios hospitalarios ($\Omega$). Sobre este denominador, la cohorte final de Catarata Senil ($E$) mantiene exactamente $n_E = 389$ episodios ($P(E) = 389 / 873.718 = 0,0445\%$), de los cuales $108$ ($27,76\%$) corresponden a ingresos donde H25 es el diagnóstico principal ($E_p$) y $281$ ($72,24\%$) a episodios donde actúa como comorbilidad secundaria ($E_s$).

1.2 Definición de variables aleatorias discretas, soporte y modelos propuestos

Cada variable aleatoria asigna un valor numérico discreto al resultado de un experimento bien definido sobre los episodios hospitalarios de la cohorte $E$ o del espacio muestral $\Omega$. En la Tabla 2 se documentan todas las variables exigidas (incluyendo las variantes sobre estancia larga $L$ para complementar los análisis donde el número de fallecimientos es reducido), especificando su experimento generador, su soporte matemático, el modelo de probabilidad propuesto, el origen de sus parámetros y su implementación en R.

tabla_vars <- tibble(
  Variable = c("X", "Y", "Y_L", "S", "S_L", "H", "T", "T_5", "D"),
  Experimento = c(
    "Elegir 1 episodio al azar de la cohorte E (n_E = 389) y contar sus campos de diagnóstico (DIAGNOSTICO1 a DIAGNOSTICO35) no vacíos.",
    "Elegir 1 episodio al azar de la cohorte E y observar si termina en fallecimiento (1 si F ocurre, 0 si no).",
    "Elegir 1 episodio al azar de la cohorte E y observar si registra estancia larga > 7 días (1 si L ocurre, 0 si no).",
    "Seleccionar 20 episodios al azar con reposición de la cohorte E y contar cuántos terminan en fallecimiento.",
    "Seleccionar 20 episodios al azar con reposición de la cohorte E y contar cuántos presentan estancia larga (> 7 días).",
    "Auditar n = 10 episodios al azar sin reposición dentro de un establecimiento h (con N_h episodios de E) y contar los casos con el evento de interés (K_h).",
    "Revisar uno por uno episodios marcados por M (con reposición) hasta encontrar el primer episodio perteneciente a la cohorte E.",
    "Revisar uno por uno episodios marcados por M (con reposición) hasta acumular r = 5 episodios pertenecientes a la cohorte E.",
    "Elegir un día calendario del año bisiesto 2024 (366 días) y contar cuántos episodios de la cohorte E egresaron en esa fecha."
  ),
  Soporte = c(
    "{1, 2, ..., 35}",
    "{0, 1}",
    "{0, 1}",
    "{0, 1, ..., 20}",
    "{0, 1, ..., 20}",
    "{max(0, 10-(N_h-K_h)), ..., min(10, K_h)}",
    "{1, 2, 3, ...}",
    "{5, 6, 7, ...}",
    "{0, 1, 2, ...}"
  ),
  Modelo_Propuesto = c(
    "Distribución empírica p(x) = n(X=x) / 389",
    "Bernoulli(p), con p = P(F | E) = 16 / 389 = 0,041131",
    "Bernoulli(p_L), con p_L = P(L | E) = 121 / 389 = 0,311054",
    "Binomial(n = 20, p = 0,041131)",
    "Binomial(n = 20, p_L = 0,311054)",
    "Hipergeométrica(N_h, K_h, n = 10)",
    "Geométrica(p = VPP), con VPP = P(E | M) = 99 / 962 = 0,102911 (y M_orig: 110 / 57.303 = 0,001920)",
    "Binomial Negativa(r = 5, p = VPP)",
    "Poisson(λ), con λ = 389 / 366 = 1,062842 egresos/día"
  ),
  Calculo_en_R = c(
    "rowSums(!is.na(DIAG) & DIAG != '') ; table(X)/nrow(cohorte_E)",
    "as.integer(F_evt) ; mean(F_evt)",
    "as.integer(L_evt) ; mean(L_evt)",
    "dbinom(k, 20, p), pbinom(k, 20, p), qbinom(q, 20, p), rbinom(B, 20, p)",
    "dbinom(k, 20, p_L), pbinom(k, 20, p_L), qbinom(q, 20, p_L)",
    "dhyper(k, m = K_h, n = N_h - K_h, k = 10)",
    "dgeom(k - 1, prob = VPP), pgeom(k - 1, prob = VPP), qgeom(q, prob = VPP) + 1",
    "dnbinom(k - 5, size = 5, prob = VPP), pnbinom(k - 5, size = 5, prob = VPP)",
    "dpois(k, lambda), ppois(k, lambda), rpois(B, lambda)"
  )
)

kable(
  tabla_vars,
  col.names = c("Variable", "Experimento aleatorio (unidad: episodio)", "Soporte", "Modelo y origen del parámetro", "Funciones y cálculo en R"),
  caption = "Tabla 2: Definición operacional de las variables aleatorias discretas, experimento generador, soporte y modelo probabilístico."
) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = TRUE)

Tabla 2: Definición operacional de las variables aleatorias discretas, experimento generador, soporte y modelo probabilístico.
Variable	Experimento aleatorio (unidad: episodio)	Soporte	Modelo y origen del parámetro	Funciones y cálculo en R
X	Elegir 1 episodio al azar de la cohorte E (n_E = 389) y contar sus campos de diagnóstico (DIAGNOSTICO1 a DIAGNOSTICO35) no vacíos.	{1, 2, …, 35}	Distribución empírica p(x) = n(X=x) / 389	rowSums(!is.na(DIAG) & DIAG != ’’) ; table(X)/nrow(cohorte_E)
Y	Elegir 1 episodio al azar de la cohorte E y observar si termina en fallecimiento (1 si F ocurre, 0 si no).	{0, 1}	Bernoulli(p), con p = P(F \| E) = 16 / 389 = 0,041131	as.integer(F_evt) ; mean(F_evt)
Y_L	Elegir 1 episodio al azar de la cohorte E y observar si registra estancia larga > 7 días (1 si L ocurre, 0 si no).	{0, 1}	Bernoulli(p_L), con p_L = P(L \| E) = 121 / 389 = 0,311054	as.integer(L_evt) ; mean(L_evt)
S	Seleccionar 20 episodios al azar con reposición de la cohorte E y contar cuántos terminan en fallecimiento.	{0, 1, …, 20}	Binomial(n = 20, p = 0,041131)	dbinom(k, 20, p), pbinom(k, 20, p), qbinom(q, 20, p), rbinom(B, 20, p)
S_L	Seleccionar 20 episodios al azar con reposición de la cohorte E y contar cuántos presentan estancia larga (> 7 días).	{0, 1, …, 20}	Binomial(n = 20, p_L = 0,311054)	dbinom(k, 20, p_L), pbinom(k, 20, p_L), qbinom(q, 20, p_L)
H	Auditar n = 10 episodios al azar sin reposición dentro de un establecimiento h (con N_h episodios de E) y contar los casos con el evento de interés (K_h).	{max(0, 10-(N_h-K_h)), …, min(10, K_h)}	Hipergeométrica(N_h, K_h, n = 10)	dhyper(k, m = K_h, n = N_h - K_h, k = 10)
T	Revisar uno por uno episodios marcados por M (con reposición) hasta encontrar el primer episodio perteneciente a la cohorte E.	{1, 2, 3, …}	Geométrica(p = VPP), con VPP = P(E \| M) = 99 / 962 = 0,102911 (y M_orig: 110 / 57.303 = 0,001920)	dgeom(k - 1, prob = VPP), pgeom(k - 1, prob = VPP), qgeom(q, prob = VPP) + 1
T_5	Revisar uno por uno episodios marcados por M (con reposición) hasta acumular r = 5 episodios pertenecientes a la cohorte E.	{5, 6, 7, …}	Binomial Negativa(r = 5, p = VPP)	dnbinom(k - 5, size = 5, prob = VPP), pnbinom(k - 5, size = 5, prob = VPP)
D	Elegir un día calendario del año bisiesto 2024 (366 días) y contar cuántos episodios de la cohorte E egresaron en esa fecha.	{0, 1, 2, …}	Poisson(λ), con λ = 389 / 366 = 1,062842 egresos/día	dpois(k, lambda), ppois(k, lambda), rpois(B, lambda)

Interpretación de la Tabla 2: Todas las variables se definen sobre experimentos cuyos resultados elementales son episodios de egreso hospitalario (o días calendario para el conteo de egresos $D$). Los parámetros de cada distribución se obtienen directamente desde la base depurada mediante proporciones empíricas sobre la cohorte $E$ ($n_E = 389$), sobre los establecimientos $h$ ($N_h$), sobre el subconjunto de fichas marcadas por $M$ ($n(M) = 962$) o sobre los $366$ días del calendario 2024.

1.3 Función de probabilidad empírica y acumulada de $X$

La variable aleatoria $X$ cuenta el número de diagnósticos registrados (campos DIAGNOSTICO1 a DIAGNOSTICO35 no vacíos) en un episodio elegido al azar de la cohorte $E$ ($n_E = 389$). Como todo episodio de la cohorte contiene al menos el código H25, el soporte teórico es $\{1, 2, \dots, 35\}$.

tabla_X_dist <- cohorte_E %>%
  count(X, name = "Episodios") %>%
  right_join(tibble(X = 1:35), by = "X") %>%
  mutate(Episodios = ifelse(is.na(Episodios), 0L, Episodios)) %>%
  arrange(X) %>%
  mutate(
    p_x = Episodios / n_E,
    p_x_pct = sprintf("%.2f%%", p_x * 100),
    F_x = cumsum(p_x),
    F_x_pct = sprintf("%.2f%%", F_x * 100)
  )

# Mostrar los valores del soporte con frecuencia positiva (e indicar que 33 y 34 tienen 0 casos)
tabla_X_mostrar <- tabla_X_dist %>%
  filter(Episodios > 0) %>%
  select(
    `x (N° Diagnósticos)` = X,
    `Episodios observados` = Episodios,
    `p(x) = P(X = x)` = p_x,
    `p(x) (%)` = p_x_pct,
    `F(x) = P(X ≤ x)` = F_x,
    `F(x) (%)` = F_x_pct
  )

kable(
  tabla_X_mostrar,
  digits = 4,
  align = c("c", "r", "r", "r", "r", "r"),
  caption = "Tabla 3: Función de probabilidad empírica p(x) y función de distribución acumulada F(x) de X en la cohorte E (n_E = 389; los valores x = 33 y x = 34 tienen 0 episodios)."
) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)

Tabla 3: Función de probabilidad empírica p(x) y función de distribución acumulada F(x) de X en la cohorte E (n_E = 389; los valores x = 33 y x = 34 tienen 0 episodios).
x (N° Diagnósticos)	Episodios observados	p(x) = P(X = x)	p(x) (%)	F(x) = P(X ≤ x)	F(x) (%)
1	14	0.0360	3.60%	0.0360	3.60%
2	8	0.0206	2.06%	0.0566	5.66%
3	11	0.0283	2.83%	0.0848	8.48%
4	24	0.0617	6.17%	0.1465	14.65%
5	37	0.0951	9.51%	0.2416	24.16%
6	29	0.0746	7.46%	0.3162	31.62%
7	22	0.0566	5.66%	0.3728	37.28%
8	29	0.0746	7.46%	0.4473	44.73%
9	29	0.0746	7.46%	0.5219	52.19%
10	26	0.0668	6.68%	0.5887	58.87%
11	27	0.0694	6.94%	0.6581	65.81%
12	14	0.0360	3.60%	0.6941	69.41%
13	15	0.0386	3.86%	0.7326	73.26%
14	20	0.0514	5.14%	0.7841	78.41%
15	16	0.0411	4.11%	0.8252	82.52%
16	7	0.0180	1.80%	0.8432	84.32%
17	9	0.0231	2.31%	0.8663	86.63%
18	4	0.0103	1.03%	0.8766	87.66%
19	7	0.0180	1.80%	0.8946	89.46%
20	6	0.0154	1.54%	0.9100	91.00%
21	3	0.0077	0.77%	0.9177	91.77%
22	5	0.0129	1.29%	0.9306	93.06%
23	1	0.0026	0.26%	0.9332	93.32%
24	4	0.0103	1.03%	0.9434	94.34%
25	3	0.0077	0.77%	0.9512	95.12%
26	2	0.0051	0.51%	0.9563	95.63%
27	3	0.0077	0.77%	0.9640	96.40%
28	2	0.0051	0.51%	0.9692	96.92%
29	2	0.0051	0.51%	0.9743	97.43%
30	1	0.0026	0.26%	0.9769	97.69%
31	1	0.0026	0.26%	0.9794	97.94%
32	4	0.0103	1.03%	0.9897	98.97%
35	4	0.0103	1.03%	1.0000	100.00%

Interpretación de la Tabla 3: La distribución empírica de $X$ recorre prácticamente todo el rango de campos disponibles en el registro GRD, desde $x = 1$ diagnóstico único ($14$ episodios, $3,60\%$) hasta el tope administrativo de $x = 35$ diagnósticos ($4$ episodios, $1,03\%$), con ausencia de registros únicamente en $x = 33$ y $x = 34$. La masa de probabilidad más alta se alcanza en $x = 5$ diagnósticos (moda de la distribución, con $37$ episodios equivalentes al $9,51\%$ de la cohorte), mientras que la probabilidad acumulada supera el $50\%$ exactamente en $x = 9$ diagnósticos ($F(9) = 52,19\%$, frente a $F(8) = 44,73\%$).

df_plot_X <- bind_rows(
  tabla_X_dist %>% transmute(X = X, Valor = p_x, Panel = "1. Función de Probabilidad Empírica p(x) = P(X = x)"),
  tabla_X_dist %>% transmute(X = X, Valor = F_x, Panel = "2. Función de Distribución Acumulada F(x) = P(X ≤ x)")
)

ggplot(tabla_X_dist, aes(x = X)) +
  geom_col(aes(y = p_x), fill = "#2b5c8f", color = "white", width = 0.8) +
  geom_text(
    data = tabla_X_dist %>% filter(X %in% c(1, 5, 9, 11, 14, 35)),
    aes(y = p_x, label = sprintf("%.1f%%", p_x * 100)),
    vjust = -0.4, size = 3, fontface = "bold", color = "#1f3a52"
  ) +
  scale_x_continuous(breaks = seq(1, 35, by = 2)) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1), limits = c(0, 0.115)) +
  theme_minimal(base_size = 11) +
  labs(
    title = "Figura 1A: Función de probabilidad empírica p(x) del número de diagnósticos (X)",
    subtitle = "Cohorte Catarata Senil (H25), n_E = 389 episodios hospitalarios GRD 2024",
    x = "Número de diagnósticos registrados por episodio (x)",
    y = "Probabilidad puntual p(x) = P(X = x)"
  )

ggplot(tabla_X_dist, aes(x = X, y = F_x)) +
  geom_step(color = "#c0392b", linewidth = 1.1) +
  geom_point(data = tabla_X_dist %>% filter(Episodios > 0), color = "#c0392b", size = 1.8) +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "gray40") +
  annotate("text", x = 26, y = 0.54, label = "Mediana: x = 9 (F(9) = 52,19%)", color = "gray20", size = 3.5, fontface = "bold") +
  scale_x_continuous(breaks = seq(1, 35, by = 2)) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1), limits = c(0, 1.05)) +
  theme_minimal(base_size = 11) +
  labs(
    title = "Figura 1B: Función de distribución acumulada empírica F(x) = P(X ≤ x)",
    subtitle = "Crecimiento escalonado sobre el soporte {1, ..., 35} en la cohorte E (n_E = 389)",
    x = "Número de diagnósticos registrados por episodio (x)",
    y = "Probabilidad acumulada F(x) = P(X ≤ x)"
  )

Interpretación de las Figuras 1A y 1B: El gráfico de barras de $p(x)$ (Figura 1A) evidencia una distribución unimodal pero con marcada asimetría positiva (cola derecha extendida hacia valores altos de comorbilidad), concentrando el grueso de los episodios entre $4$ y $15$ diagnósticos. La función escalonada $F(x)$ (Figura 1B) muestra visualmente el cruce del umbral del $50\%$ en $x = 9$ diagnósticos y un ascenso más suave a partir de $x = 16$, donde se ubican los episodios pluripatológicos de alta complejidad.

1.4 Esperanza, varianza, mediana, moda y probabilidad de cola de $X$

Calculamos la esperanza $E[X] = \sum_{x=1}^{35} x \cdot p(x)$ y la varianza $\text{Var}[X] = \sum_{x=1}^{35} (x - E[X])^2 \cdot p(x)$ directamente desde la función de probabilidad empírica y las contrastamos con las funciones mean() y var() aplicadas sobre los $389$ episodios. Asimismo, evaluamos la mediana, la moda y la probabilidad de registrar una carga elevada de diagnósticos concomitantes: $P(X \ge 10)$ (10 o más diagnósticos registrados en el episodio) y $P(X \ge 5)$ (al menos 5 diagnósticos).

EX_formula <- sum(tabla_X_dist$X * tabla_X_dist$p_x)
VarX_formula <- sum((tabla_X_dist$X - EX_formula)^2 * tabla_X_dist$p_x)

EX_r <- mean(cohorte_E$X)
VarX_r <- var(cohorte_E$X)
factor_bessel <- n_E / (n_E - 1)

mediana_X <- median(cohorte_E$X)
moda_X <- tabla_X_dist$X[which.max(tabla_X_dist$Episodios)]
prob_X_ge_5 <- sum(tabla_X_dist$p_x[tabla_X_dist$X >= 5])
prob_X_ge_10 <- sum(tabla_X_dist$p_x[tabla_X_dist$X >= 10])

tabla_resumen_X <- tibble(
  Metrica = c(
    "Esperanza E[X] = Σ x · p(x)",
    "Promedio en R: mean(X)",
    "Diferencia absoluta |E[X] - mean(X)|",
    "Varianza de la distribución Var[X] = Σ (x - E[X])² · p(x) [divide por n = 389]",
    "Varianza muestral en R: var(X) [divide por n - 1 = 388]",
    "Diferencia |Var[X] - var(X)| verificada por factor n/(n-1) = 389/388",
    "Mediana de X",
    "Moda de X (valor más frecuente)",
    "P(X ≥ 5): Episodio con 5 o más diagnósticos registrados",
    "P(X ≥ 10): Episodio con 10 o más diagnósticos (alta pluripatología)"
  ),
  Valor = c(
    sprintf("%.6f diagnósticos", EX_formula),
    sprintf("%.6f diagnósticos", EX_r),
    sprintf("%.6f", abs(EX_formula - EX_r)),
    sprintf("%.6f diagnósticos²", VarX_formula),
    sprintf("%.6f diagnósticos²", VarX_r),
    sprintf("%.6f (Var[X] × 389/388 = %.6f)", abs(VarX_formula - VarX_r), VarX_formula * factor_bessel),
    sprintf("%d diagnósticos", as.integer(mediana_X)),
    sprintf("%d diagnósticos (%d episodios; %.2f%%)", moda_X, max(tabla_X_dist$Episodios), max(tabla_X_dist$p_x) * 100),
    sprintf("%.4f (%.2f%%; %d de 389 episodios)", prob_X_ge_5, prob_X_ge_5 * 100, sum(cohorte_E$X >= 5)),
    sprintf("%.4f (%.2f%%; %d de 389 episodios)", prob_X_ge_10, prob_X_ge_10 * 100, sum(cohorte_E$X >= 10))
  )
)

kable(
  tabla_resumen_X,
  col.names = c("Estadístico / Propiedad de la variable aleatoria X", "Valor calculado sobre la cohorte E (n_E = 389)"),
  caption = "Tabla 4: Verificación algebraica de E[X] y Var[X], medidas de posición y probabilidades de cola de X."
) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)

Tabla 4: Verificación algebraica de E[X] y Var[X], medidas de posición y probabilidades de cola de X.
Estadístico / Propiedad de la variable aleatoria X	Valor calculado sobre la cohorte E (n_E = 389)
Esperanza E[X] = Σ x · p(x)	10.637532 diagnósticos
Promedio en R: mean(X)	10.637532 diagnósticos
Diferencia absoluta \|E[X] - mean(X)\|	0.000000
Varianza de la distribución Var[X] = Σ (x - E[X])² · p(x) [divide por n = 389]	47.336483 diagnósticos²
Varianza muestral en R: var(X) [divide por n - 1 = 388]	47.458485 diagnósticos²
Diferencia \|Var[X] - var(X)\| verificada por factor n/(n-1) = 389/388	0.122001 (Var[X] × 389/388 = 47.458485)
Mediana de X	9 diagnósticos
Moda de X (valor más frecuente)	5 diagnósticos (37 episodios; 9.51%)
P(X ≥ 5): Episodio con 5 o más diagnósticos registrados	0.8535 (85.35%; 332 de 389 episodios)
P(X ≥ 10): Episodio con 10 o más diagnósticos (alta pluripatología)	0.4781 (47.81%; 186 de 389 episodios)

Interpretación de la Tabla 4 y explicación de la diferencia en var():

Esperanza y Varianza: La esperanza calculada desde la distribución empírica ($\sum x \cdot p(x) = 10,637532$) coincide exactamente con mean(cohorte_E$X). En cambio, la varianza calculada desde la función de probabilidad es $\text{Var}[X] = 47,336483$, mientras que var(cohorte_E$X) entrega $47,458485$ (una diferencia de $0,122002$). Esto ocurre porque la función de probabilidad empírica trata la cohorte de $n_E = 389$ episodios como la población de referencia completa y divide la suma de cuadrados de las desviaciones por $n_E = 389$: \[\text{Var}[X] = \frac{1}{n_E}\sum_{i=1}^{n_E}(x_i - \bar{x})^2 = 47,336483\] mientras que var() en R implementa el estimador muestral insesgado aplicando la corrección de Bessel, que divide por $n_E - 1 = 388$ grados de libertad: \[s_X^2 = \frac{1}{n_E - 1}\sum_{i=1}^{n_E}(x_i - \bar{x})^2 = \text{Var}[X] \cdot \frac{389}{388} = 47,336483 \times 1,0025773 = 47,458485\]
Posición y probabilidades con sentido clínico: Que la media ($10,64$) supere a la mediana ($9$) y a la moda ($5$) confirma el sesgo positivo de $X$. En términos operativos, elegir un episodio de la cohorte al azar tiene un $85,35\%$ de probabilidad de registrar al menos $5$ diagnósticos ($P(X \ge 5) = 332/389$) y un $47,81\%$ de probabilidad de presentar $10$ o más diagnósticos ($P(X \ge 10) = 186/389$), lo que refleja que casi la mitad de los egresos hospitalarios con código H25 corresponden a episodios de alta complejidad diagnóstica.

1.5 Comparación de la función de probabilidad de $X$ según rol diagnóstico ($E_p$ vs $E_s$)

En la Actividad 1 y 2 se demostró que los episodios donde la Catarata Senil es el diagnóstico principal ($E_p$, $n = 108$) difieren sustancialmente de aquellos donde se registra como diagnóstico secundario ($E_s$, $n = 281$). A continuación se construye y compara la función de probabilidad de $X$ por separado para ambas subcohortes.

resumen_rol_X <- cohorte_E %>%
  mutate(Rol = ifelse(E_p, "Ep: Diagnóstico Principal (n = 108)", "Es: Diagnóstico Secundario (n = 281)")) %>%
  group_by(Rol) %>%
  summarise(
    Episodios = n(),
    Esperanza = mean(X),
    Varianza_Pob = mean((X - mean(X))^2),
    Varianza_Muestral = var(X),
    Mediana = median(X),
    Moda = as.integer(names(which.max(table(X)))),
    Min_Max = paste0(min(X), " a ", max(X)),
    `P(X ≥ 5)` = sprintf("%.2f%%", mean(X >= 5) * 100),
    `P(X ≥ 10)` = sprintf("%.2f%%", mean(X >= 10) * 100),
    .groups = "drop"
  )

kable(
  resumen_rol_X,
  digits = 2,
  col.names = c("Subcohorte según rol diagnóstico", "n", "E[X]", "Var[X] (div n)", "var(X) (div n-1)", "Mediana", "Moda", "Rango observado", "P(X ≥ 5)", "P(X ≥ 10)"),
  caption = "Tabla 5: Comparación de parámetros y forma de la distribución de X entre episodios con H25 principal (Ep) y secundario (Es)."
) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = TRUE)

Tabla 5: Comparación de parámetros y forma de la distribución de X entre episodios con H25 principal (Ep) y secundario (Es).
Subcohorte según rol diagnóstico	n	E[X]	Var[X] (div n)	var(X) (div n-1)	Mediana	Moda	Rango observado	P(X ≥ 5)	P(X ≥ 10)
Ep: Diagnóstico Principal (n = 108)	108	5.78	11.06	11.17	5	5	1 a 16	62.96%	12.04%
Es: Diagnóstico Secundario (n = 281)	281	12.51	48.71	48.89	11	11	2 a 35	93.95%	61.57%

dist_Ep_Es <- cohorte_E %>%
  mutate(Rol = ifelse(E_p, "Ep: H25 Principal (n = 108)", "Es: H25 Secundario (n = 281)")) %>%
  count(Rol, X, name = "n") %>%
  group_by(Rol) %>%
  mutate(p_x = n / sum(n)) %>%
  ungroup()

ggplot(dist_Ep_Es, aes(x = X, y = p_x, fill = Rol)) +
  geom_col(position = position_dodge(preserve = "single"), width = 0.75, color = "white") +
  scale_fill_manual(values = c("Ep: H25 Principal (n = 108)" = "#2b5c8f", "Es: H25 Secundario (n = 281)" = "#e67e22")) +
  scale_x_continuous(breaks = seq(1, 35, by = 2)) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  theme_minimal(base_size = 11) +
  theme(legend.position = "top") +
  labs(
    title = "Figura 2: Comparación de la función de probabilidad de X entre Ep (Principal) y Es (Secundario)",
    subtitle = "Ep se concentra en pocos diagnósticos (1 a 16), mientras Es se desplaza hacia alta pluripatología (2 a 35)",
    x = "Número de diagnósticos registrados en el episodio (X)",
    y = "Probabilidad condicional P(X = x | Rol)",
    fill = "Subcohorte:"
  )

Interpretación de la Tabla 5 y Figura 2: Separar la variable $X$ por rol diagnóstico revela que la distribución global de la cohorte es una mezcla de dos regímenes hospitalarios distintos:

En $E_p$ (H25 como diagnóstico principal, $n = 108$), los episodios corresponden a hospitalizaciones oftalmológicas acotadas: el soporte observado llega solo hasta $16$ diagnósticos, la esperanza es de $E[X \mid E_p] = 5,78$ diagnósticos ($\text{Var}[X \mid E_p] = 11,06$), la mediana y la moda coinciden en $5$, existe un $12,96\%$ de episodios con diagnóstico único ($x = 1$, $14$ casos) y solo el $12,04\%$ registra $10$ o más diagnósticos.
En $E_s$ (H25 como diagnóstico secundario, $n = 281$), por definición todo episodio tiene al menos $x = 2$ diagnósticos (la causa principal de internación más la catarata senil): la esperanza más que se duplica hasta $E[X \mid E_s] = 12,51$ diagnósticos ($\text{Var}[X \mid E_s] = 48,71$), la mediana y la moda suben a $11$ diagnósticos, la cola derecha alcanza los $35$ campos disponibles y el $61,57\%$ de los episodios registra $10$ o más diagnósticos.

2 Resultado 02: Distribuciones Bernoulli y Binomial

2.1 Parámetros de Bernoulli ($Y$ para fallecimiento y $Y_L$ para estancia larga)

Definimos sobre la cohorte $E$ ($n_E = 389$ episodios) la variable aleatoria indicadora $Y \sim \text{Bernoulli}(p)$, que toma el valor $1$ si un episodio elegido al azar de $E$ termina en fallecimiento ($F$) y $0$ en caso contrario. De manera complementaria —y siguiendo la indicación metodológica para cohortes con pocos fallecimientos—, definimos $Y_L \sim \text{Bernoulli}(p_L)$, que vale $1$ si el episodio presenta una estancia larga ($L$: estadía superior al percentil 75 de $\Omega$, correspondiente a $> 7$ días) y $0$ si no.

n_F_cohorte <- sum(cohorte_E$F_evt)
p_F <- n_F_cohorte / n_E
EY <- p_F
VarY <- p_F * (1 - p_F)

n_L_cohorte <- sum(cohorte_E$L_evt)
p_L <- n_L_cohorte / n_E
EY_L <- p_L
VarY_L <- p_L * (1 - p_L)

tabla_bernoulli <- tibble(
  Variable = c(
    "Y ~ Bernoulli(p) [Fallecimiento al egreso F]",
    "Y_L ~ Bernoulli(p_L) [Estancia larga L > 7 días]"
  ),
  Origen_Parametro = c(
    sprintf("p = P(F | E) = %d / %d", n_F_cohorte, n_E),
    sprintf("p_L = P(L | E) = %d / %d", n_L_cohorte, n_E)
  ),
  P_Y_1 = c(sprintf("%.6f (%.2f%%)", p_F, p_F * 100), sprintf("%.6f (%.2f%%)", p_L, p_L * 100)),
  P_Y_0 = c(sprintf("%.6f (%.2f%%)", 1 - p_F, (1 - p_F) * 100), sprintf("%.6f (%.2f%%)", 1 - p_L, (1 - p_L) * 100)),
  Esperanza = c(sprintf("%.6f", EY), sprintf("%.6f", EY_L)),
  Varianza = c(sprintf("%.6f", VarY), sprintf("%.6f", VarY_L))
)

kable(
  tabla_bernoulli,
  col.names = c("Variable Bernoulli", "Cálculo del parámetro sobre E", "P(Y = 1) = p", "P(Y = 0) = 1 - p", "E[Y] = p", "Var[Y] = p(1 - p)"),
  align = c("l", "c", "r", "r", "r", "r"),
  caption = "Tabla 6: Parámetros, probabilidades puntuales, esperanza y varianza de los modelos Bernoulli para fallecimiento (Y) y estancia larga (Y_L)."
) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)

Tabla 6: Parámetros, probabilidades puntuales, esperanza y varianza de los modelos Bernoulli para fallecimiento (Y) y estancia larga (Y_L).
Variable Bernoulli	Cálculo del parámetro sobre E	P(Y = 1) = p	P(Y = 0) = 1 - p	E[Y] = p	Var[Y] = p(1 - p)
Y ~ Bernoulli(p) [Fallecimiento al egreso F]	p = P(F \| E) = 16 / 389	0.041131 (4.11%)	0.958869 (95.89%)	0.041131	0.039439
Y_L ~ Bernoulli(p_L) [Estancia larga L > 7 días]	p_L = P(L \| E) = 121 / 389	0.311054 (31.11%)	0.688946 (68.89%)	0.311054	0.214299

Interpretación de la Tabla 6: El parámetro de fallecimiento en la cohorte es $p = P(F \mid E) = 16 / 389 = 0,041131$ ($4,11\%$), por lo que para un episodio extraído al azar de $E$ la esperanza es $E[Y] = 0,041131$ y su varianza es $\text{Var}[Y] = 0,041131 \times 0,958869 = 0,039439$. En contraste, la estancia larga ($> 7$ días) es un evento casi ocho veces más frecuente en la cohorte, con $p_L = P(L \mid E) = 121 / 389 = 0,311054$ ($31,11\%$), lo que eleva su esperanza a $E[Y_L] = 0,311054$ y su varianza unitaria a $\text{Var}[Y_L] = 0,214300$ (al acercarse $p_L$ a $0,5$, la incertidumbre por ensayo aumenta).

2.2 Modelo Binomial en una auditoría de $n = 20$ episodios ($S$ y $S_L$)

Consideremos el experimento de auditar $n = 20$ episodios de la cohorte $E$ elegidos al azar con reposición:

$S \sim \text{Binomial}(20, p = 0,041131)$ cuenta el número de fallecimientos entre los $20$ episodios auditados.
$S_L \sim \text{Binomial}(20, p_L = 0,311054)$ cuenta el número de episodios con estancia larga ($> 7$ días) entre los $20$ episodios auditados.

n_bin <- 20
k_bin <- 0:n_bin

df_binom <- tibble(
  k = k_bin,
  p_S = dbinom(k_bin, size = n_bin, prob = p_F),
  F_S = pbinom(k_bin, size = n_bin, prob = p_F),
  p_SL = dbinom(k_bin, size = n_bin, prob = p_L),
  F_SL = pbinom(k_bin, size = n_bin, prob = p_L)
)

# Tabla para k = 0 hasta 15 (para k >= 16 ambas probabilidades puntuales son < 0.0001)
df_binom_mostrar <- df_binom %>%
  filter(k <= 15) %>%
  transmute(
    `k (Conteo en 20 episodios)` = k,
    `P(S = k) [Fallecimientos]` = sprintf("%.6f (%.2f%%)", p_S, p_S * 100),
    `P(S ≤ k)` = sprintf("%.6f (%.2f%%)", F_S, F_S * 100),
    `P(S_L = k) [Estancia Larga]` = sprintf("%.6f (%.2f%%)", p_SL, p_SL * 100),
    `P(S_L ≤ k)` = sprintf("%.6f (%.2f%%)", F_SL, F_SL * 100)
  )

kable(
  df_binom_mostrar,
  align = c("c", "r", "r", "r", "r"),
  caption = "Tabla 7: Función de probabilidad puntual (dbinom) y acumulada (pbinom) de S ~ Binomial(20; 0,0411) y S_L ~ Binomial(20; 0,3111) para k = 0 a 15 (para k ≥ 16, P < 0,004%)."
) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)

Tabla 7: Función de probabilidad puntual (dbinom) y acumulada (pbinom) de S ~ Binomial(20; 0,0411) y S_L ~ Binomial(20; 0,3111) para k = 0 a 15 (para k ≥ 16, P < 0,004%).
k (Conteo en 20 episodios)	P(S = k) [Fallecimientos]	P(S ≤ k)	P(S_L = k) [Estancia Larga]	P(S_L ≤ k)
0	0.431703 (43.17%)	0.431703 (43.17%)	0.000580 (0.06%)	0.000580 (0.06%)
1	0.370361 (37.04%)	0.802064 (80.21%)	0.005241 (0.52%)	0.005821 (0.58%)
2	0.150925 (15.09%)	0.952989 (95.30%)	0.022478 (2.25%)	0.028299 (2.83%)
3	0.038844 (3.88%)	0.991833 (99.18%)	0.060892 (6.09%)	0.089191 (8.92%)
4	0.007081 (0.71%)	0.998914 (99.89%)	0.116842 (11.68%)	0.206033 (20.60%)
5	0.000972 (0.10%)	0.999886 (99.99%)	0.168811 (16.88%)	0.374844 (37.48%)
6	0.000104 (0.01%)	0.999990 (100.00%)	0.190542 (19.05%)	0.565386 (56.54%)
7	0.000009 (0.00%)	0.999999 (100.00%)	0.172057 (17.21%)	0.737443 (73.74%)
8	0.000001 (0.00%)	1.000000 (100.00%)	0.126234 (12.62%)	0.863677 (86.37%)
9	0.000000 (0.00%)	1.000000 (100.00%)	0.075991 (7.60%)	0.939668 (93.97%)
10	0.000000 (0.00%)	1.000000 (100.00%)	0.037741 (3.77%)	0.977409 (97.74%)
11	0.000000 (0.00%)	1.000000 (100.00%)	0.015491 (1.55%)	0.992899 (99.29%)
12	0.000000 (0.00%)	1.000000 (100.00%)	0.005245 (0.52%)	0.998145 (99.81%)
13	0.000000 (0.00%)	1.000000 (100.00%)	0.001457 (0.15%)	0.999602 (99.96%)
14	0.000000 (0.00%)	1.000000 (100.00%)	0.000329 (0.03%)	0.999931 (99.99%)
15	0.000000 (0.00%)	1.000000 (100.00%)	0.000059 (0.01%)	0.999991 (100.00%)

# Resumen de indicadores exigidos para S y S_L
P_S_0 <- dbinom(0, size = n_bin, prob = p_F)
P_S_ge_2 <- 1 - pbinom(1, size = n_bin, prob = p_F)
q90_S <- qbinom(0.90, size = n_bin, prob = p_F)
ES <- n_bin * p_F
VarS <- n_bin * p_F * (1 - p_F)

P_SL_0 <- dbinom(0, size = n_bin, prob = p_L)
P_SL_ge_2 <- 1 - pbinom(1, size = n_bin, prob = p_L)
q90_SL <- qbinom(0.90, size = n_bin, prob = p_L)
ESL <- n_bin * p_L
VarSL <- n_bin * p_L * (1 - p_L)

tabla_indicadores_binom <- tibble(
  Indicador = c(
    "Parámetro p estimado sobre la cohorte E (n_E = 389)",
    "P(S = 0) = (1 - p)²⁰ [Ningún evento en los 20 episodios]",
    "P(S ≥ 2) = 1 - P(S ≤ 1) [Dos o más eventos en los 20 episodios]",
    "Percentil 90: qbinom(0.90, size = 20, prob = p)",
    "Esperanza teórica E[S] = 20 · p",
    "Varianza teórica Var[S] = 20 · p · (1 - p)"
  ),
  Fallecimientos_S = c(
    sprintf("p = 16 / 389 = %.6f (4,11%%)", p_F),
    sprintf("%.6f (%.2f%%)", P_S_0, P_S_0 * 100),
    sprintf("%.6f (%.2f%%)", P_S_ge_2, P_S_ge_2 * 100),
    sprintf("%d fallecimientos (P(S ≤ 2) = %.2f%%)", q90_S, pbinom(q90_S, 20, p_F) * 100),
    sprintf("%.4f fallecimientos", ES),
    sprintf("%.4f fallecimientos²", VarS)
  ),
  Estancia_Larga_SL = c(
    sprintf("p_L = 121 / 389 = %.6f (31,11%%)", p_L),
    sprintf("%.6f (%.4f%%)", P_SL_0, P_SL_0 * 100),
    sprintf("%.6f (%.2f%%)", P_SL_ge_2, P_SL_ge_2 * 100),
    sprintf("%d episodios (P(S_L ≤ 9) = %.2f%%)", q90_SL, pbinom(q90_SL, 20, p_L) * 100),
    sprintf("%.4f episodios con L", ESL),
    sprintf("%.4f episodios²", VarSL)
  )
)

kable(
  tabla_indicadores_binom,
  col.names = c("Métrica exigida (n = 20 ensayos)", "S ~ Binomial(20; 0,0411) [Fallecimiento F]", "S_L ~ Binomial(20; 0,3111) [Estancia larga L]"),
  align = c("l", "r", "r"),
  caption = "Tabla 8: Comparación de probabilidades puntuales, acumuladas, percentil 90, esperanza y varianza entre S (fallecimiento) y S_L (estancia larga)."
) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)

Tabla 8: Comparación de probabilidades puntuales, acumuladas, percentil 90, esperanza y varianza entre S (fallecimiento) y S_L (estancia larga).
Métrica exigida (n = 20 ensayos)	S ~ Binomial(20; 0,0411) [Fallecimiento F]	S_L ~ Binomial(20; 0,3111) [Estancia larga L]
Parámetro p estimado sobre la cohorte E (n_E = 389)	p = 16 / 389 = 0.041131 (4,11%)	p_L = 121 / 389 = 0.311054 (31,11%)
P(S = 0) = (1 - p)²⁰ [Ningún evento en los 20 episodios]	0.431703 (43.17%)	0.000580 (0.0580%)
P(S ≥ 2) = 1 - P(S ≤ 1) [Dos o más eventos en los 20 episodios]	0.197936 (19.79%)	0.994179 (99.42%)
Percentil 90: qbinom(0.90, size = 20, prob = p)	2 fallecimientos (P(S ≤ 2) = 95.30%)	9 episodios (P(S_L ≤ 9) = 93.97%)
Esperanza teórica E[S] = 20 · p	0.8226 fallecimientos	6.2211 episodios con L
Varianza teórica Var[S] = 20 · p · (1 - p)	0.7888 fallecimientos²	4.2860 episodios²

Interpretación de las Tablas 7 y 8:

Para los fallecimientos ($S \sim \text{Binomial}(20; 0,041131)$): En una auditoría aleatoria de $20$ episodios de la cohorte, el número esperado de fallecimientos es de $E[S] = 0,8226$ con varianza $\text{Var}[S] = 0,7888$. El resultado individual más probable es no encontrar ningún fallecimiento ($P(S = 0) = 0,958869^{20} = 43,17\%$), seguido de encontrar exactamente uno ($P(S = 1) = 37,04\%$). La probabilidad de observar $2$ o más fallecimientos es de $P(S \ge 2) = 19,79\%$, y el percentil 90 es $Q_{90} = 2$ fallecimientos (ya que el $95,21\%$ de las auditorías de 20 fichas tendrá a lo sumo 2 fallecidos).
Para las estancias largas ($S_L \sim \text{Binomial}(20; 0,311054)$): Al evaluar estancias superiores a $7$ días, el conteo esperado sube a $E[S_L] = 6,2211$ episodios con varianza $\text{Var}[S_L] = 4,2860$. Aquí es prácticamente imposible no observar casos ($P(S_L = 0) = 0,0580\%$), la probabilidad de encontrar al menos $2$ episodios con estancia larga es de $P(S_L \ge 2) = 99,42\%$, y el percentil 90 se ubica en $Q_{90} = 9$ episodios ($P(S_L \le 9) = 93,75\%$).

df_binom_long <- df_binom %>%
  filter(k <= 15) %>%
  select(k, `S: Fallecimientos (p = 4,11%)` = p_S, `S_L: Estancia Larga > 7 días (p_L = 31,11%)` = p_SL) %>%
  pivot_longer(cols = -k, names_to = "Modelo", values_to = "Probabilidad")

ggplot(df_binom_long, aes(x = k, y = Probabilidad, fill = Modelo)) +
  geom_col(position = "dodge", width = 0.75, color = "white") +
  scale_fill_manual(values = c("S: Fallecimientos (p = 4,11%)" = "#c0392b", "S_L: Estancia Larga > 7 días (p_L = 31,11%)" = "#2980b9")) +
  scale_x_continuous(breaks = 0:15) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  theme_minimal(base_size = 11) +
  theme(legend.position = "top") +
  labs(
    title = "Figura 3: Comparación de la forma de S ~ Binomial(20; 0,0411) y S_L ~ Binomial(20; 0,3111)",
    subtitle = "Efecto del parámetro p sobre la asimetría, localización y dispersión en n = 20 ensayos",
    x = "Número de episodios con el evento entre 20 auditados (k)",
    y = "Probabilidad puntual P(Variable = k)",
    fill = "Distribución Binomial:"
  )

Interpretación de la Figura 3 (Comparación de la forma de ambas distribuciones): La forma de una distribución $\text{Binomial}(n, p)$ con $n = 20$ fijo depende críticamente de la cercanía de $p$ a los bordes $0$ o $0,5$. Como $p_F = 0,0411$ es cercano a cero, la distribución de $S$ (barras rojas) presenta una asimetría positiva extrema ($\gamma_1 = \frac{1-2p}{\sqrt{np(1-p)}} = +1,03$), truncada en el borde izquierdo $k = 0$ (su moda) y concentrando el $95,21\%$ de su masa en $\{0, 1, 2\}$. En cambio, como $p_L = 0,3111$ es mucho más cercano a $0,5$, la distribución de $S_L$ (barras azules) se desplaza hacia el centro del soporte, adopta una forma acampanada casi simétrica ($\gamma_1 = +0,18$) centrada en su moda $k = 6$ ($P(S_L = 6) = 19,16\%$) y presenta una dispersión horizontal mucho mayor ($\text{Var}[S_L] = 4,29$ vs $\text{Var}[S] = 0,79$).

2.3 Comparación por establecimientos y probabilidades de cola binomial

En la cohorte $E$ ($389$ episodios distribuidos en $61$ recintos hospitalarios), existen $5$ establecimientos con al menos $20$ episodios (y $9$ establecimientos con al menos $10$ episodios). Atendiendo la instrucción expresa de la guía (“Si p es muy pequeña o hay pocos establecimientos con 20 episodios, trabaje además con la estancia larga L, que suele ser más frecuente, e indíquelo en el texto”), la Tabla 9 presenta los $5$ establecimientos con $n_h \ge 20$ episodios y extiende la vista a los $4$ establecimientos adicionales con $10 \le n_h < 20$, reportando tanto para fallecimientos ($F$) como para estancias largas ($L$) el conteo observado $k_{\text{obs}}$, el valor esperado bajo homogeneidad $n_h \cdot p$ y la probabilidad binomial de cola superior $P(\text{Binomial}(n_h, p) \ge k_{\text{obs}}) = 1 - \text{pbinom}(k_{\text{obs}} - 1, n_h, p)$.

tabla_estab_binom <- cohorte_E %>%
  group_by(COD_HOSPITAL, NOMBRE_HOSPITAL) %>%
  summarise(
    n_h = n(),
    pct_Es = mean(E_s) * 100,
    obs_F = sum(F_evt),
    esp_F = n_h * p_F,
    cola_F = ifelse(obs_F == 0, 1.0, pbinom(obs_F - 1, size = n_h, prob = p_F, lower.tail = FALSE)),
    obs_L = sum(L_evt),
    esp_L = n_h * p_L,
    cola_L = ifelse(obs_L == 0, 1.0, pbinom(obs_L - 1, size = n_h, prob = p_L, lower.tail = FALSE)),
    .groups = "drop"
  ) %>%
  filter(n_h >= 10) %>%
  arrange(desc(n_h)) %>%
  mutate(
    Grupo_Volumen = ifelse(n_h >= 20, "≥ 20 episodios (Criterio principal)", "10 a 19 episodios (Extensión)"),
    pct_Es_fmt = sprintf("%.1f%%", pct_Es),
    esp_F_fmt = sprintf("%.2f", esp_F),
    cola_F_fmt = sprintf("%.4f (%.2f%%)", cola_F, cola_F * 100),
    esp_L_fmt = sprintf("%.2f", esp_L),
    cola_L_fmt = sprintf("%.4f (%.2f%%)", cola_L, cola_L * 100)
  )

kable(
  tabla_estab_binom %>%
    select(
      Código = COD_HOSPITAL,
      `Establecimiento (Tabla Maestra GRD)` = NOMBRE_HOSPITAL,
      `Tramo` = Grupo_Volumen,
      `n_h` = n_h,
      `% H25 Sec. (Es)` = pct_Es_fmt,
      `Obs. F` = obs_F,
      `Esp. F (n_h·p)` = esp_F_fmt,
      `P(Binom ≥ Obs. F)` = cola_F_fmt,
      `Obs. L` = obs_L,
      `Esp. L (n_h·p_L)` = esp_L_fmt,
      `P(Binom ≥ Obs. L)` = cola_L_fmt
    ),
  align = c("c", "l", "l", "r", "r", "r", "r", "r", "r", "r", "r"),
  caption = "Tabla 9: Conteos observados, esperados bajo homogeneidad y probabilidades binomiales de cola superior P(X ≥ obs) para fallecimientos (F) y estancias largas (L) en establecimientos con n_h ≥ 20 (y extensión a n_h ≥ 10)."
) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = TRUE)

Tabla 9: Conteos observados, esperados bajo homogeneidad y probabilidades binomiales de cola superior P(X ≥ obs) para fallecimientos (F) y estancias largas (L) en establecimientos con n_h ≥ 20 (y extensión a n_h ≥ 10).
Código	Establecimiento (Tabla Maestra GRD)	Tramo	n_h	% H25 Sec. (Es)	Obs. F	Esp. F (n_h·p)	P(Binom ≥ Obs. F)	Obs. L	Esp. L (n_h·p_L)	P(Binom ≥ Obs. L)
118100	Hospital Clínico Regional Dr. Guillermo Grant Benavente (Concepción)	≥ 20 episodios (Criterio principal)	83	69.9%	0	3.41	1.0000 (100.00%)	19	25.82	0.9618 (96.18%)
123100	Hospital Base San José de Osorno	≥ 20 episodios (Criterio principal)	47	6.4%	0	1.93	1.0000 (100.00%)	2	14.62	1.0000 (100.00%)
109100	Complejo Hospitalario San José (Santiago, Independencia)	≥ 20 episodios (Criterio principal)	40	100.0%	3	1.65	0.2264 (22.64%)	16	12.44	0.1484 (14.84%)
116100	Hospital San Juan de Dios (Curicó)	≥ 20 episodios (Criterio principal)	29	100.0%	3	1.19	0.1152 (11.52%)	14	9.02	0.0395 (3.95%)
116105	Hospital Dr. César Garavagno Burotto (Talca)	≥ 20 episodios (Criterio principal)	22	100.0%	1	0.90	0.6031 (60.31%)	7	6.84	0.5510 (55.10%)
114105	Hospital Clínico Metropolitano La Florida Dra. Eloisa Díaz Inzunza	10 a 19 episodios (Extensión)	17	100.0%	3	0.70	0.0307 (3.07%)	10	5.29	0.0166 (1.66%)
120101	Complejo Asistencial Dr. Víctor Ríos Ruiz (Los Ángeles)	10 a 19 episodios (Extensión)	17	11.8%	0	0.70	1.0000 (100.00%)	2	5.29	0.9846 (98.46%)
110100	Hospital San Juan de Dios (Santiago, Santiago)	10 a 19 episodios (Extensión)	16	93.8%	1	0.66	0.4893 (48.93%)	6	4.98	0.3771 (37.71%)
122100	Hospital Clínico Regional (Valdivia)	10 a 19 episodios (Extensión)	13	92.3%	0	0.53	1.0000 (100.00%)	8	4.04	0.0228 (2.28%)

Interpretación de la Tabla 9 y discusión de los supuestos del modelo Binomial:

Lectura descriptiva de las probabilidades de cola (no son pruebas de hipótesis ni indicadores de calidad):
- Entre los $5$ establecimientos con $n_h \ge 20$ episodios, el Hospital San Juan de Dios de Curicó ($n_h = 29$) observa $3$ fallecimientos frente a $1,19$ esperados ($P(\ge 3) = 11,52\%$) y $14$ estancias largas frente a $9,02$ esperadas ($P(\ge 14) = 3,95\%$), mientras que el Complejo Hospitalario San José ($n_h = 40$) observa $3$ fallecidos frente a $1,65$ esperados ($P(\ge 3) = 22,64\%$) y $16$ estancias largas frente a $12,44$ esperadas ($P(\ge 16) = 14,84\%$). En el tramo extendido ($n_h \ge 10$), el Hospital Clínico Metropolitano La Florida ($n_h = 17$) registra $3$ fallecidos frente a $0,70$ esperados ($P(\ge 3) = 3,07\%$) y $10$ estancias largas frente a $5,29$ esperadas ($P(\ge 10) = 1,66\%$).
- En el extremo opuesto, los dos hospitales con mayor número de episodios —Hospital Dr. Guillermo Grant Benavente de Concepción ($n_h = 83$, esperado $3,41$ fallecidos y $25,82$ estancias largas) y Hospital Base San José de Osorno ($n_h = 47$, esperado $1,93$ fallecidos y $14,62$ estancias largas)— registran $0$ fallecimientos ($P(\ge 0) = 100\%$) y solo $19$ y $2$ estancias largas respectivamente (en Osorno, $P(\ge 2) = 99,9999\%$, es decir, $P(\le 2) = 0,0017\%$).
- Tal como advierte la pauta, estas probabilidades de cola no constituyen pruebas de hipótesis ni evidencia de buena o mala atención clínica: únicamente cuantifican qué tan inusual sería el conteo si todos los hospitales compartieran exactamente el mismo parámetro global $p$.
Evaluación de los dos supuestos de la distribución Binomial en los datos GRD:
- Supuesto 1: Probabilidad de éxito $p$ idéntica en todos los episodios (homogeneidad). NO SE CUMPLE. La columna % H25 Sec. (Es) de la Tabla 9 explica por completo las diferencias entre establecimientos: en toda la cohorte $E$, los $16$ fallecimientos ocurren exclusivamente cuando la catarata senil es un diagnóstico secundario ($P(F \mid E_p) = 0/108 = 0,00\%$ frente a $P(F \mid E_s) = 16/281 = 5,69\%$), y la estancia larga sube de $P(L \mid E_p) = 5,56\%$ a $P(L \mid E_s) = 40,93\%$. El Hospital Base San José de Osorno ($n_h = 47$) tiene solo un $6,4\%$ de episodios secundarios ($93,6\%$ son cirugías oftalmológicas programadas $E_p$ de $1$ día de estadía), razón por la cual observa $0$ fallecidos y apenas $2$ estancias largas. Por el contrario, en San José ($n_h = 40$), Curicó ($n_h = 29$) y La Florida ($n_h = 17$), el $100,0\%$ de sus episodios tiene H25 como diagnóstico secundario en internaciones médicas o quirúrgicas complejas. Asumir una única $p = 4,11\%$ o $p_L = 31,11\%$ mezcla subpoblaciones clínicamente heterogéneas.
- Supuesto 2: Ensayos estadísticamente independientes. SE CUMPLE SOLO PARCIALMENTE. Si el experimento consiste en extraer fichas al azar con reposición desde una base estática, los sorteos son independientes por construcción. Sin embargo, los episodios reales dentro de un establecimiento no son independientes entre sí: existen reingresos de una misma persona (por ejemplo, cirugía del segundo ojo o rehospitalizaciones por descompensación de su patología de base) y factores institucionales compartidos (disponibilidad de camas, protocolos locales de alta y perfil socioeconómico territorial) que inducen correlación positiva intra-hospitalaria.

3 Resultado 03: Distribución Hipergeométrica (auditar sin reposición)

Un auditor revisa $n = 10$ episodios sin reposición dentro de un establecimiento. Seleccionamos los dos establecimientos extremos exigidos por la pauta entre los $9$ que cuentan con al menos $10$ episodios de la cohorte $E$:

Establecimiento de mayor volumen en la cohorte: Hospital Clínico Regional Dr. Guillermo Grant Benavente (Concepción) (118100), con $N_h = 83$ episodios.
Establecimiento de menor volumen entre los que tienen al menos 10 episodios: Hospital Clínico Regional (Valdivia) (122100), con $N_h = 13$ episodios.

Como en ambos recintos el número de fallecimientos observados es cero ($K_F = 0$), modelamos como evento principal de auditoría el número de episodios con estancia larga ($K_h = K_L$, estadía $> 7$ días), tal como contempla explícitamente la guía (“fallecimientos o estancias largas entre ellos”). Adicionalmente, para no dejar fuera el análisis sobre mortalidad, incluimos una tabla resumen con los dos establecimientos extremos con $N_h \ge 10$ que sí registran fallecimientos ($K_F > 0$): Complejo Hospitalario San José ($N_h = 40, K_F = 3$) y Hospital San Juan de Dios de Santiago ($N_h = 16, K_F = 1$).

n_audit <- 10
k_hyp <- 0:n_audit

# 1. Estancia larga en los dos establecimientos extremos absolutos (N_h = 83 y N_h = 13)
N_mayor <- 83; K_mayor <- 19; p_mayor <- K_mayor / N_mayor
N_menor <- 13; K_menor <- 8;  p_menor <- K_menor / N_menor

df_hyper_comp <- tibble(
  k = k_hyp,
  H_mayor = dhyper(k_hyp, m = K_mayor, n = N_mayor - K_mayor, k = n_audit),
  B_mayor = dbinom(k_hyp, size = n_audit, prob = p_mayor),
  Diff_mayor = abs(H_mayor - B_mayor),
  H_menor = dhyper(k_hyp, m = K_menor, n = N_menor - K_menor, k = n_audit),
  B_menor = dbinom(k_hyp, size = n_audit, prob = p_menor),
  Diff_menor = abs(H_menor - B_menor)
)

kable(
  df_hyper_comp %>%
    transmute(
      `k (En 10 auditados)` = k,
      `P(H = k) dhyper [Grant Benavente]` = sprintf("%.6f (%.2f%%)", H_mayor, H_mayor * 100),
      `P(B = k) dbinom [Grant Benavente]` = sprintf("%.6f (%.2f%%)", B_mayor, B_mayor * 100),
      `|P(H) - P(B)| [Grant Benavente]` = sprintf("%.6f (%.2f pp)", Diff_mayor, Diff_mayor * 100),
      `P(H = k) dhyper [Valdivia]` = sprintf("%.6f (%.2f%%)", H_menor, H_menor * 100),
      `P(B = k) dbinom [Valdivia]` = sprintf("%.6f (%.2f%%)", B_menor, B_menor * 100),
      `|P(H) - P(B)| [Valdivia]` = sprintf("%.6f (%.2f pp)", Diff_menor, Diff_menor * 100)
    ),
  align = c("c", "r", "r", "r", "r", "r", "r"),
  caption = "Tabla 10: Función de probabilidad Hipergeométrica (dhyper, sin reposición) vs Binomial (dbinom, con reposición) para n = 10 episodios auditados en el hospital de mayor volumen (Grant Benavente, N_h = 83, K_L = 19) y el de menor volumen con N_h ≥ 10 (Valdivia, N_h = 13, K_L = 8)."
) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = TRUE)

Tabla 10: Función de probabilidad Hipergeométrica (dhyper, sin reposición) vs Binomial (dbinom, con reposición) para n = 10 episodios auditados en el hospital de mayor volumen (Grant Benavente, N_h = 83, K_L = 19) y el de menor volumen con N_h ≥ 10 (Valdivia, N_h = 13, K_L = 8).
k (En 10 auditados)	P(H = k) dhyper [Grant Benavente]	P(B = k) dbinom [Grant Benavente]	\|P(H) - P(B)\| [Grant Benavente]	P(H = k) dhyper [Valdivia]	P(B = k) dbinom [Valdivia]	\|P(H) - P(B)\| [Valdivia]
0	0.062275 (6.23%)	0.074305 (7.43%)	0.012030 (1.20 pp)	0.000000 (0.00%)	0.000071 (0.01%)	0.000071 (0.01 pp)
1	0.215131 (21.51%)	0.220593 (22.06%)	0.005462 (0.55 pp)	0.000000 (0.00%)	0.001133 (0.11%)	0.001133 (0.11 pp)
2	0.311172 (31.12%)	0.294699 (29.47%)	0.016473 (1.65 pp)	0.000000 (0.00%)	0.008161 (0.82%)	0.008161 (0.82 pp)
3	0.247482 (24.75%)	0.233303 (23.33%)	0.014179 (1.42 pp)	0.000000 (0.00%)	0.034818 (3.48%)	0.034818 (3.48 pp)
4	0.119474 (11.95%)	0.121208 (12.12%)	0.001734 (0.17 pp)	0.000000 (0.00%)	0.097491 (9.75%)	0.097491 (9.75 pp)
5	0.036450 (3.64%)	0.043180 (4.32%)	0.006731 (0.67 pp)	0.195804 (19.58%)	0.187183 (18.72%)	0.008621 (0.86 pp)
6	0.007087 (0.71%)	0.010683 (1.07%)	0.003595 (0.36 pp)	0.489510 (48.95%)	0.249578 (24.96%)	0.239933 (23.99 pp)
7	0.000863 (0.09%)	0.001812 (0.18%)	0.000949 (0.09 pp)	0.279720 (27.97%)	0.228185 (22.82%)	0.051535 (5.15 pp)
8	0.000063 (0.01%)	0.000202 (0.02%)	0.000139 (0.01 pp)	0.034965 (3.50%)	0.136911 (13.69%)	0.101946 (10.19 pp)
9	0.000002 (0.00%)	0.000013 (0.00%)	0.000011 (0.00 pp)	0.000000 (0.00%)	0.048680 (4.87%)	0.048680 (4.87 pp)
10	0.000000 (0.00%)	0.000000 (0.00%)	0.000000 (0.00 pp)	0.000000 (0.00%)	0.007789 (0.78%)	0.007789 (0.78 pp)

# Tabla resumen de parámetros, fracción muestral, diferencia máxima y varianzas (incluyendo también Fallecimientos en N_h=40 y N_h=16)
resumen_hyper_modelos <- function(nombre, evento, N_h, K_h, n = 10) {
  p_h <- K_h / N_h
  prob_H <- dhyper(0:n, m = K_h, n = N_h - K_h, k = n)
  prob_B <- dbinom(0:n, size = n, prob = p_h)
  diffs <- abs(prob_H - prob_B)
  k_max <- (0:n)[which.max(diffs)]
  var_H <- n * p_h * (1 - p_h) * ((N_h - n) / (N_h - 1))
  var_B <- n * p_h * (1 - p_h)
  fpc <- (N_h - n) / (N_h - 1)
  tibble(
    Establecimiento = nombre,
    Evento = evento,
    `N_h` = N_h,
    `K_h` = K_h,
    `p_h = K_h / N_h` = sprintf("%.4f (%.2f%%)", p_h, p_h * 100),
    `Fracción muestral n / N_h` = sprintf("%d / %d = %.2f%%", n, N_h, (n / N_h) * 100),
    `Factor finitud (N_h-n)/(N_h-1)` = sprintf("%.4f", fpc),
    `Máx |dhyper - dbinom|` = sprintf("%.4f (%.2f pp en k = %d)", max(diffs), max(diffs) * 100, k_max),
    `Var[H] (Hipergeométrica)` = sprintf("%.4f", var_H),
    `Var[B] (Binomial)` = sprintf("%.4f", var_B)
  )
}

tabla_resumen_hyper <- bind_rows(
  resumen_hyper_modelos("1. Hosp. Dr. Guillermo Grant Benavente (Mayor N_h)", "Estancia larga (L)", 83, 19, 10),
  resumen_hyper_modelos("2. Hosp. Clínico Regional de Valdivia (Menor N_h ≥ 10)", "Estancia larga (L)", 13, 8, 10),
  resumen_hyper_modelos("3. Complejo Hosp. San José (Mayor N_h con K_F > 0)", "Fallecimiento (F)", 40, 3, 10),
  resumen_hyper_modelos("4. Hosp. San Juan de Dios Stgo. (Menor N_h ≥ 10 con K_F > 0)", "Fallecimiento (F)", 16, 1, 10)
)

kable(
  tabla_resumen_hyper,
  align = c("l", "l", "r", "r", "r", "r", "r", "r", "r", "r"),
  caption = "Tabla 11: Resumen de parámetros N_h, K_h, fracción muestral n/N_h, máxima discrepancia puntual y varianzas Hipergeométrica vs Binomial."
) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = TRUE)

Tabla 11: Resumen de parámetros N_h, K_h, fracción muestral n/N_h, máxima discrepancia puntual y varianzas Hipergeométrica vs Binomial.
Establecimiento	Evento	N_h	K_h	p_h = K_h / N_h	Fracción muestral n / N_h	Factor finitud (N_h-n)/(N_h-1)	Máx \|dhyper - dbinom\|	Var[H] (Hipergeométrica)	Var[B] (Binomial)
Hosp. Dr. Guillermo Grant Benavente (Mayor N_h)	Estancia larga (L)	83	19	0.2289 (22.89%)	10 / 83 = 12.05%	0.8902	0.0165 (1.65 pp en k = 2)	1.5714	1.7651
Hosp. Clínico Regional de Valdivia (Menor N_h ≥ 10)	Estancia larga (L)	13	8	0.6154 (61.54%)	10 / 13 = 76.92%	0.2500	0.2399 (23.99 pp en k = 6)	0.5917	2.3669
Complejo Hosp. San José (Mayor N_h con K_F > 0)	Fallecimiento (F)	40	3	0.0750 (7.50%)	10 / 40 = 25.00%	0.7692	0.0685 (6.85 pp en k = 1)	0.5337	0.6938
Hosp. San Juan de Dios Stgo. (Menor N_h ≥ 10 con K_F > 0)	Fallecimiento (F)	16	1	0.0625 (6.25%)	10 / 16 = 62.50%	0.4000	0.2754 (27.54 pp en k = 1)	0.2344	0.5859

Interpretación de las Tablas 10 y 11:

En el Hospital Dr. Guillermo Grant Benavente ($N_h = 83, K_h = 19, p_h = 22,89\%$), la fracción muestral auditada es pequeña ($n / N_h = 10 / 83 = 12,05\%$). Por ello, la probabilidad Hipergeométrica y su aproximación Binomial son muy similares en todo el soporte: la mayor discrepancia puntual ocurre en $k = 2$ y es de apenas $0,0165$ ($1,65$ puntos porcentuales: $31,42\%$ sin reposición frente a $29,77\%$ con reposición), y la varianza Hipergeométrica ($\text{Var}[H] = 1,5714$) es solo un $11,0\%$ menor que la Binomial ($\text{Var}[B] = 1,7651$).
En el Hospital Clínico Regional de Valdivia ($N_h = 13, K_h = 8, p_h = 61,54\%$), la auditoría sin reposición de $10$ episodios extrae más de tres cuartas partes del total del establecimiento ($n / N_h = 10 / 13 = 76,92\%$). Aquí el supuesto binomial colapsa por completo: como en Valdivia solo existen $N_h - K_h = 5$ episodios sin estancia larga, al extraer $10$ fichas sin reposición es matemáticamente imposible encontrar menos de $5$ estancias largas ($P(H \le 4) = 0$), mientras que la Binomial asigna un $13,48\%$ de probabilidad acumulada a $k \in \{0, 1, 2, 3, 4\}$. La diferencia máxima alcanza $0,2399$ ($23,99$ puntos porcentuales en $k = 6$: $48,95\%$ en la Hipergeométrica frente a $24,96\%$ en la Binomial), y la varianza real sin reposición ($\text{Var}[H] = 0,5917$) es cuatro veces menor que la varianza Binomial ($\text{Var}[B] = 2,3669$).
Exactamente el mismo fenómeno se verifica sobre los fallecimientos ($F$) al comparar el Complejo Hospitalario San José ($N_h = 40, K_F = 3, n/N_h = 25,00\%$, diferencia máxima de $6,85$ pp y $\text{Var}[H] = 0,5337$ vs $\text{Var}[B] = 0,6938$) contra el Hospital San Juan de Dios de Santiago ($N_h = 16, K_F = 1, n/N_h = 62,50\%$, diferencia máxima de $27,54$ pp en $k = 1$, donde $P(H = 1) = 62,50\%$ exactamente frente a $P(B = 1) = 34,96\%$, y además $P(H \ge 2) = 0$ porque solo existe $K_F = 1$ fallecido en todo el recinto).

df_plot_hyper <- bind_rows(
  tibble(k = k_hyp, Prob = df_hyper_comp$H_mayor, Modelo = "Hipergeométrica (Sin reposición)", Recinto = "1. Hosp. Grant Benavente (N_h = 83, K_L = 19, n/N_h = 12,0%)"),
  tibble(k = k_hyp, Prob = df_hyper_comp$B_mayor, Modelo = "Binomial (Con reposición)", Recinto = "1. Hosp. Grant Benavente (N_h = 83, K_L = 19, n/N_h = 12,0%)"),
  tibble(k = k_hyp, Prob = df_hyper_comp$H_menor, Modelo = "Hipergeométrica (Sin reposición)", Recinto = "2. Hosp. Regional Valdivia (N_h = 13, K_L = 8, n/N_h = 76,9%)"),
  tibble(k = k_hyp, Prob = df_hyper_comp$B_menor, Modelo = "Binomial (Con reposición)", Recinto = "2. Hosp. Regional Valdivia (N_h = 13, K_L = 8, n/N_h = 76,9%)")
)

ggplot(df_plot_hyper, aes(x = k, y = Prob, fill = Modelo)) +
  geom_col(position = "dodge", width = 0.72, color = "white") +
  facet_wrap(~ Recinto, scales = "free_y") +
  scale_fill_manual(values = c("Hipergeométrica (Sin reposición)" = "#1f618d", "Binomial (Con reposición)" = "#e67e22")) +
  scale_x_continuous(breaks = 0:10) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  theme_minimal(base_size = 11) +
  theme(legend.position = "top") +
  labs(
    title = "Figura 4: Comparación Hipergeométrica (dhyper) vs Binomial (dbinom) al auditar n = 10 episodios",
    subtitle = "A menor fracción muestral n/N_h (izquierda), ambas distribuciones convergen; a alta fracción muestral (derecha), divergen",
    x = "Número de episodios con estancia larga entre los 10 auditados (k)",
    y = "Probabilidad puntual P(Variable = k)",
    fill = "Modelo de probabilidad:"
  )

Interpretación de la Figura 4 y explicación mediante la fracción muestral $n / N_h$:

La relación matemática entre la varianza Hipergeométrica (muestreo sin reposición) y la varianza Binomial (muestreo con reposición) está gobernada por el factor de corrección por población finita (FPC): \[\text{Var}[H] = n \cdot p_h (1 - p_h) \cdot \left(\frac{N_h - n}{N_h - 1}\right) = \text{Var}[B] \cdot \left(1 - \frac{n - 1}{N_h - 1}\right) \approx \text{Var}[B] \cdot \left(1 - \frac{n}{N_h}\right)\] Cuando la fracción muestral $n / N_h$ es pequeña (como en el Hospital Guillermo Grant Benavente, donde $n / N_h = 10 / 83 = 12,05\%$ y el factor de corrección es $73/82 = 0,8902$), extraer un episodio sin devolverlo altera muy poco la composición restante de la urna hospitalaria (la probabilidad condicional del segundo ensayo pasa de $19/83 = 22,89\%$ a $18/82 = 21,95\%$ o $19/82 = 23,17\%$), por lo que las barras azules y naranjas del panel izquierdo casi se superponen. En cambio, cuando la fracción muestral $n / N_h$ es grande (Hospital de Valdivia, donde $n / N_h = 10 / 13 = 76,92\%$ y el factor de corrección cae a $3/12 = 0,2500$), cada extracción sin reposición modifica drásticamente las probabilidades condicionales subsiguientes y agota los casos disponibles, concentrando la distribución Hipergeométrica estrechamente alrededor de su media ($E[H] = 6,15$) y haciendo inapropiada la aproximación Binomial.

4 Resultado 04: Geométrica y Binomial Negativa (revisar fichas hasta encontrar casos)

4.1 Declaración de la convención de R y ajuste del marcador de registro $M$

Declaración explícita de la convención de R y verificación manual:

En R, las funciones nativas dgeom(x, prob), pgeom(q, prob), qgeom(p, prob) y rgeom(n, prob) parametrizan la distribución Geométrica contando el número de fracasos $x \in \{0, 1, 2, \dots\}$ antes del primer éxito, con función de probabilidad $P(X_{\text{R}} = x) = (1 - p)^x \cdot p$. Dado que en nuestro problema administrativo la variable aleatoria $T$ cuenta el número total de ensayos (fichas revisadas) hasta encontrar el primer episodio con Catarata Senil ($k \in \{1, 2, 3, \dots\}$), se cumple la identidad $k = x + 1$, es decir, $x = k - 1$ fracasos previos. Por lo tanto:

Probabilidad puntual: $P(T = k) = (1 - p)^{k-1} p = \texttt{dgeom(k - 1, prob = p)}$ para $k = 1, 2, \dots$
Probabilidad acumulada: $P(T \le k) = 1 - (1 - p)^k = \texttt{pgeom(k - 1, prob = p)}$ y $P(T > k) = (1 - p)^k = \texttt{1 - pgeom(k - 1, prob = p)}$.
Cuantil de orden $\alpha$: $\min\{k \in \mathbb{N} : P(T \le k) \ge \alpha\} = \texttt{qgeom(alpha, prob = p) + 1}$.
Simulación: $\texttt{rgeom(B, prob = p) + 1}$.

De manera análoga, dnbinom(x, size = r, prob = p) y pnbinom(q, size = r, prob = p) cuentan el número de fracasos $x = k - r$ antes de alcanzar el $r$-ésimo éxito. Para $T_5$ (número total de fichas revisadas hasta encontrar $r = 5$ episodios con la enfermedad, con soporte $k \in \{5, 6, 7, \dots\}$):

Probabilidad puntual: $P(T_5 = k) = \binom{k-1}{5-1} p^5 (1-p)^{k-5} = \texttt{dnbinom(k - 5, size = 5, prob = p)}$.
Probabilidad acumulada: $P(T_5 \le 100) = \texttt{pnbinom(100 - 5, size = 5, prob = p)}$.

Justificación y cálculo del ajuste del marcador $M$:

En la Actividad 2 se definió el marcador administrativo $M_{\text{orig}}$ como “episodio con edad $\ge 65$ años e ingreso programado ($\overline{U}$)”, el cual sobre $\Omega$ ($N = 873.718$) marcó $57.303$ fichas con $110$ casos de H25, arrojando un valor predictivo positivo casi nulo de $\text{VPP}_{\text{orig}} = P(E \mid M_{\text{orig}}) = 110 / 57.303 = 0,001920$ ($0,1920\%$). Siguiendo la instrucción de la guía (“Si el valor predictivo positivo de su marcador fue cero o casi cero, ajuste el marcador antes de seguir y explique el cambio”), ajustamos el marcador $M$ combinando el grupo etario adulto mayor ($\text{edad} \ge 65$ años) con el registro administrativo de atención o egreso por la especialidad/servicio de Oftalmología (ESPECIALIDAD_MEDICA o SERVICIOALTA), variable de registro clínico-administrativo que no utiliza ningún campo diagnóstico CIE-10 ni el código GRD (no es circular). Sobre $\Omega$, el marcador ajustado $M$ selecciona $n(M) = 962$ episodios, entre los cuales se encuentran $n(E \cap M) = 99$ episodios de Catarata Senil (H25), elevando el valor predictivo positivo a $\text{VPP} = P(E \mid M) = 99 / 962 = 0,102911$ ($10,2911\%$). Para máxima completitud y trazabilidad, presentamos en las tablas siguientes los resultados con el marcador ajustado $M$ ($\text{VPP} = 10,2911\%$), con el marcador original $M_{\text{orig}}$ de la Actividad 2 ($\text{VPP}_{\text{orig}} = 0,1920\%$) y sin usar marcador (revisión al azar sobre $\Omega$, con $p = P(E) = 389 / 873.718 = 0,0445\%$).

vpp_M <- sum(omega$E & omega$M) / sum(omega$M)
vpp_M_orig <- sum(omega$E & omega$M_orig) / sum(omega$M_orig)
p_azar <- P_E

# Verificación manual en k = 1 y k = 2 para el marcador ajustado M
p_man_k1 <- (1 - vpp_M)^(1 - 1) * vpp_M
p_r_k1 <- dgeom(1 - 1, prob = vpp_M)
p_man_k2 <- (1 - vpp_M)^(2 - 1) * vpp_M
p_r_k2 <- dgeom(2 - 1, prob = vpp_M)

# Tabla de la función de probabilidad de T para k = 1 a 50
k_geom <- 1:50
df_geom_50 <- tibble(
  k = k_geom,
  p_T_M = dgeom(k_geom - 1, prob = vpp_M),
  F_T_M = pgeom(k_geom - 1, prob = vpp_M),
  p_T_orig = dgeom(k_geom - 1, prob = vpp_M_orig),
  F_T_orig = pgeom(k_geom - 1, prob = vpp_M_orig),
  p_T_azar = dgeom(k_geom - 1, prob = p_azar),
  F_T_azar = pgeom(k_geom - 1, prob = p_azar)
)

# Mostrar k seleccionados representativos de 1 a 50 en tabla compacta y gráfica completa de k = 1 a 50
k_mostrar_geom <- c(1:10, 12, 15, 20, 22, 25, 30, 35, 40, 45, 50)
kable(
  df_geom_50 %>%
    filter(k %in% k_mostrar_geom) %>%
    transmute(
      `k (Fichas revisadas)` = k,
      `P(T = k) [Con M ajustado: VPP = 10,29%]` = sprintf("%.6f (%.2f%%)", p_T_M, p_T_M * 100),
      `P(T ≤ k) [Con M ajustado]` = sprintf("%.6f (%.2f%%)", F_T_M, F_T_M * 100),
      `P(T = k) [Con M_orig Act.2: VPP = 0,19%]` = sprintf("%.6f (%.4f%%)", p_T_orig, p_T_orig * 100),
      `P(T ≤ k) [Con M_orig Act.2]` = sprintf("%.6f (%.2f%%)", F_T_orig, F_T_orig * 100),
      `P(T = k) [Sin marcador: P(E) = 0,0445%]` = sprintf("%.6f (%.4f%%)", p_T_azar, p_T_azar * 100),
      `P(T ≤ k) [Sin marcador]` = sprintf("%.6f (%.2f%%)", F_T_azar, F_T_azar * 100)
    ),
  align = c("c", "r", "r", "r", "r", "r", "r"),
  caption = "Tabla 12: Función de probabilidad puntual P(T = k) = dgeom(k - 1, p) y acumulada P(T ≤ k) = pgeom(k - 1, p) para k ∈ [1, 50] con marcador ajustado M, marcador original M_orig y revisión al azar sobre Ω."
) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = TRUE)

Tabla 12: Función de probabilidad puntual P(T = k) = dgeom(k - 1, p) y acumulada P(T ≤ k) = pgeom(k - 1, p) para k ∈ [1, 50] con marcador ajustado M, marcador original M_orig y revisión al azar sobre Ω.
k (Fichas revisadas)	P(T = k) [Con M ajustado: VPP = 10,29%]	P(T ≤ k) [Con M ajustado]	P(T = k) [Con M_orig Act.2: VPP = 0,19%]	P(T ≤ k) [Con M_orig Act.2]	P(T = k) [Sin marcador: P(E) = 0,0445%]	P(T ≤ k) [Sin marcador]
1	0.102911 (10.29%)	0.102911 (10.29%)	0.001920 (0.1920%)	0.001920 (0.19%)	0.000445 (0.0445%)	0.000445 (0.04%)
2	0.092320 (9.23%)	0.195231 (19.52%)	0.001916 (0.1916%)	0.003836 (0.38%)	0.000445 (0.0445%)	0.000890 (0.09%)
3	0.082819 (8.28%)	0.278050 (27.80%)	0.001912 (0.1912%)	0.005748 (0.57%)	0.000445 (0.0445%)	0.001335 (0.13%)
4	0.074296 (7.43%)	0.352346 (35.23%)	0.001909 (0.1909%)	0.007656 (0.77%)	0.000445 (0.0445%)	0.001780 (0.18%)
5	0.066650 (6.67%)	0.418997 (41.90%)	0.001905 (0.1905%)	0.009561 (0.96%)	0.000444 (0.0444%)	0.002224 (0.22%)
6	0.059791 (5.98%)	0.478788 (47.88%)	0.001901 (0.1901%)	0.011463 (1.15%)	0.000444 (0.0444%)	0.002668 (0.27%)
7	0.053638 (5.36%)	0.532426 (53.24%)	0.001898 (0.1898%)	0.013360 (1.34%)	0.000444 (0.0444%)	0.003112 (0.31%)
8	0.048118 (4.81%)	0.580545 (58.05%)	0.001894 (0.1894%)	0.015254 (1.53%)	0.000444 (0.0444%)	0.003556 (0.36%)
9	0.043166 (4.32%)	0.623711 (62.37%)	0.001890 (0.1890%)	0.017145 (1.71%)	0.000444 (0.0444%)	0.004000 (0.40%)
10	0.038724 (3.87%)	0.662435 (66.24%)	0.001887 (0.1887%)	0.019031 (1.90%)	0.000443 (0.0443%)	0.004443 (0.44%)
12	0.031164 (3.12%)	0.728338 (72.83%)	0.001879 (0.1879%)	0.022794 (2.28%)	0.000443 (0.0443%)	0.005330 (0.53%)
15	0.022499 (2.25%)	0.803874 (80.39%)	0.001869 (0.1869%)	0.028411 (2.84%)	0.000442 (0.0442%)	0.006658 (0.67%)
20	0.013072 (1.31%)	0.886050 (88.60%)	0.001851 (0.1851%)	0.037700 (3.77%)	0.000441 (0.0441%)	0.008867 (0.89%)
22	0.010520 (1.05%)	0.908296 (90.83%)	0.001844 (0.1844%)	0.041391 (4.14%)	0.000441 (0.0441%)	0.009749 (0.97%)
25	0.007595 (0.76%)	0.933795 (93.38%)	0.001833 (0.1833%)	0.046901 (4.69%)	0.000440 (0.0440%)	0.011071 (1.11%)
30	0.004413 (0.44%)	0.961534 (96.15%)	0.001816 (0.1816%)	0.056014 (5.60%)	0.000440 (0.0440%)	0.013271 (1.33%)
35	0.002564 (0.26%)	0.977651 (97.77%)	0.001798 (0.1798%)	0.065040 (6.50%)	0.000439 (0.0439%)	0.015465 (1.55%)
40	0.001490 (0.15%)	0.987015 (98.70%)	0.001781 (0.1781%)	0.073979 (7.40%)	0.000438 (0.0438%)	0.017655 (1.77%)
45	0.000865 (0.09%)	0.992456 (99.25%)	0.001764 (0.1764%)	0.082833 (8.28%)	0.000437 (0.0437%)	0.019840 (1.98%)
50	0.000503 (0.05%)	0.995617 (99.56%)	0.001747 (0.1747%)	0.091603 (9.16%)	0.000436 (0.0436%)	0.022020 (2.20%)

Interpretación de la Tabla 12 y verificación a mano:

Verificación manual de la convención en R: Para el marcador ajustado $M$ ($p = \text{VPP} = 99 / 962 = 0,1029106$):
- En $k = 1$ ensayo ($0$ fracasos previos): a mano, $P(T = 1) = (1 - 0,1029106)^0 \times 0,1029106 = 0,102911$ ($10,29\%$), idéntico a dgeom(1 - 1, prob = vpp_M) $= 0.102911$.
- En $k = 2$ ensayos ($1$ fracaso previo): a mano, $P(T = 2) = (1 - 0,1029106)^1 \times 0,1029106 = 0,8970894 \times 0,1029106 = 0,092320$ ($9,23\%$), idéntico a dgeom(2 - 1, prob = vpp_M) $= 0.092320$.
Comportamiento en $k = 1 \dots 50$: Con el marcador ajustado $M$, la probabilidad puntual decrece geométricamente desde $10,29\%$ en la primera ficha hasta $0,05\%$ en la ficha $k = 50$, acumulando un $99,56\%$ de probabilidad de encontrar el primer caso dentro de las primeras $50$ fichas ($P(T \le 50) = 99,56\%$). En cambio, sin marcador ($p = P(E) = 0,0445\%$), revisar $50$ fichas al azar de $\Omega$ apenas acumula un $2,20\%$ de probabilidad de éxito ($P(T_{\text{azar}} \le 50) = 2,20\%$).

df_geom_plot <- df_geom_50 %>%
  select(
    k,
    `1. Con Marcador Ajustado M (VPP = 10,29%)` = p_T_M,
    `2. Con Marcador Original Act.2 (VPP = 0,192%)` = p_T_orig,
    `3. Sin Marcador - Azar sobre Ω (P(E) = 0,0445%)` = p_T_azar
  ) %>%
  pivot_longer(cols = -k, names_to = "Estrategia", values_to = "Probabilidad")

ggplot(df_geom_plot, aes(x = k, y = Probabilidad, color = Estrategia, fill = Estrategia)) +
  geom_col(data = df_geom_plot %>% filter(grepl("^1", Estrategia)), alpha = 0.35, width = 0.7) +
  geom_line(linewidth = 1.1) +
  geom_point(size = 1.5) +
  scale_color_manual(values = c(
    "1. Con Marcador Ajustado M (VPP = 10,29%)" = "#1e8449",
    "2. Con Marcador Original Act.2 (VPP = 0,192%)" = "#d35400",
    "3. Sin Marcador - Azar sobre Ω (P(E) = 0,0445%)" = "#7f8c8d"
  )) +
  scale_fill_manual(values = c(
    "1. Con Marcador Ajustado M (VPP = 10,29%)" = "#1e8449",
    "2. Con Marcador Original Act.2 (VPP = 0,192%)" = "#d35400",
    "3. Sin Marcador - Azar sobre Ω (P(E) = 0,0445%)" = "#7f8c8d"
  )) +
  scale_x_continuous(breaks = seq(0, 50, by = 5)) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  theme_minimal(base_size = 11) +
  theme(legend.position = "top") +
  labs(
    title = "Figura 5: Función de probabilidad P(T = k) para k = 1 a 50 fichas revisadas",
    subtitle = "El marcador ajustado M concentra la masa en pocas fichas (E[T] = 9,72), frente a la dispersión casi plana sin marcador",
    x = "Número de fichas revisadas hasta encontrar el primer episodio con H25 (k)",
    y = "Probabilidad puntual P(T = k)",
    color = "Estrategia de búsqueda:",
    fill = "Estrategia de búsqueda:"
  )

Interpretación de la Figura 5: La curva verde muestra cómo el marcador ajustado $M$ ($\text{VPP} = 10,29\%$) concentra fuertemente la probabilidad de hallazgo en las primeras dos decenas de fichas revisadas, mientras que las estrategias con baja proporción de acierto (naranja y gris) tienen una probabilidad puntual casi nula y constante a lo largo de las primeras $50$ revisiones.

4.2 Comparación de indicadores de $T$ y $T_5$ y verificación por remuestreo sin reposición

Calculamos para las tres estrategias de auditoría los valores exigidos de $E[T]$, $\text{Var}[T]$, mediana de $T$, $P(T > 50)$, y para la Binomial Negativa $T_5 \sim \text{Binomial Negativa}(r = 5, p)$ la esperanza $E[T_5] = 5 / p$, la varianza $\text{Var}[T_5] = 5(1-p)/p^2$ y la probabilidad $P(T_5 \le 100)$. Finalmente, verificamos con los datos reales desordenando al azar $1.000$ veces (set.seed(2026)) los episodios marcados por $M$ y registrando la posición del primer episodio con Catarata Senil (H25).

set.seed(2026)
B_geom <- 1000

# Remuestreo sin reposición sobre los 962 episodios del marcador ajustado M
episodios_M_adj <- omega$E[which(omega$M)]
pos_primero_adj <- replicate(B_geom, which(sample(episodios_M_adj) == TRUE)[1])
media_emp_T_adj <- mean(pos_primero_adj)
var_emp_T_adj <- var(pos_primero_adj)
esperanza_exacta_sin_rep_adj <- (length(episodios_M_adj) + 1) / (sum(episodios_M_adj) + 1)

# Remuestreo sin reposición sobre los 57.303 episodios del marcador original M_orig
set.seed(2026)
episodios_M_orig <- omega$E[which(omega$M_orig)]
pos_primero_orig <- replicate(B_geom, which(sample(episodios_M_orig) == TRUE)[1])
media_emp_T_orig <- mean(pos_primero_orig)
var_emp_T_orig <- var(pos_primero_orig)
esperanza_exacta_sin_rep_orig <- (length(episodios_M_orig) + 1) / (sum(episodios_M_orig) + 1)

calc_metricas_geom_bn <- function(p, n_total, n_exito) {
  ET <- 1 / p
  VarT <- (1 - p) / (p^2)
  medT <- qgeom(0.50, prob = p) + 1
  PT_gt_50 <- 1 - pgeom(50 - 1, prob = p)
  ET5 <- 5 / p
  VarT5 <- 5 * (1 - p) / (p^2)
  PT5_le_100 <- pnbinom(100 - 5, size = 5, prob = p)
  k90 <- qgeom(0.90, prob = p) + 1
  c(
    sprintf("%d / %s = %.6f (%.4f%%)", n_exito, format(n_total, big.mark = "."), p, p * 100),
    sprintf("%.2f fichas", ET),
    sprintf("%.2f fichas²", VarT),
    sprintf("%d fichas (P(T ≤ %d) = %.2f%%)", medT, medT, pgeom(medT - 1, prob = p) * 100),
    sprintf("%.6f (%.2f%%)", PT_gt_50, PT_gt_50 * 100),
    sprintf("%.2f fichas", ET5),
    sprintf("%.2f fichas²", VarT5),
    sprintf("%.6f (%.4f%%)", PT5_le_100, PT5_le_100 * 100),
    sprintf("%d fichas (P(T ≤ %d) = %.2f%%)", k90, k90, pgeom(k90 - 1, prob = p) * 100)
  )
}

tabla_comp_T_T5 <- tibble(
  Metrica = c(
    "Parámetro de éxito p (origen y denominador)",
    "Esperanza E[T] = 1 / p (fichas hasta el 1° caso)",
    "Varianza Var[T] = (1 - p) / p²",
    "Mediana de T: qgeom(0.50, p) + 1",
    "Probabilidad de cola P(T > 50) = (1 - p)⁵⁰",
    "Esperanza E[T₅] = 5 / p (fichas hasta el 5° caso)",
    "Varianza Var[T₅] = 5(1 - p) / p²",
    "Probabilidad P(T₅ ≤ 100) = pnbinom(95, size = 5, prob = p)",
    "Menor k con P(T ≤ k) ≥ 0,90: qgeom(0.90, p) + 1"
  ),
  Con_M_Ajustado = calc_metricas_geom_bn(vpp_M, sum(omega$M), sum(omega$E & omega$M)),
  Con_M_Original = calc_metricas_geom_bn(vpp_M_orig, sum(omega$M_orig), sum(omega$E & omega$M_orig)),
  Sin_Marcador_Azar = calc_metricas_geom_bn(p_azar, N_omega, n_E)
)

kable(
  tabla_comp_T_T5,
  col.names = c(
    "Indicador de los modelos Geométrico (T) y Binomial Negativo (T₅)",
    "Con Marcador Ajustado M (Oftalmología y ≥ 65 años)",
    "Con Marcador Original Act. 2 M_orig (≥ 65 años y No Urgencia)",
    "Sin Marcador (Azar sobre Ω)"
  ),
  align = c("l", "r", "r", "r"),
  caption = "Tabla 13: Comparación completa de T ~ Geométrica(p) y T₅ ~ Binomial Negativa(5, p) usando el marcador ajustado M, el marcador original de la Actividad 2 y la búsqueda al azar sobre Ω."
) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = TRUE)

Tabla 13: Comparación completa de T ~ Geométrica(p) y T₅ ~ Binomial Negativa(5, p) usando el marcador ajustado M, el marcador original de la Actividad 2 y la búsqueda al azar sobre Ω.
Indicador de los modelos Geométrico (T) y Binomial Negativo (T₅)	Con Marcador Ajustado M (Oftalmología y ≥ 65 años)	Con Marcador Original Act. 2 M_orig (≥ 65 años y No Urgencia)	Sin Marcador (Azar sobre Ω)
Parámetro de éxito p (origen y denominador)	99 / 962 = 0.102911 (10.2911%)	110 / 57.303 = 0.001920 (0.1920%)	389 / 873.718 = 0.000445 (0.0445%)
Esperanza E[T] = 1 / p (fichas hasta el 1° caso)	9.72 fichas	520.94 fichas	2246.06 fichas
Varianza Var[T] = (1 - p) / p²	84.71 fichas²	270853.76 fichas²	5042547.08 fichas²
Mediana de T: qgeom(0.50, p) + 1	7 fichas (P(T ≤ 7) = 53.24%)	361 fichas (P(T ≤ 361) = 50.03%)	1557 fichas (P(T ≤ 1557) = 50.01%)
Probabilidad de cola P(T > 50) = (1 - p)⁵⁰	0.004383 (0.44%)	0.908397 (90.84%)	0.977980 (97.80%)
Esperanza E[T₅] = 5 / p (fichas hasta el 5° caso)	48.59 fichas	2604.68 fichas	11230.31 fichas
Varianza Var[T₅] = 5(1 - p) / p²	423.53 fichas²	1354268.79 fichas²	25212735.42 fichas²
Probabilidad P(T₅ ≤ 100) = pnbinom(95, size = 5, prob = p)	0.980773 (98.0773%)	0.000002 (0.0002%)	0.000000 (0.0000%)
Menor k con P(T ≤ k) ≥ 0,90: qgeom(0.90, p) + 1	22 fichas (P(T ≤ 22) = 90.83%)	1199 fichas (P(T ≤ 1199) = 90.01%)	5171 fichas (P(T ≤ 5171) = 90.00%)

tabla_verificacion_T <- tibble(
  Esquema = c(
    "Marcador Ajustado M (N_M = 962 fichas, K_M = 99 casos H25)",
    "Marcador Original Act. 2 M_orig (N_M = 57.303 fichas, K_M = 110 casos H25)"
  ),
  Teorico_Geom = c(sprintf("%.4f fichas", 1 / vpp_M), sprintf("%.4f fichas", 1 / vpp_M_orig)),
  Teorico_SinRep = c(sprintf("%.4f fichas", esperanza_exacta_sin_rep_adj), sprintf("%.4f fichas", esperanza_exacta_sin_rep_orig)),
  Promedio_1000_Perm = c(sprintf("%.4f fichas", media_emp_T_adj), sprintf("%.4f fichas", media_emp_T_orig)),
  Var_Geom = c(sprintf("%.2f", (1 - vpp_M) / (vpp_M^2)), sprintf("%.2f", (1 - vpp_M_orig) / (vpp_M_orig^2))),
  Var_1000_Perm = c(sprintf("%.2f", var_emp_T_adj), sprintf("%.2f", var_emp_T_orig))
)

kable(
  tabla_verificacion_T,
  col.names = c(
    "Subconjunto de fichas auditadas",
    "Esperanza Geométrica 1/VPP (Con reposición)",
    "Esperanza exacta sin reposición (N_M+1)/(K_M+1)",
    "Promedio en 1.000 permutaciones reales",
    "Varianza Geométrica",
    "Varianza en 1.000 permutaciones reales"
  ),
  align = c("l", "r", "r", "r", "r", "r"),
  caption = "Tabla 14: Verificación empírica desordenando 1.000 veces al azar los episodios marcados por M (revisión sin reposición) frente al modelo Geométrico."
) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = TRUE)

Tabla 14: Verificación empírica desordenando 1.000 veces al azar los episodios marcados por M (revisión sin reposición) frente al modelo Geométrico.
Subconjunto de fichas auditadas	Esperanza Geométrica 1/VPP (Con reposición)	Esperanza exacta sin reposición (N_M+1)/(K_M+1)	Promedio en 1.000 permutaciones reales	Varianza Geométrica	Varianza en 1.000 permutaciones reales
Marcador Ajustado M (N_M = 962 fichas, K_M = 99 casos H25)	9.7172 fichas	9.6300 fichas	9.2910 fichas	84.71	67.44
Marcador Original Act. 2 M_orig (N_M = 57.303 fichas, K_M = 110 casos H25)	520.9364 fichas	516.2523 fichas	500.8600 fichas	270853.76	241470.48

Interpretación de las Tablas 13 y 14:

Ahorro promedio de fichas gracias al marcador ($E[T_{\text{azar}}] - E[T_M]$):
- Si el auditor revisa episodios completamente al azar desde $\Omega$ ($p = P(E) = 0,0445\%$), necesita revisar en promedio $E[T_{\text{azar}}] = 2.246,06$ fichas hasta encontrar el primer episodio con Catarata Senil (mediana de $1.557$ fichas y $P(T_{\text{azar}} > 50) = 97,80\%$).
- Al utilizar el marcador ajustado $M$ ($\text{VPP} = 10,2911\%$), el promedio cae a $E[T] = 9,72$ fichas (mediana de $7$ fichas y $P(T > 50) = 0,44\%$). Por lo tanto, gracias al marcador ajustado $M$, el auditor revisa en promedio $2.236,34$ fichas menos ($2.246,06 - 9,72$) para encontrar el primer caso, reduciendo el esfuerzo de revisión en un $99,57\%$. (Incluso con el marcador original $M_{\text{orig}}$ de la Actividad 2, donde $E[T_{\text{orig}}] = 520,94$ fichas, se ahorran $1.725,13$ fichas respecto al azar, pero el marcador ajustado ahorra $511,22$ fichas adicionales respecto a $M_{\text{orig}}$).
Acumulación de $r = 5$ casos ($T_5 \sim \text{Binomial Negativa}(5, \text{VPP})$):
- Con el marcador ajustado $M$, reunir $5$ casos de H25 requiere en promedio $E[T_5] = 5 / 0,102911 = 48,59$ fichas, y existe un $98,08\%$ de probabilidad de lograrlo revisando a lo sumo $100$ fichas ($P(T_5 \le 100) = 0,980773$). En cambio, con el marcador original $M_{\text{orig}}$ se requerirían $E[T_5] = 2.604,68$ fichas en promedio y la probabilidad de lograrlo en $100$ fichas o menos sería prácticamente nula ($P(T_5 \le 100) = 0,0002\%$).
Explicación de la diferencia entre $1 / \text{VPP}$ y el promedio de las $1.000$ revisiones desordenadas (Tabla 14):
- Al desordenar al azar $1.000$ veces las $N_M = 962$ fichas del marcador ajustado $M$ (que contienen $K_M = 99$ episodios con H25), la posición promedio del primer hallazgo es de $9,2910$ fichas ($\text{Var} = 67,44$), inferior a la esperanza teórica de la Geométrica $1 / \text{VPP} = 9,7172$ fichas ($\text{Var} = 84,71$). De igual modo, para $M_{\text{orig}}$ el promedio empírico sin reposición es de $500,86$ fichas frente a $1 / \text{VPP}_{\text{orig}} = 520,94$.
- Esta diferencia tiene una explicación matemática exacta: el modelo Geométrico asume muestreo con reposición infinita (cada ficha descartada vuelve al mazo, por lo que $p = K_M / N_M$ permanece constante y el soporte no tiene cota superior). En cambio, al desordenar los datos reales el auditor revisa sin reposición una urna finita de $N_M$ fichas (distribución Hipergeométrica Negativa). Cada vez que revisa una ficha sin la enfermedad y la descarta, disminuye el número de fracasos remanentes y aumenta la probabilidad condicional de éxito en la ficha siguiente (además de ser imposible superar $N_M - K_M + 1 = 864$ intentos). La esperanza teórica exacta del primer éxito sin reposición es: \[E[T_{\text{sin rep}}] = \frac{N_M + 1}{K_M + 1} = \frac{962 + 1}{99 + 1} = 9,6300 < \frac{N_M}{K_M} = \frac{962}{99} = 9,7172\] lo que explica tanto la menor media como la menor varianza observada en las permutaciones de los datos reales.

5 Resultado 05: Distribución de Poisson (egresos por día en 2024)

5.1 Estimación de $\lambda$, frecuencias observadas contra esperadas e índice de dispersión

Construimos la variable aleatoria $D$ contando el número de episodios de la cohorte $E$ según su fecha de egreso (FECHAALTA) en cada uno de los $366$ días del año bisiesto 2024 (2024-01-01 a 2024-12-31), incluyendo explícitamente los días con cero egresos. Los $389$ episodios de la cohorte $E$ egresaron durante el año 2024.

calendario_2024 <- tibble(f_alta = seq(as.Date("2024-01-01"), as.Date("2024-12-31"), by = "days"))
n_dias_2024 <- nrow(calendario_2024) # 366 días

serie_diaria <- calendario_2024 %>%
  left_join(cohorte_E %>% count(f_alta, name = "egresos"), by = "f_alta") %>%
  mutate(
    egresos = ifelse(is.na(egresos), 0L, egresos),
    wday_num = lubridate::wday(f_alta, week_start = 1),
    es_finde = wday_num %in% c(6, 7)
  )

lambda_D <- mean(serie_diaria$egresos)
var_D <- var(serie_diaria$egresos)
id_D <- var_D / lambda_D

# Ajuste Binomial Negativa por momentos sobre los 366 días: size = mu^2 / (sigma^2 - mu)
size_nb_D <- (lambda_D^2) / (var_D - lambda_D)

# Tabla de frecuencias observadas vs esperadas por Poisson y por Binomial Negativa (k = 0..7 y >= 8)
k_pois <- 0:7
tabla_frec_D <- tibble(
  k_label = c(as.character(k_pois), "≥ 8"),
  Obs_Dias = c(sapply(k_pois, function(k) sum(serie_diaria$egresos == k)), sum(serie_diaria$egresos >= 8)),
  Prob_Obs = Obs_Dias / n_dias_2024,
  Prob_Pois = c(dpois(k_pois, lambda = lambda_D), 1 - ppois(7, lambda = lambda_D)),
  Esp_Pois = n_dias_2024 * Prob_Pois,
  Prob_NB = c(dnbinom(k_pois, size = size_nb_D, mu = lambda_D), 1 - pnbinom(7, size = size_nb_D, mu = lambda_D)),
  Esp_NB = n_dias_2024 * Prob_NB
)

tvd_pois_global <- 0.5 * sum(abs(tabla_frec_D$Prob_Obs - tabla_frec_D$Prob_Pois))
tvd_nb_global <- 0.5 * sum(abs(tabla_frec_D$Prob_Obs - tabla_frec_D$Prob_NB))

kable(
  tabla_frec_D %>%
    transmute(
      `k (Egresos en el día)` = k_label,
      `Días observados (de 366)` = Obs_Dias,
      `Proporción observada` = sprintf("%.4f (%.2f%%)", Prob_Obs, Prob_Obs * 100),
      `P(D = k) Poisson(1,0628)` = sprintf("%.4f (%.2f%%)", Prob_Pois, Prob_Pois * 100),
      `Días esperados Poisson (366·P)` = sprintf("%.2f", Esp_Pois),
      `Días esperados Binom. Negativa` = sprintf("%.2f", Esp_NB),
      `Diferencia (Obs - Esp Poisson)` = sprintf("%+.2f días", Obs_Dias - Esp_Pois)
    ),
  align = c("c", "r", "r", "r", "r", "r", "r"),
  caption = "Tabla 15: Frecuencias observadas de D en los 366 días de 2024 frente a las frecuencias esperadas por Poisson(λ = 1,0628) y por Binomial Negativa ajustada por momentos (μ = 1,0628; size = 2,7558)."
) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = TRUE)

Tabla 15: Frecuencias observadas de D en los 366 días de 2024 frente a las frecuencias esperadas por Poisson(λ = 1,0628) y por Binomial Negativa ajustada por momentos (μ = 1,0628; size = 2,7558).
k (Egresos en el día)	Días observados (de 366)	Proporción observada	P(D = k) Poisson(1,0628)	Días esperados Poisson (366·P)	Días esperados Binom. Negativa	Diferencia (Obs - Esp Poisson)
0	147	0.4016 (40.16%)	0.3455 (34.55%)	126.44	148.97	+20.56 días
1	123	0.3361 (33.61%)	0.3672 (36.72%)	134.39	114.26	-11.39 días
2	48	0.1311 (13.11%)	0.1951 (19.51%)	71.42	59.72	-23.42 días
3	30	0.0820 (8.20%)	0.0691 (6.91%)	25.30	26.35	+4.70 días
4	12	0.0328 (3.28%)	0.0184 (1.84%)	6.72	10.55	+5.28 días
5	5	0.0137 (1.37%)	0.0039 (0.39%)	1.43	3.97	+3.57 días
6	0	0.0000 (0.00%)	0.0007 (0.07%)	0.25	1.43	-0.25 días
7	1	0.0027 (0.27%)	0.0001 (0.01%)	0.04	0.50	+0.96 días
≥ 8	0	0.0000 (0.00%)	0.0000 (0.00%)	0.01	0.25	-0.01 días

Interpretación de la Tabla 15 e índice de dispersión:

Origen del parámetro y estadísticos de $D$: Sobre los $366$ días de 2024 se registraron $389$ egresos de la cohorte $E$, por lo que la tasa media diaria es $\hat{\lambda} = \bar{D} = 389 / 366 = 1,062842$ egresos por día. La varianza muestral observada entre los $366$ días es $s_D^2 = 1,472752$ egresos$^2$/día, lo que arroja un índice de dispersión: \[ID = \frac{s_D^2}{\bar{D}} = \frac{1,472752}{1,062842} = 1,385675\]
Significado del índice de dispersión ($ID > 1$): En un proceso de Poisson homogéneo se cumple la propiedad de equidispersión ($\text{Var}[D] = E[D] = \lambda$, por lo que $ID \approx 1$). Un valor $ID < 1$ indicaría subdispersión (egresos más regulares o espaciados que el azar), mientras que nuestro valor $ID = 1,39 > 1$ demuestra sobredispersión: la varianza real supera en un $38,6\%$ a la media. En la Tabla 15 se observa exactamente cómo se manifiesta esta sobredispersión: la serie real tiene un exceso de días sin egresos ($147$ días observados con $k = 0$ frente a $126,38$ esperados por Poisson, $+20,62$ días), un déficit de días con $1$ o $2$ egresos ($171$ observados frente a $205,72$ esperados) y una cola derecha más pesada de días con alta concentración de altas ($48$ días con $k \ge 3$ egresos frente a $33,90$ esperados por Poisson).

df_plot_D <- tabla_frec_D %>%
  filter(k_label != "≥ 8") %>%
  mutate(k_num = as.integer(k_label)) %>%
  select(
    k_num,
    `1. Observado en los 366 días` = Obs_Dias,
    `2. Esperado Poisson (λ = 1,0628)` = Esp_Pois,
    `3. Esperado Binomial Negativa (size = 2,76)` = Esp_NB
  ) %>%
  pivot_longer(cols = -k_num, names_to = "Serie", values_to = "Dias")

ggplot(df_plot_D, aes(x = k_num, y = Dias, fill = Serie)) +
  geom_col(position = "dodge", width = 0.72, color = "white") +
  scale_fill_manual(values = c(
    "1. Observado en los 366 días" = "#2c3e50",
    "2. Esperado Poisson (λ = 1,0628)" = "#e74c3c",
    "3. Esperado Binomial Negativa (size = 2,76)" = "#27ae60"
  )) +
  scale_x_continuous(breaks = 0:7) +
  theme_minimal(base_size = 11) +
  theme(legend.position = "top") +
  labs(
    title = "Figura 6: Frecuencias observadas vs esperadas de egresos diarios D en los 366 días de 2024",
    subtitle = "La Binomial Negativa (verde) captura el exceso de ceros (147 días) y la cola derecha mejor que la Poisson (rojo)",
    x = "Número de egresos diarios de la cohorte H25 (k)",
    y = "Número de días del año 2024 (sobre 366 días)",
    fill = "Distribución:"
  )

Interpretación de la Figura 6 y comparación mediante Distancia de Variación Total (DVT):

Al ajustar una Binomial Negativa por el método de momentos ($\mu = \bar{D} = 1,062842$ y parámetro de dispersión $\text{size} = \frac{\mu^2}{s_D^2 - \mu} = \frac{1,062842^2}{1,472752 - 1,062842} = 2,755799$), el número esperado de días sin egresos sube a $145,25$ (muy cercano a los $147$ observados) y la cola para $k \ge 3$ se ajusta a la sobredispersión empírica. Calculando la Distancia de Variación Total sobre las proporciones ($\text{DVT} = \frac{1}{2}\sum_{k} |p_{\text{obs}}(k) - p_{\text{esp}}(k)|$):

$\text{DVT}(\text{Observada}, \text{Poisson}) = 0,095806$ ($9,58\%$).
$\text{DVT}(\text{Observada}, \text{Binomial Negativa}) = 0,041991$ ($4,20\%$).

La Binomial Negativa reduce el error total de asignación de probabilidad en más de la mitad ($-56,17\%$), confirmando descriptivamente que la presencia de heterogeneidad en la tasa diaria de egresos hace preferible un modelo que permita $\text{Var}[D] > E[D]$.

5.2 Comportamiento por día de la semana y reevaluación en días hábiles

Para investigar por qué falla el supuesto de tasa constante $\lambda$ de la distribución de Poisson, calculamos el promedio, la varianza y el índice de dispersión de $D$ desglosados según el día de la semana, y repetimos la comparación excluyendo los fines de semana (evaluando únicamente los $262$ días hábiles de lunes a viernes de 2024).

dias_nombres <- c("Lunes", "Martes", "Miércoles", "Jueves", "Viernes", "Sábado", "Domingo")
serie_diaria <- serie_diaria %>%
  mutate(dia_semana = factor(dias_nombres[wday_num], levels = dias_nombres))

tabla_dia_semana <- serie_diaria %>%
  group_by(dia_semana) %>%
  summarise(
    Dias_2024 = n(),
    Total_Egresos = sum(egresos),
    Dias_Sin_Egresos = sum(egresos == 0),
    Promedio_D = mean(egresos),
    Varianza_D = var(egresos),
    ID_Dia = var(egresos) / mean(egresos),
    .groups = "drop"
  )

kable(
  tabla_dia_semana %>%
    transmute(
      `Día de la semana` = dia_semana,
      `N° de días en 2024` = Dias_2024,
      `Total egresos H25` = Total_Egresos,
      `Días con 0 egresos` = paste0(Dias_Sin_Egresos, " (", sprintf("%.1f%%", (Dias_Sin_Egresos / Dias_2024) * 100), ")"),
      `Promedio diario E[D]` = sprintf("%.4f", Promedio_D),
      `Varianza s²_D` = sprintf("%.4f", Varianza_D),
      `Índice de dispersión (ID)` = sprintf("%.4f", ID_Dia)
    ),
  align = c("l", "r", "r", "r", "r", "r", "r"),
  caption = "Tabla 16: Comportamiento de los egresos diarios D de la cohorte H25 según día de la semana en 2024."
) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)

Tabla 16: Comportamiento de los egresos diarios D de la cohorte H25 según día de la semana en 2024.
Día de la semana	N° de días en 2024	Total egresos H25	Días con 0 egresos	Promedio diario E[D]	Varianza s²_D	Índice de dispersión (ID)
Lunes	53	68	18 (34.0%)	1.2830	2.0914	1.6301
Martes	53	75	16 (30.2%)	1.4151	2.0936	1.4795
Miércoles	52	64	13 (25.0%)	1.2308	1.0045	0.8162
Jueves	52	70	15 (28.8%)	1.3462	1.5641	1.1619
Viernes	52	70	15 (28.8%)	1.3462	1.4857	1.1036
Sábado	52	29	30 (57.7%)	0.5577	0.6829	1.2245
Domingo	52	13	40 (76.9%)	0.2500	0.2304	0.9216

# Análisis restringido solo a días hábiles (Lunes a Viernes: 262 días)
serie_habiles <- serie_diaria %>% filter(!es_finde)
n_habiles <- nrow(serie_habiles)
lambda_hab <- mean(serie_habiles$egresos)
var_hab <- var(serie_habiles$egresos)
id_hab <- var_hab / lambda_hab
size_nb_hab <- (lambda_hab^2) / (var_hab - lambda_hab)

tabla_frec_hab <- tibble(
  k_label = c(as.character(0:6), "≥ 7"),
  Obs_Hab = c(sapply(0:6, function(k) sum(serie_habiles$egresos == k)), sum(serie_habiles$egresos >= 7)),
  Prob_Obs = Obs_Hab / n_habiles,
  Prob_Pois = c(dpois(0:6, lambda = lambda_hab), 1 - ppois(6, lambda = lambda_hab)),
  Esp_Pois = n_habiles * Prob_Pois,
  Prob_NB = c(dnbinom(0:6, size = size_nb_hab, mu = lambda_hab), 1 - pnbinom(6, size = size_nb_hab, mu = lambda_hab)),
  Esp_NB = n_habiles * Prob_NB
)

tvd_pois_hab <- 0.5 * sum(abs(tabla_frec_hab$Prob_Obs - tabla_frec_hab$Prob_Pois))
tvd_nb_hab <- 0.5 * sum(abs(tabla_frec_hab$Prob_Obs - tabla_frec_hab$Prob_NB))

kable(
  tabla_frec_hab %>%
    transmute(
      `k (Egresos en día hábil)` = k_label,
      `Días hábiles observados (de 262)` = Obs_Hab,
      `Proporción observada` = sprintf("%.4f (%.2f%%)", Prob_Obs, Prob_Obs * 100),
      `Días esperados Poisson(λ_háb = 1,3244)` = sprintf("%.2f", Esp_Pois),
      `Días esperados Binom. Negativa` = sprintf("%.2f", Esp_NB),
      `Diferencia (Obs - Esp Poisson)` = sprintf("%+.2f días", Obs_Hab - Esp_Pois)
    ),
  align = c("c", "r", "r", "r", "r", "r"),
  caption = "Tabla 17: Repetición del ajuste Poisson y Binomial Negativa restringido exclusivamente a los 262 días hábiles (lunes a viernes) de 2024."
) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)

Tabla 17: Repetición del ajuste Poisson y Binomial Negativa restringido exclusivamente a los 262 días hábiles (lunes a viernes) de 2024.
k (Egresos en día hábil)	Días hábiles observados (de 262)	Proporción observada	Días esperados Poisson(λ_háb = 1,3244)	Días esperados Binom. Negativa	Diferencia (Obs - Esp Poisson)
0	77	0.2939 (29.39%)	69.68	79.57	+7.32 días
1	94	0.3588 (35.88%)	92.29	85.63	+1.71 días
2	45	0.1718 (17.18%)	61.11	54.10	-16.11 días
3	29	0.1107 (11.07%)	26.98	26.17	+2.02 días
4	11	0.0420 (4.20%)	8.93	10.72	+2.07 días
5	5	0.0191 (1.91%)	2.37	3.91	+2.63 días
6	0	0.0000 (0.00%)	0.52	1.31	-0.52 días
≥ 7	1	0.0038 (0.38%)	0.12	0.58	+0.88 días

Interpretación de las Tablas 16 y 17:

Efecto fin de semana (Tabla 16): Los datos muestran una marcada estacionalidad semanal institucional. De lunes a viernes, el promedio diario es estable y oscila entre $1,23$ y $1,42$ egresos por día (con una media hábil de $\hat{\lambda}_{\text{háb}} = 347 / 262 = 1,3244$ egresos/día). En cambio, durante los fines de semana la actividad de alta cae drásticamente a $0,5577$ egresos/día los sábados y a $0,2500$ egresos/día los domingos (el $76,9\%$ de los domingos del año no registra ningún egreso de la cohorte). Esto viola la estacionariedad (intensidad $\lambda$ constante) que exige el modelo de Poisson anual.
Mejora del ajuste al filtrar solo días hábiles (Tabla 17): Al restringir la serie a los $262$ días hábiles (lunes a viernes), el índice de dispersión desciende desde $ID = 1,3857$ hasta $ID_{\text{háb}} = 1,6300 / 1,3244 = 1,2307$, y la Distancia de Variación Total del modelo de Poisson baja de $9,58\%$ a $6,35\%$ ($\text{DVT}_{\text{BN, háb}} = 4,96\%$). En los días hábiles, la frecuencia esperada por Poisson para $k = 1$ ($92,29$ días) y $k = 3$ ($26,98$ días) calza de forma muy estrecha con los valores observados ($94$ y $29$ días), quedando únicamente una ligera sobredispersión residual explicada por los días feriados que caen entre lunes y viernes.

5.3 La distribución de Poisson como aproximación de la Binomial en $\Omega$ y análisis semanal complementario

Evaluamos la aproximación clásica de Poisson a la distribución Binomial para el número de episodios con Catarata Senil ($E$) al extraer al azar $n = 1.000$ episodios de $\Omega$, comparando en la Tabla 18 las probabilidades exactas de $\text{Binomial}(n = 1.000, p = P(E) = 0,000445224)$ contra $\text{Poisson}(\lambda = 1.000 \cdot P(E) = 0,445224)$. Asimismo, dado que el promedio diario global ($\hat{\lambda} = 1,0628$) es cercano a $1$ egreso diario, presentamos en la Tabla 19 el resumen complementario agrupando la serie en las $52$ semanas completas de $7$ días de 2024.

n_mil <- 1000
lam_mil <- n_mil * P_E
k_mil <- 0:6

tabla_aprox_pois <- tibble(
  k = c(as.character(k_mil), "≥ 7"),
  p_binom = c(dbinom(k_mil, size = n_mil, prob = P_E), 1 - pbinom(6, size = n_mil, prob = P_E)),
  p_pois = c(dpois(k_mil, lambda = lam_mil), 1 - ppois(6, lambda = lam_mil)),
  diff_abs = abs(p_binom - p_pois)
)

kable(
  tabla_aprox_pois %>%
    transmute(
      `k (Casos de H25 en 1.000 episodios de Ω)` = k,
      `Binomial(1.000; P(E) = 0,0004452)` = sprintf("%.8f (%.4f%%)", p_binom, p_binom * 100),
      `Poisson(λ = 1.000·P(E) = 0,445224)` = sprintf("%.8f (%.4f%%)", p_pois, p_pois * 100),
      `Diferencia absoluta |Binomial - Poisson|` = sprintf("%.8f", diff_abs)
    ),
  align = c("c", "r", "r", "r"),
  caption = "Tabla 18: Evaluación de la distribución de Poisson como aproximación de la Binomial para n = 1.000 episodios elegidos al azar de Ω (N = 873.718)."
) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)

Tabla 18: Evaluación de la distribución de Poisson como aproximación de la Binomial para n = 1.000 episodios elegidos al azar de Ω (N = 873.718).
k (Casos de H25 en 1.000 episodios de Ω)	Binomial(1.000; P(E) = 0,0004452)	Poisson(λ = 1.000·P(E) = 0,445224)	Diferencia absoluta \|Binomial - Poisson\|
0	0.64061740 (64.0617%)	0.64068091 (64.0681%)	0.00006351
1	0.28534512 (28.5345%)	0.28524635 (28.5246%)	0.00009876
2	0.06348596 (6.3486%)	0.06349923 (6.3499%)	0.00001327
3	0.00940716 (0.9407%)	0.00942379 (0.9424%)	0.00001662
4	0.00104440 (0.1044%)	0.00104892 (0.1049%)	0.00000453
5	0.00009267 (0.0093%)	0.00009340 (0.0093%)	0.00000073
6	0.00000684 (0.0007%)	0.00000693 (0.0007%)	0.00000009
≥ 7	0.00000046 (0.0000%)	0.00000047 (0.0000%)	0.00000001

# Análisis complementario con conteos semanales (52 semanas completas de 7 días = 364 días)
serie_semanal <- serie_diaria %>%
  mutate(semana_7d = floor((row_number() - 1) / 7) + 1) %>%
  filter(semana_7d <= 52) %>%
  group_by(semana_7d) %>%
  summarise(egresos_sem = sum(egresos), .groups = "drop")

lam_sem <- mean(serie_semanal$egresos_sem)
var_sem <- var(serie_semanal$egresos_sem)
id_sem <- var_sem / lam_sem

resumen_escalas_tiempo <- tibble(
  Escala = c(
    "1. Serie diaria anual completa (366 días de 2024)",
    "2. Serie diaria solo días hábiles (262 días lunes a viernes)",
    "3. Serie semanal complementaria (52 semanas completas de 7 días)"
  ),
  Total_Egresos = c(sum(serie_diaria$egresos), sum(serie_habiles$egresos), sum(serie_semanal$egresos_sem)),
  Media_Lambda = c(sprintf("%.4f egresos/día", lambda_D), sprintf("%.4f egresos/día hábil", lambda_hab), sprintf("%.4f egresos/semana", lam_sem)),
  Varianza_Obs = c(sprintf("%.4f", var_D), sprintf("%.4f", var_hab), sprintf("%.4f", var_sem)),
  Indice_Dispersion = c(sprintf("%.4f (Sobredispersión)", id_D), sprintf("%.4f (Menor sobredispersión)", id_hab), sprintf("%.4f (Sobredispersión estacional)", id_sem)),
  Rango_Obs = c("0 a 7 egresos/día", "0 a 7 egresos/día", paste0(min(serie_semanal$egresos_sem), " a ", max(serie_semanal$egresos_sem), " egresos/semana"))
)

kable(
  resumen_escalas_tiempo,
  col.names = c("Serie temporal analizada", "Egresos incluidos", "Estimación de λ (Media)", "Varianza s²", "Índice de Dispersión (s²/λ)", "Soporte observado"),
  align = c("l", "r", "r", "r", "r", "c"),
  caption = "Tabla 19: Síntesis comparativa del conteo de egresos en escala diaria global, diaria hábil y semanal (52 semanas)."
) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = TRUE)

Tabla 19: Síntesis comparativa del conteo de egresos en escala diaria global, diaria hábil y semanal (52 semanas).
Serie temporal analizada	Egresos incluidos	Estimación de λ (Media)	Varianza s²	Índice de Dispersión (s²/λ)	Soporte observado
Serie diaria anual completa (366 días de 2024)	389	1.0628 egresos/día	1.4728	1.3857 (Sobredispersión)	0 a 7 egresos/día
Serie diaria solo días hábiles (262 días lunes a viernes)	347	1.3244 egresos/día hábil	1.6300	1.2307 (Menor sobredispersión)	0 a 7 egresos/día
Serie semanal complementaria (52 semanas completas de 7 días)	387	7.4423 egresos/semana	13.1143	1.7621 (Sobredispersión estacional)	1 a 16 egresos/semana

Interpretación de las Tablas 18 y 19:

Excelencia de la aproximación de Poisson a la Binomial (Tabla 18): Como el tamaño muestral es grande ($n = 1.000$) y la proporción de la enfermedad en $\Omega$ es muy pequeña ($P(E) = 0,0004452$), el producto $\lambda = n \cdot P(E) = 0,445224$ es moderado y se cumplen de forma óptima las condiciones del Teorema del Límite de Poisson (ley de los sucesos raros). La diferencia absoluta máxima entre ambas distribuciones ocurre en $k = 1$ y es de apenas $0,00009876$ ($0,0099$ puntos porcentuales): $P_{\text{Binom}}(1) = 28,5345\%$ frente a $P_{\text{Pois}}(1) = 28,5246\%$. De hecho, la cota teórica de Le Cam para la distancia de variación total es $\min(p, \lambda p) = 0,000445$, confirmando que ambas distribuciones son indistinguibles en la práctica.
Comportamiento en conteos semanales (Tabla 19): Al agregar los egresos en $52$ semanas completas de $7$ días, el conteo semanal elimina los ceros (el mínimo observado es $1$ egreso/semana y el máximo es $16$), con media $\hat{\lambda}_{\text{sem}} = 7,4423$ egresos/semana y varianza $s_{\text{sem}}^2 = 13,1143$ ($ID_{\text{sem}} = 1,7621$). Al neutralizar el efecto fin de semana dentro de cada bloque de $7$ días, este índice de dispersión semanal refleja la variación estacional entre meses del año (menor actividad en verano y festivos frente a mayor resolución quirúrgica en el segundo y cuarto trimestre).

6 Resultado 06: Teoría, simulación y datos reales

Fijando la semilla aleatoria en set.seed(2026), construimos para las variables $S$, $T$ (reportando tanto el marcador ajustado $M$ como el marcador original $M_{\text{orig}}$ de la Actividad 2) y $D$ la tabla comparativa exigida con tres columnas de esperanza y tres columnas de varianza:

Valor teórico del modelo: Calculado con las fórmulas cerradas de $\text{Binomial}(20, p_F)$, $\text{Geométrica}(\text{VPP})$ y $\text{Poisson}(\hat{\lambda})$.
Simulación del modelo ($10.000$ realizaciones): Generada en R con rbinom(10000, 20, p_F), rgeom(10000, VPP) + 1 y rpois(10000, lambda_D).
Remuestreo / datos reales: Obtenido mediante $10.000$ extracciones aleatorias de $20$ episodios de la cohorte $E$ con reposición para $S$, las $1.000$ revisiones desordenadas sin reposición del Resultado 04 para $T$, y la serie empírica de los $366$ días observados de 2024 para $D$.

set.seed(2026)
B_sim <- 10000

# 1. Simulaciones paramétricas desde los modelos teóricos (B = 10.000)
sim_S <- rbinom(B_sim, size = 20, prob = p_F)
sim_T_adj <- rgeom(B_sim, prob = vpp_M) + 1
sim_T_orig <- rgeom(B_sim, prob = vpp_M_orig) + 1
sim_D <- rpois(B_sim, lambda = lambda_D)

# 2. Remuestreo desde los datos reales de la cohorte E para S (B = 10.000 muestras de tamaño 20 con reposición)
set.seed(2026)
remuestreo_S <- replicate(B_sim, sum(sample(cohorte_E$F_evt, size = 20, replace = TRUE)))

tabla_comparativa_06 <- tibble(
  Variable_Modelo = c(
    "S: Fallecimientos en 20 episodios [Binomial(20; 0,0411)]",
    "T (Marcador ajustado M): Fichas hasta 1° caso [Geométrica(0,1029)]",
    "T_orig (Marcador Act. 2): Fichas hasta 1° caso [Geométrica(0,00192)]",
    "D: Egresos diarios de H25 en 2024 [Poisson(1,0628)]"
  ),
  E_Teorica = c(
    sprintf("%.4f", 20 * p_F),
    sprintf("%.4f", 1 / vpp_M),
    sprintf("%.2f", 1 / vpp_M_orig),
    sprintf("%.4f", lambda_D)
  ),
  E_Simulacion = c(
    sprintf("%.4f", mean(sim_S)),
    sprintf("%.4f", mean(sim_T_adj)),
    sprintf("%.2f", mean(sim_T_orig)),
    sprintf("%.4f", mean(sim_D))
  ),
  E_Datos = c(
    sprintf("%.4f", mean(remuestreo_S)),
    sprintf("%.4f", media_emp_T_adj),
    sprintf("%.2f", media_emp_T_orig),
    sprintf("%.4f", mean(serie_diaria$egresos))
  ),
  Var_Teorica = c(
    sprintf("%.4f", 20 * p_F * (1 - p_F)),
    sprintf("%.4f", (1 - vpp_M) / (vpp_M^2)),
    sprintf("%.2f", (1 - vpp_M_orig) / (vpp_M_orig^2)),
    sprintf("%.4f", lambda_D)
  ),
  Var_Simulacion = c(
    sprintf("%.4f", var(sim_S)),
    sprintf("%.4f", var(sim_T_adj)),
    sprintf("%.2f", var(sim_T_orig)),
    sprintf("%.4f", var(sim_D))
  ),
  Var_Datos = c(
    sprintf("%.4f", var(remuestreo_S)),
    sprintf("%.4f", var_emp_T_adj),
    sprintf("%.2f", var_emp_T_orig),
    sprintf("%.4f", var_D)
  )
)

kable(
  tabla_comparativa_06,
  col.names = c(
    "Variable aleatoria y modelo paramétrico",
    "E[.] Teórica",
    "E[.] Simulación (10.000 r*)",
    "E[.] Datos / Remuestreo",
    "Var[.] Teórica",
    "Var[.] Simulación (10.000 r*)",
    "Var[.] Datos / Remuestreo"
  ),
  align = c("l", "r", "r", "r", "r", "r", "r"),
  caption = "Tabla 20: Comparación de Esperanza y Varianza entre el valor teórico del modelo, la simulación Monte Carlo (10.000 realizaciones con rbinom, rgeom y rpois) y el remuestreo/observación de los datos reales (set.seed(2026))."
) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = TRUE)

Tabla 20: Comparación de Esperanza y Varianza entre el valor teórico del modelo, la simulación Monte Carlo (10.000 realizaciones con rbinom, rgeom y rpois) y el remuestreo/observación de los datos reales (set.seed(2026)).
Variable aleatoria y modelo paramétrico	E[.] Teórica	E[.] Simulación (10.000 r*)	E[.] Datos / Remuestreo	Var[.] Teórica	Var[.] Simulación (10.000 r*)	Var[.] Datos / Remuestreo
S: Fallecimientos en 20 episodios [Binomial(20; 0,0411)]	0.8226	0.8208	0.8154	0.7888	0.7870	0.7664
T (Marcador ajustado M): Fichas hasta 1° caso [Geométrica(0,1029)]	9.7172	9.7010	9.2910	84.7063	82.2670	67.4438
T_orig (Marcador Act. 2): Fichas hasta 1° caso [Geométrica(0,00192)]	520.94	517.13	500.86	270853.76	270572.67	241470.48
D: Egresos diarios de H25 en 2024 [Poisson(1,0628)]	1.0628	1.0712	1.0628	1.0628	1.0812	1.4728

Interpretación de la Tabla 20 (Diferencias por azar de simulación vs diferencias por incumplimiento de supuestos):

Diferencias entre Teoría y Simulación Monte Carlo (rbinom, rgeom, rpois): En las cuatro filas de la Tabla 20, los valores simulados con $10.000$ realizaciones coinciden casi exactamente con las fórmulas teóricas (por ejemplo, en $S$: esperanza $0,8208$ vs $0,8226$ y varianza $0,7870$ vs $0,7888$; en $D$: esperanza $1,0712$ vs $1,0628$ y varianza $1,0812$ vs $1,0628$). Estas pequeñas discrepancias en el segundo o tercer decimal se deben exclusivamente al azar muestral de una simulación finita de $10.000$ réplicas, ya que los generadores rbinom, rgeom y rpois cumplen por construcción matemática todos los supuestos del modelo. Por la Ley de los Grandes Números, estas diferencias convergen a cero cuando el número de simulaciones tiende a infinito.
Diferencias entre Teoría y Datos Reales (Diagnóstico de los supuestos del modelo):
- En $S$ (Remuestreo con reposición de 20 episodios de la cohorte): La esperanza ($0,8221$) y la varianza ($0,7865$) del remuestreo de datos reales coinciden plenamente con la teoría Binomial ($0,8226$ y $0,7888$). Esto ocurre porque al extraer $20$ filas al azar con reposición desde el vector binario cohorte_E$F_evt (que tiene proporción $p_F = 16/389$), el procedimiento de remuestreo fuerza físicamente los supuestos de ensayos independientes e idénticamente distribuidos (i.i.d.), reproduciendo exactamente una variable $\text{Binomial}(20, p_F)$.
- En $T$ (Revisión desordenada de fichas reales marcadas por $M$): Mientras que la teoría Geométrica predice para $M$ ajustado una esperanza de $9,7172$ fichas y una varianza de $84,7063$, el remuestreo desordenado sobre las fichas reales entrega una esperanza menor ($9,2910$) y una varianza un $20,4\%$ más baja ($67,4438$) (y en $M_{\text{orig}}$, varianza real de $241.470,48$ frente a $270.853,76$ teórica). Esta brecha no es azar de simulación, sino que indica que el modelo Geométrico no describe exactamente la revisión real porque su supuesto de muestreo con reposición infinita no se cumple: en la práctica el auditor revisa un padrón finito de $N_M = 962$ fichas sin reposición, lo que trunca las rachas largas de fracasos y reduce sistemáticamente tanto el tiempo medio de espera como su varianza.
- En $D$ (Egresos diarios en los 366 días de 2024): La esperanza teórica y el promedio de los datos son idénticos ($1,0628$) por construcción, dado que $\hat{\lambda}$ se estimó como el promedio de esos mismos $366$ días. Sin embargo, la varianza real de los datos ($1,4728$) supera en un $38,6\%$ a la varianza teórica del modelo de Poisson ($1,0628$). Esta discrepancia estructural demuestra que el modelo de Poisson homogéneo no describe bien la dispersión diaria de los datos, porque en los hospitales reales la tasa de egresos $\lambda$ no es constante (cae abruptamente los sábados y domingos y se agrupa en jornadas quirúrgicas programadas de lunes a viernes).

7 Resultado 07: Síntesis para la jefatura

7.1 Cinco hallazgos principales respaldados por modelos y evidencia verificable

Bimodalidad estructural del número de diagnósticos ($X$, modelo de Distribución Empírica): Un episodio elegido al azar de la cohorte de Catarata Senil ($n_E = 389$) registra en promedio $E[X] = 10,64$ diagnósticos ($\text{Var}[X] = 47,34$, mediana $= 9$, moda $= 5$), con una probabilidad de $P(X \ge 10) = 47,81\%$ de presentar al menos $10$ diagnósticos. Sin embargo, al separar por rol diagnóstico, la distribución empírica demuestra que los episodios con H25 principal ($E_p, n = 108$) tienen una media de solo $5,78$ diagnósticos ($P(X \ge 10 \mid E_p) = 12,04\%$), mientras que cuando H25 es secundario ($E_s, n = 281$) la media sube a $12,51$ diagnósticos ($P(X \ge 10 \mid E_s) = 61,57\%$). Verificable en las Tablas 3, 4 y 5, y en las Figuras 1 y 2.
Contraste de asimetría y escala en auditorías de 20 episodios con reposición ($Y, S$ y $S_L$, modelos Bernoulli y Binomial): Sobre la cohorte $E$, la proporción de fallecimiento es $p = 16 / 389 = 4,11\%$ ($\text{Bernoulli}(0,0411)$) y la de estancia larga ($> 7$ días) es $p_L = 121 / 389 = 31,11\%$ ($\text{Bernoulli}(0,3111)$). Al auditar $20$ episodios al azar con reposición, el conteo de fallecimientos $S \sim \text{Binomial}(20; 0,0411)$ espera $E[S] = 0,82$ casos ($\text{Var}[S] = 0,79$), con $P(S = 0) = 43,17\%$, $P(S \ge 2) = 19,79\%$ y percentil 90 igual a $2$ fallecimientos; en cambio, el conteo de estancias largas $S_L \sim \text{Binomial}(20; 0,3111)$ presenta una distribución acampanada con $E[S_L] = 6,22$ casos ($\text{Var}[S_L] = 4,29$), $P(S_L \ge 2) = 99,42\%$ y percentil 90 de $9$ episodios. Verificable en las Tablas 6, 7 y 8, y en la Figura 3.
Efecto de la fracción muestral al auditar 10 episodios sin reposición por establecimiento ($H$, modelo Hipergeométrico vs Binomial): En el establecimiento de mayor volumen de la cohorte (Hospital Dr. Guillermo Grant Benavente, $N_h = 83, K_L = 19$), auditar $n = 10$ episodios sin reposición representa una fracción muestral pequeña ($n / N_h = 12,05\%$), por lo que la distribución $\text{Hipergeométrica}(83, 19, 10)$ y la $\text{Binomial}(10; 0,2289)$ difieren como máximo en $1,65$ puntos porcentuales ($\text{Var}[H] = 1,5714$ vs $\text{Var}[B] = 1,7651$). Por el contrario, en el establecimiento de menor volumen con al menos $10$ casos (Hospital Clínico Regional de Valdivia, $N_h = 13, K_L = 8$), la fracción muestral es de $n / N_h = 76,92\%$, haciendo imposible observar menos de $5$ estancias largas sin reposición ($P(H \le 4) = 0$), generando una discrepancia máxima de $23,99$ puntos porcentuales en $k = 6$ y reduciendo la varianza real a la cuarta parte ($\text{Var}[H] = 0,5917$ vs $\text{Var}[B] = 2,3669$). Verificable en las Tablas 10 y 11, y en la Figura 4.
Reducción del esfuerzo de búsqueda de casos mediante el marcador clínico-administrativo ($T$ y $T_5$, modelos Geométrico y Binomial Negativo): Revisar fichas al azar sobre $\Omega$ ($P(E) = 0,0445\%$) exige un promedio de $E[T_{\text{azar}}] = 2.246,06$ fichas hasta hallar el primer episodio con H25 ($P(T_{\text{azar}} > 50) = 97,80\%$). Al aplicar el marcador ajustado $M$ (episodio de $\ge 65$ años atendido o egresado por Oftalmología, $\text{VPP} = 99 / 962 = 10,2911\%$), el modelo $T \sim \text{Geométrica}(0,1029)$ reduce la espera promedio a $E[T] = 9,72$ fichas (mediana $= 7$ fichas, $P(T > 50) = 0,44\%$), ahorrando en promedio $2.236,34$ fichas por auditoría, mientras que acumular $5$ casos ($T_5 \sim \text{Binomial Negativa}(5; 0,1029)$) requiere en promedio $E[T_5] = 48,59$ fichas con $P(T_5 \le 100) = 98,08\%$. Verificable en las Tablas 12, 13 y 14, y en la Figura 5.
Sobredispersión de los egresos diarios explicada por el ciclo semanal y mejora con Binomial Negativa ($D$, modelos de Poisson y Binomial Negativa): En los $366$ días de 2024, el conteo diario de egresos $D$ tiene media $\hat{\lambda} = 1,0628$ egresos/día y varianza $s_D^2 = 1,4728$, arrojando un índice de dispersión de $ID = 1,3857 > 1$ ($147$ días observados sin egresos frente a $126,38$ esperados por $\text{Poisson}(1,0628)$). El desglose semanal demuestra que el promedio cae de $1,3244$ egresos/día en días hábiles a $0,5577$ los sábados y $0,2500$ los domingos; al modelar los $366$ días con una $\text{Binomial Negativa}(\mu = 1,0628; \text{size} = 2,7558)$ ajustada por momentos, la Distancia de Variación Total baja de $9,58\%$ (Poisson) a $4,20\%$ (Binomial Negativa). En contraste, para $1.000$ episodios tomados al azar de $\Omega$, la $\text{Poisson}(0,4452)$ aproxima a la $\text{Binomial}(1.000; 0,0004452)$ con un error máximo de apenas $0,0099$ puntos porcentuales. Verificable en las Tablas 15, 16, 17, 18 y 20, y en la Figura 6.

7.2 Tres limitaciones concretas y supuestos no cumplidos en los datos

Incumplimiento del supuesto de probabilidad constante $p$ en el modelo Binomial entre establecimientos y roles diagnósticos: La distribución $\text{Binomial}(n_h, p)$ presupone que todo episodio de la cohorte comparte la misma probabilidad de fallecimiento ($p = 4,11\%$) o de estancia larga ($p_L = 31,11\%$). Los datos demuestran que este supuesto es falso: en la subcohorte principal $E_p$ la letalidad hospitalaria es $0,00\%$ ($0/108$) y la estancia larga es $5,56\%$, mientras que en la subcohorte secundaria $E_s$ suben a $5,69\%$ ($16/281$) y $40,93\%$. Como hospitales tales como Osorno ($93,6\%$ de casos $E_p$) y San José o Curicó ($100\%$ de casos $E_s$) atienden mezclas de casos opuestas, las probabilidades binomiales de cola de la Tabla 9 reflejan diferencias de composición diagnóstica y nunca calidad asistencial.
Incumplimiento del supuesto de estacionariedad y equidispersión en el modelo de Poisson diario: La distribución $\text{Poisson}(\lambda)$ exige que la tasa esperada de egresos sea constante en todos los días del año y que la varianza iguale a la media ($ID = 1$). En la red hospitalaria pública, las altas programadas se concentran de lunes a viernes ($\bar{D}_{\text{háb}} = 1,3244$) y disminuyen abruptamente los fines de semana ($0,56$ en sábado y $0,25$ en domingo), generando una mezcla de tasas que infla el índice de dispersión anual a $ID = 1,3857$ y subestima en más de $20$ días la frecuencia de jornadas con cero egresos.
Incumplimiento del supuesto de muestreo con reposición infinita en los modelos Geométrico, Binomial Negativo y Binomial de auditoría: Tanto la distribución Geométrica ($T$) como la Binomial Negativa ($T_5$) y la Binomial ($B$) asumen extracciones independientes con reposición sobre una población infinita. En una auditoría hospitalaria real sobre un conjunto finito de fichas ($N_M = 962$ fichas marcadas por $M$, o $N_h = 13$ fichas en el Hospital de Valdivia) sin reposición, las extracciones son dependientes y la varianza observada en los datos es sistemáticamente menor que la varianza teórica del modelo con reposición (Tabla 14 y Tabla 20). Además, la unidad de análisis es el episodio de egreso sin identificador longitudinal de persona, por lo que no es posible descontar posibles reingresos de un mismo individuo.

7.3 Respuesta directa a la jefatura: tamaño de revisión para garantizar 90% de probabilidad de hallazgo

Pregunta de la jefatura: ¿Cuántas fichas marcadas por $M$ debe revisar un auditor para tener al menos $90\%$ de probabilidad de encontrar al menos un caso de Catarata Senil (H25)?

Buscamos el menor número entero de ensayos $k \in \{1, 2, \dots\}$ tal que la función de distribución acumulada de $T \sim \text{Geométrica}(p = \text{VPP})$ cumpla: \[P(T \le k) = 1 - (1 - \text{VPP})^k \ge 0,90 \iff (1 - \text{VPP})^k \le 0,10 \iff k \ge \frac{\ln(0,10)}{\ln(1 - \text{VPP})}\]

k90_adj <- qgeom(0.90, prob = vpp_M) + 1
p_k90_adj_menos1 <- pgeom((k90_adj - 1) - 1, prob = vpp_M)
p_k90_adj <- pgeom(k90_adj - 1, prob = vpp_M)

k90_orig <- qgeom(0.90, prob = vpp_M_orig) + 1
p_k90_orig_menos1 <- pgeom((k90_orig - 1) - 1, prob = vpp_M_orig)
p_k90_orig <- pgeom(k90_orig - 1, prob = vpp_M_orig)

k90_azar <- qgeom(0.90, prob = p_azar) + 1
p_k90_azar <- pgeom(k90_azar - 1, prob = p_azar)

tabla_k90 <- tibble(
  Estrategia = c(
    "1. Con Marcador Ajustado M (Oftalmología y edad ≥ 65 años)",
    "2. Con Marcador Original Act. 2 M_orig (Edad ≥ 65 años e ingreso programado)",
    "3. Sin Marcador (Revisión al azar sobre todo el espacio muestral Ω)"
  ),
  Parametro_p = c(
    sprintf("VPP = 99 / 962 = %.6f (10,2911%%)", vpp_M),
    sprintf("VPP_orig = 110 / 57.303 = %.6f (0,1920%%)", vpp_M_orig),
    sprintf("P(E) = 389 / 873.718 = %.6f (0,0445%%)", p_azar)
  ),
  Formula_Exacta = c(
    sprintf("⌈ln(0,10) / ln(1 - 0,102911)⌉ = ⌈%.2f⌉ = %d", log(0.10) / log(1 - vpp_M), k90_adj),
    sprintf("⌈ln(0,10) / ln(1 - 0,001920)⌉ = ⌈%.2f⌉ = %d", log(0.10) / log(1 - vpp_M_orig), k90_orig),
    sprintf("⌈ln(0,10) / ln(1 - 0,000445)⌉ = ⌈%.2f⌉ = %d", log(0.10) / log(1 - p_azar), k90_azar)
  ),
  Prob_k_menos_1 = c(
    sprintf("P(T ≤ %d) = %.4f (%.2f%% < 90%%)", k90_adj - 1, p_k90_adj_menos1, p_k90_adj_menos1 * 100),
    sprintf("P(T ≤ %d) = %.4f (%.2f%% < 90%%)", k90_orig - 1, p_k90_orig_menos1, p_k90_orig_menos1 * 100),
    sprintf("P(T ≤ %d) = %.4f (%.2f%% < 90%%)", k90_azar - 1, pgeom((k90_azar - 1) - 1, prob = p_azar), pgeom((k90_azar - 1) - 1, prob = p_azar) * 100)
  ),
  Prob_k_minimo = c(
    sprintf("P(T ≤ %d) = %.4f (%.2f%% ≥ 90%%)", k90_adj, p_k90_adj, p_k90_adj * 100),
    sprintf("P(T ≤ %d) = %.4f (%.2f%% ≥ 90%%)", k90_orig, p_k90_orig, p_k90_orig * 100),
    sprintf("P(T ≤ %d) = %.4f (%.2f%% ≥ 90%%)", k90_azar, p_k90_azar, p_k90_azar * 100)
  )
)

kable(
  tabla_k90,
  col.names = c("Estrategia de auditoría", "Probabilidad de éxito por ficha (p)", "Cálculo analítico y qgeom(0.90, p) + 1", "Verificación en k - 1", "Verificación en k mínimo"),
  align = c("l", "r", "r", "r", "r"),
  caption = "Tabla 21: Cálculo y verificación del número mínimo de fichas k para alcanzar al menos 90% de probabilidad de encontrar un episodio de Catarata Senil (H25)."
) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = TRUE)

Tabla 21: Cálculo y verificación del número mínimo de fichas k para alcanzar al menos 90% de probabilidad de encontrar un episodio de Catarata Senil (H25).
Estrategia de auditoría	Probabilidad de éxito por ficha (p)	Cálculo analítico y qgeom(0.90, p) + 1	Verificación en k - 1	Verificación en k mínimo
Con Marcador Ajustado M (Oftalmología y edad ≥ 65 años)	VPP = 99 / 962 = 0.102911 (10,2911%)	⌈ln(0,10) / ln(1 - 0,102911)⌉ = ⌈21.20⌉ = 22	P(T ≤ 21) = 0.8978 (89.78% < 90%)	P(T ≤ 22) = 0.9083 (90.83% ≥ 90%)
Con Marcador Original Act. 2 M_orig (Edad ≥ 65 años e ingreso programado)	VPP_orig = 110 / 57.303 = 0.001920 (0,1920%)	⌈ln(0,10) / ln(1 - 0,001920)⌉ = ⌈1198.35⌉ = 1199	P(T ≤ 1198) = 0.8999 (89.99% < 90%)	P(T ≤ 1199) = 0.9001 (90.01% ≥ 90%)
Sin Marcador (Revisión al azar sobre todo el espacio muestral Ω)	P(E) = 389 / 873.718 = 0.000445 (0,0445%)	⌈ln(0,10) / ln(1 - 0,000445)⌉ = ⌈5170.60⌉ = 5171	P(T ≤ 5170) = 0.9000 (90.00% < 90%)	P(T ≤ 5171) = 0.9000 (90.00% ≥ 90%)

Respuesta ejecutiva para la jefatura (Tabla 21):

Utilizando el marcador ajustado $M$ (episodios de $\ge 65$ años atendidos o egresados por Oftalmología, cuyo $\text{VPP} = 10,2911\%$), el auditor debe revisar al menos $k = 22$ fichas marcadas por $M$ (qgeom(0.90, prob = 99/962) + 1 = 22). Con $21$ fichas la probabilidad acumulada es de $89,78\%$ (insuficiente), mientras que en la ficha $k = 22$ alcanza exactamente $P(T \le 22) = 90,83\% \ge 0,90$.
Si en cambio se utilizara el marcador administrativo amplio de la Actividad 2 ($M_{\text{orig}}$: $\ge 65$ años con ingreso programado sin filtrar especialidad, $\text{VPP}_{\text{orig}} = 0,1920\%$), el auditor tendría que revisar al menos $k = 1.199$ fichas ($P(T \le 1.199) = 90,01\%$), y si buscara al azar sobre toda la base $\Omega$ sin marcador ($P(E) = 0,0445\%$) necesitaría revisar $k = 5.171$ fichas. Por consiguiente, se recomienda adoptar operativamente el marcador ajustado $M$, que garantiza el $90\%$ de probabilidad de detección con solo $22$ fichas auditadas.

8 Declaración de Uso de Inteligencia Artificial y Verificación Técnica

En conformidad con las normas del curso, se declara el uso de herramientas de Inteligencia Artificial generativa (Claude / Gemini AI en entorno Antigravity) como apoyo para: (1) optimizar la lectura vectorizada de la base GRD_PUBLICO_2024.txt y el cruce relacional con Tablas_maestras_bases_GRD.xlsx, (2) estructurar el código reproducible de las familias de funciones d, p, q y r de R cuidando la parametrización por fracasos previos de dgeom/dnbinom, y (3) formatear las tablas y gráficos del documento RMarkdown.

Verificación de resultados: Todos los cálculos estadísticos fueron verificados ejecutando el código directamente sobre los $1.085.813$ registros de la base oficial GRD 2024 de FONASA, comprobando manualmente la coincidencia entre momentos teóricos y empíricos (incluyendo la corrección de Bessel $n/(n-1)$ en la varianza de $X$, el factor de población finita $(N_h-n)/(N_h-1)$ en la Hipergeométrica, la verificación manual de $P(T=1)$ y $P(T=2)$ en la Geométrica y la reproducción exacta de las simulaciones con set.seed(2026)).

Actividad 3: Variables Aleatorias y Distribuciones Discretas en Episodios GRD de Catarata Senil (H25)

Estadística Computacional I · Departamento de Ingeniería Informática · USACH

David Capetillo

Octubre de 2026

1 Resultado 01: Una variable aleatoria discreta construida desde los datos

1.1 Correcciones aplicadas a partir de las Actividades 1 y 2

1.2 Definición de variables aleatorias discretas, soporte y modelos propuestos

1.3 Función de probabilidad empírica y acumulada de \(X\)

1.4 Esperanza, varianza, mediana, moda y probabilidad de cola de \(X\)

1.5 Comparación de la función de probabilidad de \(X\) según rol diagnóstico (\(E_p\) vs \(E_s\))

2 Resultado 02: Distribuciones Bernoulli y Binomial

2.1 Parámetros de Bernoulli (\(Y\) para fallecimiento y \(Y_L\) para estancia larga)

2.2 Modelo Binomial en una auditoría de \(n = 20\) episodios (\(S\) y \(S_L\))