Versión de impresión: se muestran todas las secciones en una sola página.

# ---- rutas relativas ----
# La base GRD original no se entrega (es pública y pesa ~570 MB). Para reproducir el informe,
# se deja en la subcarpeta GRD_PUBLICO_2024/ junto a este .Rmd, igual que en las Actividades 1 y 2.
# Como respaldo se busca también en la carpeta de la Actividad 1 (../informe-1/).
buscar_archivo <- function(rel) {
  candidatos <- c(rel, file.path("..", "informe-1", rel))
  f <- candidatos[file.exists(candidatos)][1]
  if (is.na(f)) stop("No se encontró el archivo ", rel, " junto al .Rmd")
  f
}
ruta_grd <- buscar_archivo(file.path("GRD_PUBLICO_2024", "GRD_PUBLICO_2024.txt"))

cols_diag <- paste0("DIAGNOSTICO", 1:35)
cols_usar <- c("COD_HOSPITAL", "ID_BENEFICIARIO", "FECHA_NACIMIENTO", "TIPO_INGRESO",
               "ESPECIALIDAD_MEDICA", "TIPO_ACTIVIDAD", "FECHA_INGRESO", "FECHAALTA", "TIPOALTA",
               cols_diag)

base <- fread(ruta_grd, sep = "|", encoding = "Latin-1", quote = "", na.strings = "",
              select = cols_usar, colClasses = list(character = c("COD_HOSPITAL", "ID_BENEFICIARIO")))

# ---- Ω: igual que en la Actividad 2 (sin CMA y con fechas válidas en todos los egresos) ----
n0    <- nrow(base)
n_cma <- base[TIPO_ACTIVIDAD == "CIRUGÍA MAYOR AMBULATORIA (CMA)", .N]
hosp_all <- base[TIPO_ACTIVIDAD == "HOSPITALIZACIÓN"]
rm(base); invisible(gc())
n1 <- nrow(hosp_all)
fecha_valida <- function(x) grepl("^[0-9]{4}-[0-9]{2}-[0-9]{2}$", x)
ok_fechas <- fecha_valida(hosp_all$FECHA_INGRESO) & fecha_valida(hosp_all$FECHA_NACIMIENTO) &
             fecha_valida(hosp_all$FECHAALTA)
n_malas_omega <- sum(!ok_fechas)
omega <- hosp_all[ok_fechas]
rm(hosp_all, ok_fechas); invisible(gc())
N <- nrow(omega)

# ---------------- Eventos de la Actividad 2 ----------------
omega[, E  := Reduce(`|`, lapply(cols_diag, function(cc) grepl("^I50", omega[[cc]])))]
omega[, Ep := grepl("^I50", DIAGNOSTICO1)]            # I50 como diagnóstico principal
omega[, Es := E & !Ep]                                # I50 solo como diagnóstico secundario

# ---------------- Variable X: número de diagnósticos registrados ----------------
# Un campo cuenta si no está vacío. El texto "DESCONOCIDO" no es un código CIE-10, por lo que
# no se cuenta como diagnóstico (X_bruto lo cuenta, solo para documentar el efecto de la decisión).
diag_lleno  <- function(x) !is.na(x) & trimws(x) != ""
omega[, X_bruto := Reduce(`+`, lapply(cols_diag, function(cc) diag_lleno(omega[[cc]])))]
omega[, X := Reduce(`+`, lapply(cols_diag, function(cc) diag_lleno(omega[[cc]]) & omega[[cc]] != "DESCONOCIDO"))]
omega[, (cols_diag) := NULL]; invisible(gc())

omega[, Fal := TIPOALTA == "FALLECIDO"]               # evento F (se evita el nombre F, que en R es FALSE)
omega[, U   := TIPO_INGRESO == "URGENCIA"]            # evento U
omega[, estancia := as.numeric(as.Date(FECHAALTA) - as.Date(FECHA_INGRESO))]
corte_L <- unname(quantile(omega$estancia, 0.75))     # percentil 75 de la estancia en todo Ω
omega[, L := estancia > corte_L]                      # evento L: estancia larga
omega[, edad := as.numeric(as.Date(FECHA_INGRESO) - as.Date(FECHA_NACIMIENTO)) / 365.25]

# marcador M de la Actividad 2 (sin cambios)
esp_M <- c("MEDICINA INTERNA", "CARDIOLOGÍA", "GERIATRÍA")
omega[, M := ESPECIALIDAD_MEDICA %in% esp_M & edad >= 65]

# ---------------- cohorte E ----------------
coh <- omega[E == TRUE]
nE  <- nrow(coh); pE <- nE / N
nEF <- sum(coh$Fal); p <- nEF / nE                     # p = P(F | E)
nEL <- sum(coh$L);   pL <- nEL / nE                    # P(L | E)
nM  <- sum(omega$M); nME <- sum(omega$M & omega$E)
vpp <- nME / nM                                       # VPP = P(E | M)
sens <- nME / nE; espec <- sum(!omega$M & !omega$E) / (N - nE)
n_desc_E <- sum(coh$X_bruto != coh$X)                 # episodios de E con algún "DESCONOCIDO"
med_omega <- median(omega$X)

# ---------------- serie diaria D (366 días de 2024, los días sin egresos quedan en 0) ----------------
dias <- data.table(fecha = seq(as.Date("2024-01-01"), as.Date("2024-12-31"), by = "day"))
cnt_dia <- coh[, .(D = .N), by = .(fecha = as.Date(FECHAALTA))]
dias <- merge(dias, cnt_dia, by = "fecha", all.x = TRUE)
dias[is.na(D), D := 0L]
n_fuera_2024 <- sum(!(as.Date(coh$FECHAALTA) %in% dias$fecha))
lambda <- mean(dias$D)                                 # = n(E) / 366

# ---------------- establecimientos ----------------
hosp <- suppressMessages(as.data.table(read_excel("Tablas_maestras_GRD.xlsx", sheet = "Hospitales")))
setnames(hosp, c("cod_hospital", "nombre_hospital"))
hosp[, cod_hospital := as.character(as.integer(cod_hospital))]
est <- coh[, .(N_h = .N, K_F = sum(Fal), K_L = sum(L)), by = COD_HOSPITAL]
est <- merge(est, hosp, by.x = "COD_HOSPITAL", by.y = "cod_hospital", all.x = TRUE)
est[is.na(nombre_hospital), nombre_hospital := paste0("Cód. ", COD_HOSPITAL)]
n_est <- nrow(est)

# ---------------- beneficiarios con más de un episodio (para discutir la independencia) ----------------
benef <- coh[!is.na(ID_BENEFICIARIO) & ID_BENEFICIARIO != "", .N, by = ID_BENEFICIARIO]
n_benef <- nrow(benef)
pct_ep_repetidos <- benef[N >= 2, sum(N)] / sum(benef$N)
max_ep_benef <- max(benef$N)

1 Introducción

1.1 Contexto y objetivo

Este informe continúa las Actividades 1 y 2 en el rol de analista de la Unidad de Análisis Hospitalario. La Actividad 1 dimensionó y caracterizó los episodios con insuficiencia cardíaca (CIE-10 I50) en los egresos GRD 2024. La Actividad 2 expresó esas características como eventos y calculó sus probabilidades. Ahora la jefatura quiere anticipar conteos: cuántos fallecimientos esperar en una auditoría de 20 episodios, cuántas fichas revisar hasta encontrar un episodio con insuficiencia cardíaca o cuántos egresos de la cohorte esperar en un día cualquiera.

Para responder, los eventos de la Actividad 2 se convierten en variables aleatorias discretas, y cada una se modela con una distribución: Bernoulli, binomial, hipergeométrica, geométrica, binomial negativa y Poisson. Cada modelo se compara con lo que muestran los datos y con simulación en R. Los parámetros se obtienen con proporciones y promedios de la base. No se hacen pruebas de hipótesis ni ajustes por máxima verosimilitud, y la comparación entre modelo y datos es descriptiva.

Qué es una variable aleatoria en este informe Una variable aleatoria asigna un número a cada resultado de un experimento bien definido: elegir un episodio al azar, auditar 20 episodios, revisar fichas una por una o contar los egresos de un día. Cada valor se refiere a episodios de hospitalización, no a personas: una misma persona puede aportar varios episodios, y ninguna cifra describe riesgo de personas, incidencia ni prevalencia.

1.2 Declaración de uso de inteligencia artificial

Se utilizó Claude Code (Anthropic, modelo Claude Opus 5.5) como apoyo para: (i) escribir y depurar el código R de construcción de variables, cálculo de distribuciones y simulación; (ii) redactar borradores del texto interpretativo; y (iii) revisar el cumplimiento de la pauta y de la retroalimentación de las Actividades 1 y 2. Los resultados se verificaron así: (a) ninguna cifra del texto está escrita a mano: todas se insertan desde el código y se recalculan en cada compilación; (b) la convención de R para la geométrica y la binomial negativa se verifica comparando dgeom y dnbinom con el cálculo a mano de la fórmula, con una columna de diferencia (sección 2.5); (c) la esperanza y la varianza de X se calculan desde la función de probabilidad y se contrastan con mean() y var(); (d) cada modelo se compara con los datos (frecuencias observadas contra esperadas, remuestreo de episodios reales), no solo se calcula; (e) todas las simulaciones usan semilla fija y se comprobó que dos compilaciones sucesivas entregan los mismos números; y (f) las frases que comparan valores se generan a partir de los propios valores calculados, para que el texto no pueda contradecir a su tabla. La responsabilidad técnica y de contenido del informe es del autor.

2 Datos, correcciones y variables

RESULTADO 01 · PARTE 1

Espacio muestral y cohorte Ω = egresos hospitalarios GRD 2024, sin CMA y con fechas de ingreso, nacimiento y alta válidas: N = 873.719 episodios. La cohorte E (insuficiencia cardíaca en cualquiera de los 35 campos de diagnóstico) tiene n(E) = 52.312 episodios. Son exactamente Ω y E de la Actividad 2.

2.1 Fuente de datos y tablas auxiliares

Base: GRD_PUBLICO_2024.txt, base pública de egresos hospitalarios y cirugía mayor ambulatoria del mecanismo GRD, Fondo Nacional de Salud (FONASA), año 2024. Texto plano delimitado por |, codificación Latin-1. Descargada desde el portal de datos abiertos de FONASA el 21 de agosto de 2026. No se entrega por su tamaño; para reproducir el informe se ubica en la subcarpeta GRD_PUBLICO_2024/.
Tablas_maestras_GRD.xlsx (se entrega en el zip): hoja Hospitales, para el nombre de cada establecimiento.
CIE-10.xlsx (se entrega en el zip): catálogo CIE-10, para verificar los códigos I50.

cie10 <- suppressMessages(as.data.table(read_excel("CIE-10.xlsx", sheet = "CIE 10")))
tabla_html(cie10[grepl("^I50", Código), .(Código, Descripción)],
           caption = "Códigos CIE-10 que definen la cohorte E (fuente: CIE-10.xlsx)",
           label = "cie10")

Tabla 1. Códigos CIE-10 que definen la cohorte E (fuente: CIE-10.xlsx)
Código	Descripción
I50	Insuficiencia cardíaca
I50.0	Insuficiencia cardíaca congestiva
I50.1	Insuficiencia ventricular izquierda
I50.9	Insuficiencia cardíaca, no especificada

La cohorte se define con el patrón ^I50, que captura I50 y sus subcategorías (Tabla 1), buscado en los 35 campos de diagnóstico.

2.2 Correcciones a partir de la retroalimentación de las Actividades 1 y 2

cambios <- data.table(
  `N.º` = 1:5,
  `Observación` = c(
    "Actividad 2: las diferencias de las tablas de verificación aparecían en notación científica (0.0e+00) en un informe que usa coma decimal.",
    "Actividad 2: el zip no traía las tablas auxiliares que el texto decía entregar.",
    "Actividad 2: el marcador M debe revisarse si su valor predictivo positivo es cero o casi cero (instrucción de la Actividad 3).",
    "Actividad 1: hablar de episodios y no de personas, nombrar el estadístico que se compara (mediana o promedio), no proponer hipótesis explicativas y reescribir la frase sobre la fecha de extracción.",
    "Actividad 1: la depuración de fechas no válidas se aplicaba solo a los episodios con I50 (corrección propia hecha en la Actividad 2)."),
  `Cambio aplicado en la Actividad 3` = c(
    "Todas las cifras usan decimales fijos con coma decimal. Las diferencias de verificación se muestran con 6 decimales y las probabilidades pequeñas con 3 cifras significativas (por ejemplo, P(T > 50) en el Resultado 04), sin notación científica ni redondeos a cero.",
    "El zip incluye el .Rmd, el HTML, Tablas_maestras_GRD.xlsx, CIE-10.xlsx y el enlace al HTML publicado. Se revisó el contenido del zip antes de entregarlo.",
    paste0("No se requiere ajuste: el VPP de M es ", fpct(vpp, 1), " (", fnum(nME), " / ", fnum(nM), "), lejos de cero. Se usa el mismo M de la Actividad 2."),
    "Se mantiene lo aplicado en la Actividad 2: todo el texto habla de episodios, se nombra el estadístico y las diferencias entre grupos se describen sin atribuirles causas.",
    "Se mantiene: la regla de fechas válidas se aplica a todo Ω."),
  `Efecto sobre Ω o la cohorte` = c("Ninguno (formato).", "Ninguno (entrega).", "Ninguno.",
                                     "Ninguno (redacción).",
                                     paste0("Ω = ", fnum(N), " egresos (", fnum(n_malas_omega), " excluidos por fecha no válida), igual que en la Actividad 2."))
)
tabla_html(cambios, caption = "Correcciones aplicadas a partir de la retroalimentación de las Actividades 1 y 2",
           label = "cambios") |>
  column_spec(2:3, width = "40%") |> kable_styling(full_width = TRUE)

Tabla 2. Correcciones aplicadas a partir de la retroalimentación de las Actividades 1 y 2
N.º	Observación	Cambio aplicado en la Actividad 3	Efecto sobre Ω o la cohorte
1	Actividad 2: las diferencias de las tablas de verificación aparecían en notación científica (0.0e+00) en un informe que usa coma decimal.	Todas las cifras usan decimales fijos con coma decimal. Las diferencias de verificación se muestran con 6 decimales y las probabilidades pequeñas con 3 cifras significativas (por ejemplo, P(T > 50) en el Resultado 04), sin notación científica ni redondeos a cero.	Ninguno (formato).
2	Actividad 2: el zip no traía las tablas auxiliares que el texto decía entregar.	El zip incluye el .Rmd, el HTML, Tablas_maestras_GRD.xlsx, CIE-10.xlsx y el enlace al HTML publicado. Se revisó el contenido del zip antes de entregarlo.	Ninguno (entrega).
3	Actividad 2: el marcador M debe revisarse si su valor predictivo positivo es cero o casi cero (instrucción de la Actividad 3).	No se requiere ajuste: el VPP de M es 26,6% (26.127 / 98.239), lejos de cero. Se usa el mismo M de la Actividad 2.	Ninguno.
4	Actividad 1: hablar de episodios y no de personas, nombrar el estadístico que se compara (mediana o promedio), no proponer hipótesis explicativas y reescribir la frase sobre la fecha de extracción.	Se mantiene lo aplicado en la Actividad 2: todo el texto habla de episodios, se nombra el estadístico y las diferencias entre grupos se describen sin atribuirles causas.	Ninguno (redacción).
5	Actividad 1: la depuración de fechas no válidas se aplicaba solo a los episodios con I50 (corrección propia hecha en la Actividad 2).	Se mantiene: la regla de fechas válidas se aplica a todo Ω.	Ω = 873.719 egresos (53 excluidos por fecha no válida), igual que en la Actividad 2.

No hubo cambios en la cohorte: el tamaño final es de 52.312 episodios, el mismo de las Actividades 1 y 2 (Tabla 2). Las correcciones pendientes de la Actividad 2 se refieren al formato de los números y al contenido de la entrega.

2.3 Tabla de variables aleatorias

Cada variable se define por el experimento que la genera. "Al azar" significa que cada episodio del conjunto indicado tiene la misma probabilidad de ser elegido.

codigo <- function(x) paste0("<code>", gsub("`", "&#96;", htmltools::htmlEscape(x, attribute = TRUE)), "</code>")
vars <- data.table(
  Variable = c("X", "Y", "S", "S<sub>L</sub>", "H", "T", "T<sub>0</sub>", "T<sub>5</sub>", "D", "B"),
  Experimento = c(
    "Elegir un episodio al azar de la cohorte E y contar sus diagnósticos registrados (campos de diagnóstico no vacíos).",
    "Elegir un episodio al azar de E; Y vale 1 si termina en fallecimiento y 0 si no.",
    "Elegir 20 episodios de E al azar con reposición y contar cuántos terminan en fallecimiento.",
    "Elegir 20 episodios de E al azar con reposición y contar cuántos tienen estancia larga (L).",
    "Auditar 10 episodios sin reposición dentro de un establecimiento (con N episodios de E y K fallecimientos) y contar los fallecimientos.",
    "Revisar, una por una y al azar, fichas de episodios marcados por M hasta encontrar la primera con I50; T es el número de fichas revisadas.",
    "Lo mismo que T, pero revisando episodios al azar de Ω, sin usar el marcador.",
    "Revisar fichas de episodios marcados por M hasta encontrar 5 con I50; T₅ es el número de fichas revisadas.",
    "Elegir un día de 2024 y contar los episodios de E con fecha de egreso ese día.",
    "Elegir 1.000 episodios de Ω al azar y contar cuántos tienen I50."),
  Soporte = c("{1, 2, ..., 35}", "{0, 1}", "{0, 1, ..., 20}", "{0, 1, ..., 20}",
              "{max(0, 10 − (N − K)), ..., min(10, K)}", "{1, 2, 3, ...}", "{1, 2, 3, ...}",
              "{5, 6, 7, ...}", "{0, 1, 2, ...}", "{0, 1, ..., 1.000}"),
  `Modelo propuesto` = c("Distribución empírica", "Bernoulli(p)", "Binomial(20, p)",
                         "Binomial(20, P(L | E))", "Hipergeométrica(N, K, 10)", "Geométrica(VPP)",
                         "Geométrica(P(E))", "Binomial negativa(5, VPP)", "Poisson(λ)",
                         "Binomial(1.000, P(E)) ≈ Poisson(1.000 · P(E))"),
  `Cálculo en R` = codigo(c("table(coh$X) / nrow(coh)", "dbinom(k, 1, p)", "dbinom(k, 20, p)",
                            "dbinom(k, 20, pL)", "dhyper(k, K, N - K, 10)", "dgeom(k - 1, vpp)",
                            "dgeom(k - 1, pE)", "dnbinom(k - 5, 5, vpp)", "dpois(k, lambda)",
                            "dbinom(k, 1000, pE); dpois(k, 1000 * pE)"))
)
tabla_html(vars, caption = "Variables aleatorias: experimento, soporte, modelo propuesto y cálculo en R",
           label = "vars", escape = FALSE) |>
  column_spec(1, bold = TRUE) |> column_spec(2, width = "38%") |> kable_styling(full_width = TRUE)

Tabla 3. Variables aleatorias: experimento, soporte, modelo propuesto y cálculo en R
Variable	Experimento	Soporte	Modelo propuesto	Cálculo en R
X	Elegir un episodio al azar de la cohorte E y contar sus diagnósticos registrados (campos de diagnóstico no vacíos).	{1, 2, ..., 35}	Distribución empírica	`table(coh$X) / nrow(coh)`
Y	Elegir un episodio al azar de E; Y vale 1 si termina en fallecimiento y 0 si no.	{0, 1}	Bernoulli(p)	`dbinom(k, 1, p)`
S	Elegir 20 episodios de E al azar con reposición y contar cuántos terminan en fallecimiento.	{0, 1, ..., 20}	Binomial(20, p)	`dbinom(k, 20, p)`
S_L	Elegir 20 episodios de E al azar con reposición y contar cuántos tienen estancia larga (L).	{0, 1, ..., 20}	Binomial(20, P(L \| E))	`dbinom(k, 20, pL)`
H	Auditar 10 episodios sin reposición dentro de un establecimiento (con N episodios de E y K fallecimientos) y contar los fallecimientos.	{max(0, 10 − (N − K)), ..., min(10, K)}	Hipergeométrica(N, K, 10)	`dhyper(k, K, N - K, 10)`
T	Revisar, una por una y al azar, fichas de episodios marcados por M hasta encontrar la primera con I50; T es el número de fichas revisadas.	{1, 2, 3, ...}	Geométrica(VPP)	`dgeom(k - 1, vpp)`
T₀	Lo mismo que T, pero revisando episodios al azar de Ω, sin usar el marcador.	{1, 2, 3, ...}	Geométrica(P(E))	`dgeom(k - 1, pE)`
T₅	Revisar fichas de episodios marcados por M hasta encontrar 5 con I50; T₅ es el número de fichas revisadas.	{5, 6, 7, ...}	Binomial negativa(5, VPP)	`dnbinom(k - 5, 5, vpp)`
D	Elegir un día de 2024 y contar los episodios de E con fecha de egreso ese día.	{0, 1, 2, ...}	Poisson(λ)	`dpois(k, lambda)`
B	Elegir 1.000 episodios de Ω al azar y contar cuántos tienen I50.	{0, 1, ..., 1.000}	Binomial(1.000, P(E)) ≈ Poisson(1.000 · P(E))	`dbinom(k, 1000, pE); dpois(k, 1000 * pE)`

En la notación de la pauta, T₅ corresponde a la segunda variable T (binomial negativa). Se agregan S_L (pedida en el Resultado 02), T₀ (pedida en el Resultado 04) y B (pedida en el Resultado 05) para que todas las variables usadas queden definidas en un solo lugar (Tabla 3).

Decisión sobre X. En 306 episodios de la cohorte algún campo de diagnóstico contiene el texto "DESCONOCIDO", que no es un código CIE-10. Ese texto no se cuenta como diagnóstico. Con esta regla, X cambia en 306 de los 52.312 episodios (0,58%), y E[X] cambia de 12,4293 a 12,4234.

2.4 Parámetros estimados desde la base

lambda_txt <- paste0(fnum(nE), " / 366")
par_tab <- data.table(
  `Parámetro` = c("p = P(F | E)", "P(L | E)", "VPP = P(E | M)", "P(E)", "λ (por día)",
                  "N y K de cada establecimiento"),
  `Se usa en` = c("Y, S", "S<sub>L</sub>", "T, T<sub>5</sub>", "T<sub>0</sub>, B", "D", "H"),
  `Cálculo` = c(paste0("n(F ∩ E) / n(E) = ", fnum(nEF), " / ", fnum(nE)),
                paste0("n(L ∩ E) / n(E) = ", fnum(nEL), " / ", fnum(nE)),
                paste0("n(M ∩ E) / n(M) = ", fnum(nME), " / ", fnum(nM)),
                paste0("n(E) / N = ", fnum(nE), " / ", fnum(N)),
                paste0("n(E) / días de 2024 = ", lambda_txt),
                "N = episodios de E en el establecimiento; K = fallecimientos (o estancias largas) entre ellos"),
  Valor = c(paste0(f4(p), " (", fpct(p), ")"), paste0(f4(pL), " (", fpct(pL), ")"),
            paste0(f4(vpp), " (", fpct(vpp), ")"), paste0(f4(pE), " (", fpct(pE), ")"),
            fnum(lambda, 2), "Ver Resultado 03"),
  `Total sobre el que se obtuvo` = c(paste0("n(E) = ", fnum(nE), " episodios"),
                                     paste0("n(E) = ", fnum(nE), " episodios"),
                                     paste0("n(M) = ", fnum(nM), " episodios marcados"),
                                     paste0("N = ", fnum(N), " episodios de Ω"),
                                     "366 días (incluye días sin egresos)",
                                     paste0(n_est, " establecimientos con episodios de E")),
  Origen = c("Actividad 2 (P(F | E))", "Actividad 2 (evento L: estancia > percentil 75 de Ω)",
             "Actividad 2 (marcador M)", "Actividad 2 (P(E))", "Nuevo (serie diaria)", "Nuevo (por establecimiento)")
)
tabla_html(par_tab, caption = "Parámetros de los modelos: cálculo, valor, total sobre el que se obtuvieron y origen",
           label = "param", escape = FALSE) |> column_spec(1, bold = TRUE)

Tabla 4. Parámetros de los modelos: cálculo, valor, total sobre el que se obtuvieron y origen
Parámetro	Se usa en	Cálculo	Valor	Total sobre el que se obtuvo	Origen
p = P(F \| E)	Y, S	n(F ∩ E) / n(E) = 5.231 / 52.312	0,1000 (10,00%)	n(E) = 52.312 episodios	Actividad 2 (P(F \| E))
P(L \| E)	S_L	n(L ∩ E) / n(E) = 23.382 / 52.312	0,4470 (44,70%)	n(E) = 52.312 episodios	Actividad 2 (evento L: estancia > percentil 75 de Ω)
VPP = P(E \| M)	T, T₅	n(M ∩ E) / n(M) = 26.127 / 98.239	0,2660 (26,60%)	n(M) = 98.239 episodios marcados	Actividad 2 (marcador M)
P(E)	T₀, B	n(E) / N = 52.312 / 873.719	0,0599 (5,99%)	N = 873.719 episodios de Ω	Actividad 2 (P(E))
λ (por día)	D	n(E) / días de 2024 = 52.312 / 366	142,93	366 días (incluye días sin egresos)	Nuevo (serie diaria)
N y K de cada establecimiento	H	N = episodios de E en el establecimiento; K = fallecimientos (o estancias largas) entre ellos	Ver Resultado 03	72 establecimientos con episodios de E	Nuevo (por establecimiento)

Interpretación (Tabla 4): todos los parámetros son proporciones o promedios de la base. p = 10,00% es la proporción de episodios de la cohorte que terminan en fallecimiento, y P(L | E) = 44,70% la de episodios con estancia mayor que 7 días (el percentil 75 de Ω). Que P(L | E) sea mucho mayor que p hace que S_L tenga una forma muy distinta de S (Resultado 02). El VPP = 26,60% es 4,4 veces P(E), y esa diferencia es la que determina cuántas fichas ahorra el marcador (Resultado 04).

2.5 Convención de R para la geométrica y la binomial negativa

Convención declarada En este informe T y T₅ cuentan ensayos (fichas revisadas, incluida la que encuentra el caso). Las funciones dgeom, pgeom, qgeom y rgeom de R cuentan fracasos antes del primer éxito, y dnbinom cuenta fracasos antes del r-ésimo éxito. Por eso: P(T = k) = dgeom(k - 1, p), P(T ≤ k) = pgeom(k - 1, p), el cuantil de T es qgeom(q, p) + 1 y una simulación de T es rgeom(n, p) + 1. Para T₅: P(T₅ = k) = dnbinom(k - 5, size = 5, prob = p) y una simulación es rnbinom(n, 5, p) + 5.

k_g <- 3; k_nb <- 7
a_mano_g  <- (1 - vpp)^(k_g - 1) * vpp
a_mano_nb <- choose(k_nb - 1, 5 - 1) * vpp^5 * (1 - vpp)^(k_nb - 5)
conv <- data.table(
  Probabilidad = c(paste0("P(T = ", k_g, ")"), paste0("P(T₅ = ", k_nb, ")")),
  `Fórmula a mano` = c(paste0("(1 − VPP)^", k_g - 1, " · VPP = ", f4(1 - vpp), "^", k_g - 1, " × ", f4(vpp), " = ", fnum(a_mano_g, 6)),
                       paste0("C(", k_nb - 1, ", 4) · VPP^5 · (1 − VPP)^", k_nb - 5, " = ", choose(k_nb - 1, 4), " × ", f4(vpp), "^5 × ", f4(1 - vpp), "^", k_nb - 5, " = ", fnum(a_mano_nb, 6))),
  `Función de R` = codigo(c(paste0("dgeom(", k_g - 1, ", vpp)"), paste0("dnbinom(", k_nb - 5, ", size = 5, prob = vpp)"))),
  `Valor en R` = fnum(c(dgeom(k_g - 1, vpp), dnbinom(k_nb - 5, 5, vpp)), 6),
  `Diferencia` = fdif(c(a_mano_g - dgeom(k_g - 1, vpp), a_mano_nb - dnbinom(k_nb - 5, 5, vpp))),
  `Error frecuente` = codigo(c(paste0("dgeom(", k_g, ", vpp) = ", fnum(dgeom(k_g, vpp), 6)),
                               paste0("dnbinom(", k_nb, ", 5, vpp) = ", fnum(dnbinom(k_nb, 5, vpp), 6))))
)
tabla_html(conv, caption = "Verificación a mano de la convención de R (ensayos = fracasos + éxitos)", label = "conv",
           escape = FALSE)

Tabla 5. Verificación a mano de la convención de R (ensayos = fracasos + éxitos)
Probabilidad	Fórmula a mano	Función de R	Valor en R	Diferencia	Error frecuente
P(T = 3)	(1 − VPP)^2 · VPP = 0,7340^2 × 0,2660 = 0,143302	`dgeom(2, vpp)`	0,143302	0,000000	`dgeom(3, vpp) = 0,105190`
P(T₅ = 7)	C(6, 4) · VPP^5 · (1 − VPP)^2 = 15 × 0,2660^5 × 0,7340^2 = 0,010754	`dnbinom(2, size = 5, prob = vpp)`	0,010754	0,000000	`dnbinom(7, 5, vpp) = 0,050421`

La fórmula a mano coincide con dgeom(k - 1, ·) y con dnbinom(k - 5, 5, ·) (Tabla 5, diferencia 0). La última columna muestra el error de pasar k directamente: dgeom(3, vpp) entrega P(T = 4), no P(T = 3).

3 Variable aleatoria empírica: número de diagnósticos

RESULTADO 01 · PARTE 2

X es el número de diagnósticos registrados en un episodio elegido al azar de la cohorte. Su función de probabilidad se obtiene directamente de las frecuencias relativas: p(x) = n(X = x) / n(E), con n(E) = 52.312.

3.1 Función de probabilidad y función de distribución acumulada

fpX <- coh[, .(n = .N), by = X][order(X)]
fpX[, px := n / nE][, Fx := cumsum(px)]
moda_X <- fpX[which.max(px), X]
mediana_X <- fpX[Fx >= 0.5][1, X]                 # menor x con F(x) ≥ 0,5
k_X <- 2 * med_omega                              # el doble de la mediana de Ω
tabla_html(fpX[, .(`x (diagnósticos)` = X, `n(X = x)` = fnum(n), `p(x) = n / n(E)` = fprob(px),
                   `F(x) = P(X ≤ x)` = fprob_cdf(Fx))],
           caption = paste0("Función de probabilidad empírica de X y su función de distribución acumulada (n(E) = ", fnum(nE), ")"),
           label = "fpX") |>
  row_spec(which(fpX$X == moda_X), bold = TRUE) |>
  scroll_box(height = "380px")

Tabla 6. Función de probabilidad empírica de X y su función de distribución acumulada (n(E) = 52.312)
x (diagnósticos)	n(X = x)	p(x) = n / n(E)	F(x) = P(X ≤ x)
1	17	0,0325%	0,0325%
2	95	0,18%	0,21%
3	355	0,68%	0,89%
4	884	1,69%	2,58%
5	1.700	3,25%	5,83%
6	2.768	5,29%	11,12%
7	3.555	6,80%	17,92%
8	4.194	8,02%	25,94%
9	4.605	8,80%	34,74%
10	4.477	8,56%	43,30%
11	4.314	8,25%	51,54%
12	3.966	7,58%	59,13%
13	3.365	6,43%	65,56%
14	3.026	5,78%	71,34%
15	2.582	4,94%	76,28%
16	2.139	4,09%	80,37%
17	1.776	3,40%	83,76%
18	1.480	2,83%	86,59%
19	1.238	2,37%	88,96%
20	1.005	1,92%	90,88%
21	829	1,58%	92,46%
22	685	1,31%	93,77%
23	568	1,09%	94,86%
24	476	0,91%	95,77%
25	412	0,79%	96,56%
26	330	0,63%	97,19%
27	261	0,50%	97,69%
28	234	0,45%	98,13%
29	168	0,32%	98,46%
30	144	0,28%	98,73%
31	125	0,24%	98,97%
32	105	0,20%	99,17%
33	95	0,18%	99,35%
34	109	0,21%	99,56%
35	230	0,44%	100%

EX_tmp <- sum(fpX$X * fpX$px)
ggplot(fpX, aes(x = X, y = px, fill = X == moda_X)) +
  geom_col(width = 0.8, show.legend = FALSE) +
  geom_vline(xintercept = EX_tmp, linetype = "dashed", color = paleta[2], linewidth = 0.9) +
  annotate("text", x = EX_tmp + 0.5, y = max(fpX$px) * 0.98, hjust = 0, color = paleta[2], size = 3.6,
           label = paste0("E[X] = ", fnum(EX_tmp, 2))) +
  scale_fill_manual(values = c(`TRUE` = paleta[3], `FALSE` = paleta[1])) +
  scale_x_continuous(breaks = seq(0, 35, 5)) +
  scale_y_continuous(labels = lab_pct) +
  labs(title = "Número de diagnósticos registrados en los episodios con insuficiencia cardíaca",
       subtitle = paste0("Cohorte E, n = ", fnum(nE), " episodios"), x = "x: diagnósticos registrados", y = "p(x)")

Figura 1. Función de probabilidad empírica de X (barras) en la cohorte E. La línea vertical marca E[X] y la barra oscura la moda.

ggplot(fpX, aes(x = X, y = Fx)) +
  geom_step(color = paleta[3], linewidth = 1) + geom_point(color = paleta[3], size = 1.6) +
  geom_hline(yintercept = c(0.5, 0.9), linetype = "dotted", color = "grey40") +
  scale_x_continuous(breaks = seq(0, 35, 5)) +
  scale_y_continuous(labels = lab_pct, breaks = seq(0, 1, 0.1)) +
  labs(title = "Función de distribución acumulada de X", x = "x: diagnósticos registrados", y = "F(x) = P(X ≤ x)")

Figura 2. Función de distribución acumulada de X. Las líneas punteadas marcan F(x) = 0,5 (mediana) y F(x) = 0,9.

Interpretación (Tabla 6, Figura 1 y Figura 2): la distribución de X es asimétrica hacia la derecha: la moda es 9 diagnósticos (p(9) = 8,80%), la mediana es 11 y la cola se extiende hasta 35. El 90% de los episodios tiene a lo más 20 diagnósticos (Figura 2). El valor 35 concentra 230 episodios (0,44%), más que 33 o 34: la base tiene solo 35 campos de diagnóstico, por lo que X = 35 significa "35 o más" y el soporte está truncado por la estructura del registro. Solo 17 episodios tienen un único diagnóstico, que por definición de la cohorte es I50.

3.2 Esperanza, varianza y medidas de posición

EX   <- sum(fpX$X * fpX$px)                    # Σ x · p(x)
VX   <- sum((fpX$X - EX)^2 * fpX$px)            # Σ (x − E[X])² · p(x)  (divide por n)
EX_r <- mean(coh$X); VX_r <- var(coh$X)         # var() divide por n − 1
P_k  <- 1 - fpX[X <= k_X - 1, sum(px)]          # P(X ≥ k)
P_med <- fpX[X <= med_omega, sum(px)]           # P(X ≤ mediana de Ω)
mom <- data.table(
  Medida = c("E[X]", "Var[X]", "Mediana", "Moda", paste0("P(X ≥ ", k_X, ")"), paste0("P(X ≤ ", med_omega, ")")),
  `Desde la función de probabilidad` = c(paste0("Σ x · p(x) = ", fnum(EX, 6)),
                                         paste0("Σ (x − E[X])² · p(x) = ", fnum(VX, 6)),
                                         paste0("menor x con F(x) ≥ 0,5: ", mediana_X),
                                         paste0("x con mayor p(x): ", moda_X),
                                         paste0("1 − F(", k_X - 1, ") = ", fpct(P_k)),
                                         paste0("F(", med_omega, ") = ", fpct(P_med))),
  `Sobre los episodios (R)` = c(paste0("mean() = ", fnum(EX_r, 6)), paste0("var() = ", fnum(VX_r, 6)),
                                paste0("median() = ", median(coh$X)), "—",
                                paste0("mean(X >= ", k_X, ") = ", fpct(mean(coh$X >= k_X))),
                                paste0("mean(X <= ", med_omega, ") = ", fpct(mean(coh$X <= med_omega)))),
  Diferencia = c(fdif(EX - EX_r), fdif(VX - VX_r), fdif(mediana_X - median(coh$X), 0), "—",
                 fdif(P_k - mean(coh$X >= k_X)), fdif(P_med - mean(coh$X <= med_omega))))
tabla_html(mom, caption = paste0("Esperanza, varianza y medidas de posición de X: desde la función de probabilidad y verificadas sobre los ", fnum(nE), " episodios"),
           label = "momX")

Tabla 7. Esperanza, varianza y medidas de posición de X: desde la función de probabilidad y verificadas sobre los 52.312 episodios
Medida	Desde la función de probabilidad	Sobre los episodios (R)	Diferencia
E[X]	Σ x · p(x) = 12,423440	mean() = 12,423440	0,000000
Var[X]	Σ (x − E[X])² · p(x) = 32,692449	var() = 32,693074	-0,000625
Mediana	menor x con F(x) ≥ 0,5: 11	median() = 11	0
Moda	x con mayor p(x): 9	—	—
P(X ≥ 10)	1 − F(9) = 65,26%	mean(X >= 10) = 65,26%	0,000000
P(X ≤ 5)	F(5) = 5,83%	mean(X <= 5) = 5,83%	0,000000

Interpretación (Tabla 7): un episodio de la cohorte registra en promedio E[X] = 12,42 diagnósticos, con varianza 32,69 (desviación estándar 5,72). La esperanza calculada como Σ x · p(x) coincide con mean(). La varianza de var() es apenas mayor (diferencia 0,000625) porque divide por n − 1 = 52.311 en vez de n = 52.312: var() = Var[X] · n / (n − 1), y con n = 52.312 el factor n / (n − 1) es 1,000019, por lo que la diferencia es mínima. La media (12,42) es mayor que la mediana (11) y que la moda (9), lo que confirma la asimetría a la derecha.

Como k se usa el doble de la mediana de diagnósticos de todo Ω (5), es decir k = 10: 65,3% de los episodios de la cohorte registra 10 o más diagnósticos, y solo 5,8% registra 5 o menos, que es lo típico en Ω. Los episodios con insuficiencia cardíaca acumulan bastantes más diagnósticos que un egreso típico de la base.

3.3 Comparación según el rol diagnóstico: Eₚ contra Eₛ

rol_lbl <- c("Eₚ: I50 diagnóstico principal", "Eₛ: I50 solo diagnóstico secundario")
coh[, rol := factor(ifelse(Ep, rol_lbl[1], rol_lbl[2]), levels = rol_lbl)]
res_rol <- coh[, {
  t <- table(X); pxr <- as.numeric(t) / .N; xs <- as.integer(names(t))
  list(n = .N, EX = sum(xs * pxr), VX = sum((xs - sum(xs * pxr))^2 * pxr),
       mediana = xs[which(cumsum(pxr) >= 0.5)[1]], moda = xs[which.max(pxr)],
       minimo = min(xs), Pk = mean(X >= k_X), P35 = mean(X == 35))
}, by = rol][order(rol)]
tabla_html(res_rol[, .(`Rol diagnóstico` = rol, `Episodios` = fnum(n), `E[X]` = fnum(EX, 2),
                       `Var[X]` = fnum(VX, 2), Mediana = mediana, Moda = moda, `Mínimo` = minimo,
                       `P(X ≥ k)` = fpct(Pk, 1), `P(X = 35)` = fprob(P35))],
           caption = paste0("Resumen de la función de probabilidad de X por rol diagnóstico (k = ", k_X, "; cada fila se calcula sobre los episodios de su rol)"),
           label = "rolX")

Tabla 8. Resumen de la función de probabilidad de X por rol diagnóstico (k = 10; cada fila se calcula sobre los episodios de su rol)
Rol diagnóstico	Episodios	E[X]	Var[X]	Mediana	Moda	Mínimo	P(X ≥ k)	P(X = 35)
Eₚ: I50 diagnóstico principal	11.331	11,55	28,07	11	9	1	59,9%	0,20%
Eₛ: I50 solo diagnóstico secundario	40.981	12,66	33,70	12	9	2	66,7%	0,51%

fp_rol <- coh[, .(n = .N), by = .(rol, X)][, px := n / sum(n), by = rol]

fp_rol[, rol_fig := gsub("ₛ", "s", gsub("ₚ", "p", as.character(rol)))]   # la fuente del gráfico no tiene subíndices
ggplot(fp_rol, aes(x = X, y = px, color = rol_fig)) +
  geom_line(linewidth = 1) + geom_point(size = 1.6) +
  scale_color_manual(values = c(paleta[2], paleta[3]), name = NULL) +
  scale_x_continuous(breaks = seq(0, 35, 5)) +
  scale_y_continuous(labels = lab_pct) +
  labs(title = "Número de diagnósticos según el rol de I50 en el episodio",
       x = "x: diagnósticos registrados", y = "p(x) dentro del rol") +
  theme(legend.position = "bottom")

Figura 3. Función de probabilidad de X por separado para Eₚ y Eₛ. Cada curva suma 1 dentro de su rol.

Interpretación (Tabla 8 y Figura 3): los episodios con I50 solo como diagnóstico secundario (Eₛ, 40.981 episodios) registran en promedio 12,66 diagnósticos, 1,11 más que los episodios con I50 como diagnóstico principal (Eₚ, 11,55), y su mediana es 12 frente a 11. Parte de la diferencia es estructural: en Eₛ el primer campo tiene otro diagnóstico y I50 ocupa un campo secundario, por lo que X vale al menos 2, mientras que en Eₚ puede valer 1. La curva de Eₚ está desplazada a la izquierda y es algo más concentrada (varianza 28,1 frente a 33,7). La distribución de X en toda la cohorte es una mezcla de estas dos, dominada por Eₛ.

4 Bernoulli y binomial: fallecimientos en una auditoría de 20 episodios

RESULTADO 02

4.1 La variable Bernoulli Y

EY <- p; VY <- p * (1 - p)
bern <- data.table(
  Medida = c("p = P(Y = 1) = P(F | E)", "P(Y = 0) = 1 − p", "E[Y] = p", "Var[Y] = p(1 − p)"),
  `Modelo Bernoulli` = c(paste0(fnum(nEF), " / ", fnum(nE), " = ", fnum(p, 6)), fnum(1 - p, 6),
                         fnum(EY, 6), paste0(f4(p), " × ", f4(1 - p), " = ", fnum(VY, 6))),
  `Sobre los episodios (R)` = c(paste0("mean(Fal) = ", fnum(mean(coh$Fal), 6)), fnum(mean(!coh$Fal), 6),
                                paste0("mean(Fal) = ", fnum(mean(coh$Fal), 6)),
                                paste0("var(Fal) = ", fnum(var(coh$Fal), 6))),
  `Función de R` = codigo(c("dbinom(1, 1, p)", "dbinom(0, 1, p)", "—", "—")))
tabla_html(bern, caption = paste0("Variable Y ~ Bernoulli(p): parámetro, esperanza y varianza (n(E) = ", fnum(nE), ")"),
           label = "bern", escape = FALSE)

Tabla 9. Variable Y ~ Bernoulli(p): parámetro, esperanza y varianza (n(E) = 52.312)
Medida	Modelo Bernoulli	Sobre los episodios (R)	Función de R
p = P(Y = 1) = P(F \| E)	5.231 / 52.312 = 0,099996	mean(Fal) = 0,099996	`dbinom(1, 1, p)`
P(Y = 0) = 1 − p	0,900004	0,900004	`dbinom(0, 1, p)`
E[Y] = p	0,099996	mean(Fal) = 0,099996	`—`
Var[Y] = p(1 − p)	0,1000 × 0,9000 = 0,089997	var(Fal) = 0,089999	`—`

Interpretación (Tabla 9): un episodio elegido al azar de la cohorte termina en fallecimiento con probabilidad p = 10,00%, el mismo valor de P(F | E) de la Actividad 2. La esperanza de una Bernoulli es la propia proporción, y su varianza p(1 − p) = 0,0900 coincide con var() salvo el factor n / (n − 1) ya explicado para X.

4.2 S ~ Binomial(20, p) y S_L ~ Binomial(20, P(L | E))

S suma 20 variables Bernoulli: el número de fallecimientos entre 20 episodios elegidos al azar con reposición. Con reposición, cada elección es independiente y tiene la misma p. Se repite el cálculo para la estancia larga con P(L | E) = 44,70%.

kS <- 0:20
fpS <- data.table(k = kS, dS = dbinom(kS, 20, p), FS = pbinom(kS, 20, p),
                  dL = dbinom(kS, 20, pL), FL = pbinom(kS, 20, pL))
tabla_html(fpS[, .(k, `P(S = k)` = fprob(dS), `P(S ≤ k)` = fprob_cdf(FS),
                   `P(S_L = k)` = fprob(dL), `P(S_L ≤ k)` = fprob_cdf(FL))],
           caption = "Función de probabilidad (dbinom) y acumulada (pbinom) de S ~ Binomial(20, p) y de S_L ~ Binomial(20, P(L | E))",
           label = "fpS") |>
  add_header_above(c(" " = 1, "Fallecimientos: S" = 2, "Estancias largas: S_L" = 2)) |>
  scroll_box(height = "420px")

Tabla 10. Función de probabilidad (dbinom) y acumulada (pbinom) de S ~ Binomial(20, p) y de S_L ~ Binomial(20, P(L | E))
	Fallecimientos: S		Estancias largas: S_L
k	P(S = k)	P(S ≤ k)	P(S_L = k)	P(S_L ≤ k)
0	12,16%	12,16%	0,000716%	0,000716%
1	27,02%	39,18%	0,0116%	0,0123%
2	28,52%	67,69%	0,0889%	0,10%
3	19,01%	86,71%	0,43%	0,53%
4	8,98%	95,68%	1,48%	2,01%
5	3,19%	98,87%	3,83%	5,84%
6	0,89%	99,76%	7,74%	13,58%
7	0,20%	99,96%	12,51%	26,08%
8	0,0356%	99,9940%	16,42%	42,51%
9	0,00527%	99,99928%	17,70%	60,21%
10	0,000644%	99,999929%	15,74%	75,94%
11	0,0000650%	99,9999941%	11,56%	87,50%
12	0,00000542%	99,99999960%	7,01%	94,51%
13	0,000000371%	99,999999978%	3,49%	98,00%
14	0,0000000206%	99,99999999905%	1,41%	99,40%
15	0,000000000915%	99,999999999967%	0,46%	99,86%
16	< 0,0000000001%	100%	0,12%	99,97%
17	< 0,0000000001%	100%	0,0219%	99,9967%
18	< 0,0000000001%	100%	0,00295%	99,99973%
19	< 0,0000000001%	100%	0,000251%	99,999989%
20	< 0,0000000001%	100%	0,0000101%	100%

bin_res <- function(pp, nombre) {
  c(nombre = nombre, p = pp, P0 = dbinom(0, 20, pp), P0_mano = (1 - pp)^20,
    P2 = 1 - pbinom(1, 20, pp), q90 = qbinom(0.9, 20, pp), E = 20 * pp, V = 20 * pp * (1 - pp),
    asim = (1 - 2 * pp) / sqrt(20 * pp * (1 - pp)), moda = kS[which.max(dbinom(kS, 20, pp))])
}
rS <- bin_res(p, "S"); rL <- bin_res(pL, "S_L")
num <- function(r, x) as.numeric(r[[x]])
res_bin <- data.table(
  Medida = c("Parámetro p", "P(· = 0) con dbinom", "P(· = 0) a mano: (1 − p)^20",
             "P(· ≥ 2) = 1 − pbinom(1, 20, p)", "Percentil 90: qbinom(0,9, 20, p)",
             "P(· ≤ percentil 90)", "Moda", "E = 20p", "Var = 20p(1 − p)",
             "Asimetría (1 − 2p) / √(20p(1 − p))"),
  S = c(fpct(p), fprob(num(rS, "P0")), fprob(num(rS, "P0_mano")), fpct(num(rS, "P2")), num(rS, "q90"),
        fpct(pbinom(num(rS, "q90"), 20, p)), num(rS, "moda"), fnum(num(rS, "E"), 2), fnum(num(rS, "V"), 2),
        fnum(num(rS, "asim"), 3)),
  S_L = c(fpct(pL), fprob(num(rL, "P0")), fprob(num(rL, "P0_mano")), fpct(num(rL, "P2")), num(rL, "q90"),
          fpct(pbinom(num(rL, "q90"), 20, pL)), num(rL, "moda"), fnum(num(rL, "E"), 2), fnum(num(rL, "V"), 2),
          fnum(num(rL, "asim"), 3)))
setnames(res_bin, c("Medida", "S: fallecimientos entre 20", "S_L: estancias largas entre 20"))
tabla_html(res_bin, caption = "Probabilidades, percentil 90, esperanza y varianza de S y S_L", label = "resS")

Tabla 11. Probabilidades, percentil 90, esperanza y varianza de S y S_L
Medida	S: fallecimientos entre 20	S_L: estancias largas entre 20
Parámetro p	10,00%	44,70%
P(· = 0) con dbinom	12,16%	0,000716%
P(· = 0) a mano: (1 − p)^20	12,16%	0,000716%
P(· ≥ 2) = 1 − pbinom(1, 20, p)	60,82%	99,99%
Percentil 90: qbinom(0,9, 20, p)	4	12
P(· ≤ percentil 90)	95,68%	94,51%
Moda	2	9
E = 20p	2,00	8,94
Var = 20p(1 − p)	1,80	4,94
Asimetría (1 − 2p) / √(20p(1 − p))	0,596	0,048

fpS_l <- rbind(data.table(var = paste0("S: fallecimientos (p = ", fpct(p, 1), ")"), k = kS, d = fpS$dS, E = 20 * p, q90 = num(rS, "q90")),
               data.table(var = paste0("S_L: estancias largas (p = ", fpct(pL, 1), ")"), k = kS, d = fpS$dL, E = 20 * pL, q90 = num(rL, "q90")))
ggplot(fpS_l, aes(x = k, y = d)) +
  geom_col(fill = paleta[1], width = 0.75) +
  geom_vline(aes(xintercept = E), linetype = "dashed", color = paleta[2], linewidth = 0.9) +
  geom_vline(aes(xintercept = q90), linetype = "dotted", color = paleta[6], linewidth = 0.9) +
  facet_wrap(~ var, ncol = 2) +
  scale_x_continuous(breaks = seq(0, 20, 2)) +
  scale_y_continuous(labels = lab_pct) +
  labs(title = "Número de episodios con el desenlace entre 20 auditados al azar",
       x = "k (episodios entre 20)", y = "P(· = k)")

Figura 4. Función de probabilidad de S ~ Binomial(20; p) y de S_L ~ Binomial(20; P(L | E)). La línea discontinua marca la esperanza y la punteada el percentil 90.

Interpretación (Tabla 10, Tabla 11 y Figura 4): en una auditoría de 20 episodios de la cohorte elegidos al azar, lo esperable son E[S] = 2,00 fallecimientos (Var[S] = 1,80). Hay 12,2% de probabilidad de no encontrar ninguno, P(S = 0) = (1 − 0,1000)^20, y 60,8% de encontrar dos o más. El percentil 90 es 4: qbinom entrega el menor k con P(S ≤ k) ≥ 0,90, y aquí P(S ≤ 4) = 95,7%. Por eso, encontrar más de 4 fallecimientos en 20 episodios ocurriría en menos del 10% de las auditorías si todos los episodios compartieran p.

Las formas son distintas. Como p = 10,0% es pequeña, S se concentra cerca de 0 y tiene cola a la derecha (asimetría 0,60). S_L, con P(L | E) = 44,7%, está centrada en 8,9, es casi simétrica (asimetría 0,05) y es más dispersa (varianza 4,94 frente a 1,80), porque p(1 − p) es máxima cerca de p = 0,5. En S_L, P(S_L = 0) = 0,000716%: es casi imposible que una auditoría de 20 episodios no encuentre ninguna estancia larga.

4.3 Fallecimientos observados y esperados por establecimiento

Para cada establecimiento con al menos 20 episodios de la cohorte se compara el número observado de fallecimientos con el esperado n · p, y se calcula la probabilidad binomial de observar un valor igual o mayor, P(S_h ≥ observado) = 1 - pbinom(obs - 1, n, p), con S_h ~ Binomial(n_h, p). Se agrega la probabilidad de observar un valor igual o menor, para describir también los establecimientos con menos fallecimientos que lo esperado.

Probabilidad de cola no es una prueba Estas probabilidades indican qué tan inusual sería cada conteo si todos los episodios compartieran la misma p. Se usan solo para ordenar y describir establecimientos. No son una prueba de hipótesis ni evidencia sobre la calidad de la atención, porque los establecimientos atienden casos distintos (severidad, edad, rol de I50, derivaciones).

est20 <- est[N_h >= 20]
est20[, esperado := N_h * p]
est20[, cola_sup := pbinom(K_F - 1, N_h, p, lower.tail = FALSE)]   # P(S_h ≥ observado)
est20[, cola_inf := pbinom(K_F, N_h, p)]                            # P(S_h ≤ observado)
setorder(est20, cola_sup)
tabla_html(est20[, .(Establecimiento = nombre_hospital, `Episodios n` = fnum(N_h),
                     `Fallecimientos observados` = fnum(K_F), `Esperados n · p` = fnum(esperado, 1),
                     `Observados / esperados` = fnum(K_F / esperado, 2),
                     `P(S_h ≥ observado)` = fprob_cdf(cola_sup), `P(S_h ≤ observado)` = fprob_cdf(cola_inf))],
           caption = paste0("Fallecimientos observados y esperados en los ", nrow(est20), " establecimientos con al menos 20 episodios de la cohorte (p = ", fpct(p), "), ordenados por P(S_h ≥ observado)"),
           label = "est") |>
  scroll_box(height = "420px")

Tabla 12. Fallecimientos observados y esperados en los 72 establecimientos con al menos 20 episodios de la cohorte (p = 10,00%), ordenados por P(S_h ≥ observado)
Establecimiento	Episodios n	Fallecimientos observados	Esperados n · p	Observados / esperados	P(S_h ≥ observado)	P(S_h ≤ observado)
Complejo Hospitalario San José (Santiago, Independencia)	1.275	175	127,5	1,37	0,00138%	99,99905%
Hospital Clínico Metropolitano La Florida Dra. Eloisa Díaz Inzunza	962	135	96,2	1,40	0,00440%	99,9970%
Hospital Dr. Leonardo Guzmán (Antofagasta)	600	90	60,0	1,50	0,00759%	99,9953%
Hospital Presidente Carlos Ibáñez del Campo (Linares)	643	92	64,3	1,43	0,0335%	99,98%
Hospital de Castro	339	51	33,9	1,50	0,22%	99,87%
Hospital Dr. Félix Bulnes Cerda (Santiago, Quinta Normal)	1.492	181	149,2	1,21	0,42%	99,67%
Hospital Dr. Antonio Tirado Lanas (Ovalle)	189	29	18,9	1,53	1,35%	99,24%
Hospital San José (Melipilla)	356	48	35,6	1,35	2,10%	98,59%
Hospital El Pino (Santiago, San Bernardo)	1.118	133	111,8	1,19	2,15%	98,30%
Hospital Dr. Luis Tisné B. (Santiago, Peñalolén)	1.014	120	101,4	1,18	3,13%	97,51%
Hospital San José (Victoria)	230	32	23,0	1,39	3,54%	97,76%
Hospital Provincial del Huasco Monseñor Fernando Ariztía Ruiz (Vallenar)	319	41	31,9	1,29	5,80%	95,95%
Hospital Padre Alberto Hurtado (San Ramón)	943	109	94,3	1,16	6,38%	94,82%
Hospital Dr. César Garavagno Burotto (Talca)	1.452	163	145,2	1,12	6,68%	94,34%
Hospital Dr. Mauricio Heyermann (Angol)	458	56	45,8	1,22	6,84%	94,89%
Hospital Claudio Vicuña (San Antonio)	598	71	59,8	1,19	7,49%	94,18%
Hospital Clínico Metropolitano El Carmen Doctor Luis Valentín Ferrada	1.132	126	113,2	1,11	11,25%	90,48%
Hospital San Luis (Buin)	393	46	39,3	1,17	14,88%	88,53%
Hospital Dr. Juan Noé Crevanni (Arica)	479	55	47,9	1,15	15,72%	87,53%
Hospital Dr. Hernán Henríquez Aravena (Temuco)	1.020	112	102,0	1,10	16,05%	86,30%
Hospital Dr. Carlos Cisternas (Calama)	259	31	25,9	1,20	16,94%	87,55%
Hospital Dr. Gustavo Fricke (Viña del Mar)	1.198	130	119,8	1,09	17,45%	84,86%
Hospital San Juan de Dios (San Fernando)	279	33	27,9	1,18	17,80%	86,71%
Hospital Provincial Dr. Rafael Avaría (Curanilahue)	439	50	43,9	1,14	18,51%	85,29%
Hospital de Pitrufquén	311	36	31,1	1,16	20,04%	84,63%
Hospital San Juan de Dios (Cauquenes)	407	46	40,7	1,13	21,14%	83,17%
Hospital San Juan de Dios (Los Andes)	390	44	39,0	1,13	22,08%	82,44%
Hospital Dr. Eduardo Pereira Ramírez (Valparaíso)	116	14	11,6	1,21	26,92%	81,79%
Hospital Carlos Van Buren (Valparaíso)	1.287	134	128,7	1,04	32,41%	70,82%
Hospital de San Camilo (San Felipe)	565	60	56,5	1,06	33,13%	71,70%
Hospital de Quilpué	296	32	29,6	1,08	34,81%	71,89%
Hospital Adalberto Steeger (Talagante)	438	46	43,8	1,05	38,57%	67,27%
Hospital San José del Carmen (Copiapó)	391	40	39,1	1,02	46,41%	60,16%
Hospital San Juan de Dios (Curicó)	1.070	108	107,0	1,01	47,41%	56,61%
Hospital de Niños Dr. Luis Calvo Mackenna (Santiago, Providencia)	36	4	3,6	1,11	49,15%	71,08%
Hospital San Pablo (Coquimbo)	694	69	69,4	0,99	53,86%	51,19%
Hospital de Santa Cruz	192	19	19,2	0,99	55,43%	54,16%
Hospital San Juan de Dios (Santiago, Santiago)	1.006	99	100,6	0,98	58,20%	45,96%
Hospital de Urgencia Asistencia Pública Dr. Alejandro del Río	1.237	121	123,7	0,98	61,46%	42,23%
Hospital Intercultural de Nueva Imperial	396	38	39,6	0,96	62,99%	43,56%
Hospital Dr. Abraham Godoy (Lautaro)	255	24	25,5	0,94	65,32%	42,80%
Cód. 200717	514	49	51,4	0,95	65,90%	39,71%
Hospital Clínico de Magallanes Dr. Lautaro Navarro Avaria	609	58	60,9	0,95	67,17%	37,92%
Hospital Del Salvador (Santiago, Providencia)	954	91	95,4	0,95	69,78%	34,14%
Hospital Dr. Exequiel González Cortés (Santiago, San Miguel)	26	2	2,6	0,77	74,87%	51,05%
Hospital San Martín (Quillota)	595	55	59,5	0,92	74,96%	29,69%
Hospital San José (Parral)	398	36	39,8	0,90	76,03%	29,62%
Hospital Las Higueras (Talcahuano)	1.386	131	138,6	0,95	76,40%	26,50%
Hospital de San Carlos	695	64	69,5	0,92	77,37%	26,71%
Hospital Barros Luco Trudeau (Santiago, San Miguel)	1.987	189	198,7	0,95	77,58%	24,76%
Hospital San José (Coronel)	466	42	46,6	0,90	78,22%	26,75%
Hospital Clínico de Niños Dr. Roberto del Río (Santiago, Independencia)	53	4	5,3	0,75	78,94%	37,83%
Hospital Regional de Rancagua	1.611	150	161,1	0,93	83,20%	19,01%
Instituto de Neurocirugía Dr. Alfonso Asenjo	34	2	3,4	0,59	86,71%	32,55%
Hospital Dr. Ernesto Torres Galdames (Iquique)	561	48	56,1	0,86	88,88%	14,17%
Hospital de Villarrica	247	19	24,7	0,77	90,98%	13,33%
Hospital de Lota	406	33	40,6	0,81	91,31%	11,82%
Hospital San Juan de Dios (La Serena)	824	68	82,4	0,83	96,11%	5,04%
Hospital de Tomé	375	28	37,5	0,75	96,20%	5,63%
Hospital de Constitución	369	27	36,9	0,73	96,90%	4,69%
Instituto Traumatológico Dr. Teodoro Gebauer	52	2	5,2	0,38	97,17%	9,66%
Hospital Clínico Regional Dr. Guillermo Grant Benavente (Concepción)	2.078	180	207,8	0,87	98,23%	2,13%
Hospital Clínico San Borja-Arriarán (Santiago, Santiago)	458	33	45,8	0,72	98,43%	2,38%
Hospital Penco - Lirquén	293	19	29,3	0,65	98,65%	2,32%
Hospital Base San José de Osorno	1.381	113	138,1	0,82	99,07%	1,20%
Hospital Regional (Coihaique)	306	19	30,6	0,62	99,27%	1,31%
Hospital Clínico Herminda Martín (Chillán)	1.382	112	138,2	0,81	99,30%	0,91%
Hospital Clínico Regional (Valdivia)	1.042	80	104,2	0,77	99,58%	0,58%
Complejo Hospitalario Dr. Sótero del Río (Santiago, Puente Alto)	3.138	269	313,8	0,86	99,70%	0,36%
Hospital de Puerto Montt	1.475	112	147,5	0,76	99,94%	0,0832%
Complejo Asistencial Dr. Víctor Ríos Ruiz (Los Ángeles)	1.606	117	160,6	0,73	99,9934%	0,00940%
Instituto Nacional de Enfermedades Respiratorias y Cirugía Torácica	717	36	71,7	0,50	99,999958%	0,0000892%

# cuántos establecimientos tendrían una cola ≤ 5% si todos compartieran p (cálculo exacto, binomial discreta)
c_sup <- qbinom(0.95, est20$N_h, p) + 1                 # menor c con P(S ≥ c) ≤ 0,05
esp_sup <- sum(pbinom(c_sup - 1, est20$N_h, p, lower.tail = FALSE))
c_inf <- qbinom(0.05, est20$N_h, p)                     # P(S ≤ c) ≥ 0,05; se baja hasta ≤ 0,05
c_inf <- ifelse(pbinom(c_inf, est20$N_h, p) <= 0.05, c_inf, c_inf - 1)
esp_inf <- sum(ifelse(c_inf >= 0, pbinom(c_inf, est20$N_h, p), 0))
obs_sup <- sum(est20$cola_sup <= 0.05); obs_inf <- sum(est20$cola_inf <= 0.05)
# razón de dispersión: promedio de (obs − esp)² / (n p (1 − p)); vale cerca de 1 si todos comparten p
razon_disp <- mean((est20$K_F - est20$esperado)^2 / (est20$N_h * p * (1 - p)))
rango_p_h <- range(est20$K_F / est20$N_h)

banda <- data.table(n = seq(20, max(est20$N_h), by = 1))
banda[, lo := qbinom(0.05, n, p) / n][, hi := qbinom(0.95, n, p) / n]
est20[, grupo := fifelse(cola_sup <= 0.05, "Más fallecimientos que lo esperado (cola ≤ 5%)",
                 fifelse(cola_inf <= 0.05, "Menos fallecimientos que lo esperado (cola ≤ 5%)", "Dentro de lo habitual"))]
ggplot(est20, aes(x = N_h, y = K_F / N_h)) +
  geom_ribbon(data = banda, aes(x = n, ymin = lo, ymax = hi), inherit.aes = FALSE, fill = paleta[4], alpha = 0.35) +
  geom_hline(yintercept = p, linetype = "dashed", color = paleta[3]) +
  geom_point(aes(color = grupo), size = 2.6, alpha = 0.9) +
  scale_color_manual(values = c("Más fallecimientos que lo esperado (cola ≤ 5%)" = paleta[2],
                                "Menos fallecimientos que lo esperado (cola ≤ 5%)" = paleta[6],
                                "Dentro de lo habitual" = paleta[1]), name = NULL) +
  scale_x_log10(labels = function(x) fnum(x)) +
  scale_y_continuous(labels = lab_pct) +
  labs(title = "Fallecimientos por establecimiento frente a la binomial con p común",
       subtitle = paste0("p = ", fpct(p, 1), " (línea discontinua); eje x en escala logarítmica"),
       x = "Episodios de la cohorte en el establecimiento (n)", y = "Proporción de fallecimientos") +
  theme(legend.position = "bottom") + guides(color = guide_legend(ncol = 1))

Figura 5. Proporción de fallecimientos por establecimiento según su número de episodios. La banda contiene el 90% central de la binomial con la p común de la cohorte; los puntos naranjos tienen P(S_h ≥ obs) ≤ 5% y los oscuros P(S_h ≤ obs) ≤ 5%.

Interpretación (Tabla 12 y Figura 5): los 72 establecimientos con episodios de la cohorte tienen todos al menos 20, y su proporción de fallecimientos va de 3,8% a 15,3%. Si todos los episodios compartieran p = 10,0%, se esperarían unos 3,1 establecimientos con P(S_h ≥ observado) ≤ 5% y unos 3,0 con P(S_h ≤ observado) ≤ 5%. Se observan 11 y 12, respectivamente. Además, la razón de dispersión (promedio de (observado − esperado)² / (n p (1 − p))) es 3,29, cuando bajo una binomial común sería cercana a 1. Los conteos por establecimiento varían bastante más de lo que permite una p única. La banda de Figura 5 se angosta a medida que crece n: en un establecimiento grande, una diferencia de pocos puntos en la proporción ya queda fuera de lo habitual para una p común, mientras que en uno pequeño la binomial admite proporciones muy distintas. Este orden describe la heterogeneidad de los establecimientos y no permite calificar su atención.

4.4 Supuestos de la binomial en estos datos

Ensayos independientes. En el experimento definido (20 episodios de toda la cohorte, con reposición) la independencia se cumple por construcción. En los datos no es tan claro: la cohorte tiene 52.312 episodios de 39.995 beneficiarios distintos, y 39,5% de los episodios pertenece a un beneficiario con dos o más episodios en la cohorte (hasta 20). Esos episodios no son independientes: a lo más uno de ellos puede terminar en fallecimiento. Una auditoría real dentro de un establecimiento o de un período también agrupa episodios parecidos.
Misma probabilidad p en todos los ensayos. No se cumple dentro de la cohorte. En la Actividad 2, P(F | Eₛ) fue mayor que P(F | Eₚ), y la proporción de fallecimientos cambia entre niveles de severidad. Entre establecimientos va de 3,8% a 15,3% (Tabla 12), con una razón de dispersión de 3,29. La binomial con una p común describe bien una auditoría de episodios elegidos al azar de toda la cohorte (se verifica por remuestreo en el Resultado 06), pero no describe los conteos de cada establecimiento, que agrupan episodios con distinta p.

5 Hipergeométrica: auditar sin reposición en un establecimiento

RESULTADO 03

n_aud <- if (any(est$N_h >= 10)) 10 else 5
setorder(est, -N_h, COD_HOSPITAL)
h_g <- est[1]                                       # mayor número de episodios
h_c <- est[N_h >= n_aud][order(N_h, COD_HOSPITAL)][1]   # menor número entre los que tienen al menos 10

Un auditor revisa 10 episodios sin reposición dentro de un establecimiento. Se eligen el establecimiento con más episodios de la cohorte, Complejo Hospitalario Dr. Sótero del Río (Santiago, Puente Alto), y el de menos episodios entre los que tienen al menos 10, Hospital Dr. Exequiel González Cortés (Santiago, San Miguel). H cuenta los fallecimientos entre los 10 episodios auditados. Como el establecimiento pequeño tiene solo 2 fallecimientos, se repite el cálculo con K = estancias largas, que dan más valores posibles.

nk <- data.table(
  Establecimiento = c(h_g$nombre_hospital, h_c$nombre_hospital),
  `Rol en la comparación` = c("Mayor número de episodios", paste0("Menor número entre los que tienen al menos ", n_aud)),
  `N (episodios de E)` = fnum(c(h_g$N_h, h_c$N_h)),
  `K (fallecimientos)` = fnum(c(h_g$K_F, h_c$K_F)),
  `K / N` = fpct(c(h_g$K_F / h_g$N_h, h_c$K_F / h_c$N_h)),
  `K (estancias largas)` = fnum(c(h_g$K_L, h_c$K_L)),
  `Fracción muestral n / N` = fpct(n_aud / c(h_g$N_h, h_c$N_h)))
tabla_html(nk, caption = paste0("Valores de N y K en los dos establecimientos elegidos (n = ", n_aud, " episodios auditados)"), label = "NK")

Tabla 13. Valores de N y K en los dos establecimientos elegidos (n = 10 episodios auditados)
Establecimiento	Rol en la comparación	N (episodios de E)	K (fallecimientos)	K / N	K (estancias largas)	Fracción muestral n / N
Complejo Hospitalario Dr. Sótero del Río (Santiago, Puente Alto)	Mayor número de episodios	3.138	269	8,57%	1.041	0,32%
Hospital Dr. Exequiel González Cortés (Santiago, San Miguel)	Menor número entre los que tienen al menos 10	26	2	7,69%	21	38,46%

hiper_tab <- function(h, K) {
  k <- 0:n_aud; Nh <- h$N_h
  data.table(k = k, hip = dhyper(k, K, Nh - K, n_aud), bin = dbinom(k, n_aud, K / Nh))
}
hg_F <- hiper_tab(h_g, h_g$K_F); hc_F <- hiper_tab(h_c, h_c$K_F)
hg_L <- hiper_tab(h_g, h_g$K_L); hc_L <- hiper_tab(h_c, h_c$K_L)
tab_h <- data.table(k = 0:n_aud,
                    g_h = fprob(hg_F$hip), g_b = fprob(hg_F$bin), g_d = fpp_dif(hg_F$hip - hg_F$bin),
                    c_h = fprob(hc_F$hip), c_b = fprob(hc_F$bin), c_d = fpp_dif(hc_F$hip - hc_F$bin))
setnames(tab_h, c("k", "Hipergeométrica (dhyper)", "Binomial p = K/N (dbinom)", "Diferencia",
                  "Hipergeométrica (dhyper)", "Binomial p = K/N (dbinom)", "Diferencia"))
cab <- c(1, 3, 3); names(cab) <- c(" ", paste0(h_g$nombre_hospital, " (N = ", fnum(h_g$N_h), ", K = ", h_g$K_F, ")"),
                                   paste0(h_c$nombre_hospital, " (N = ", fnum(h_c$N_h), ", K = ", h_c$K_F, ")"))
tabla_html(tab_h, caption = paste0("Función de probabilidad de H (fallecimientos entre ", n_aud, " episodios auditados sin reposición) y de la binomial con p = K / N"),
           label = "hip") |> add_header_above(cab)

Tabla 14. Función de probabilidad de H (fallecimientos entre 10 episodios auditados sin reposición) y de la binomial con p = K / N
	Complejo Hospitalario Dr. Sótero del Río (Santiago, Puente Alto) (N = 3.138, K = 269)			Hospital Dr. Exequiel González Cortés (Santiago, San Miguel) (N = 26, K = 2)
k	Hipergeométrica (dhyper)	Binomial p = K/N (dbinom)	Diferencia	Hipergeométrica (dhyper)	Binomial p = K/N (dbinom)	Diferencia
0	40,76%	40,81%	-0,055 p.p.	36,92%	44,91%	-7,991 p.p.
1	38,33%	38,26%	0,069 p.p.	49,23%	37,43%	11,803 p.p.
2	16,16%	16,14%	0,014 p.p.	13,85%	14,04%	-0,189 p.p.
3	4,02%	4,04%	-0,017 p.p.	0%	3,12%	-3,119 p.p.
4	0,65%	0,66%	-0,009 p.p.	0%	0,45%	-0,455 p.p.
5	0,0726%	0,0745%	-0,002 p.p.	0%	0,0455%	-0,045 p.p.
6	0,00557%	0,00582%	0,000 p.p.	0%	0,00316%	-0,003 p.p.
7	0,000292%	0,000312%	0,000 p.p.	0%	0,000150%	0,000 p.p.
8	0,0000100%	0,0000110%	0,000 p.p.	0%	0,00000470%	0,000 p.p.
9	0,000000203%	0,000000229%	0,000 p.p.	0%	0,0000000870%	0,000 p.p.
10	0,00000000184%	0,00000000214%	0,000 p.p.	0%	0,000000000725%	0,000 p.p.

envolver <- function(s, w = 38) paste(strwrap(s, width = w), collapse = "\n")
lbl_g <- paste0(envolver(h_g$nombre_hospital), "\n(n / N = ", fpct(n_aud / h_g$N_h, 1), ")")
lbl_c <- paste0(envolver(h_c$nombre_hospital), "\n(n / N = ", fpct(n_aud / h_c$N_h, 1), ")")
fig_h <- rbind(
  melt(hg_F[, est := lbl_g][, ev := "Fallecimientos"], id.vars = c("k", "est", "ev")),
  melt(hc_F[, est := lbl_c][, ev := "Fallecimientos"], id.vars = c("k", "est", "ev")),
  melt(hg_L[, est := lbl_g][, ev := "Estancias largas"], id.vars = c("k", "est", "ev")),
  melt(hc_L[, est := lbl_c][, ev := "Estancias largas"], id.vars = c("k", "est", "ev")))
fig_h[, modelo := factor(ifelse(variable == "hip", "Hipergeométrica (sin reposición)", "Binomial (con reposición)"),
                         levels = c("Hipergeométrica (sin reposición)", "Binomial (con reposición)"))]
fig_h[, est := factor(est, levels = c(lbl_g, lbl_c))]
fig_h[, ev := factor(ev, levels = c("Fallecimientos", "Estancias largas"))]
ggplot(fig_h, aes(x = k, y = value, fill = modelo)) +
  geom_col(position = position_dodge(width = 0.8), width = 0.75) +
  facet_grid(ev ~ est) +
  scale_fill_manual(values = c(paleta[3], paleta[5]), name = NULL) +
  scale_x_continuous(breaks = 0:n_aud) +
  scale_y_continuous(labels = lab_pct) +
  labs(title = "Auditoría de 10 episodios: sin reposición contra con reposición",
       x = "k (episodios con el desenlace entre los auditados)", y = "P(· = k)") +
  theme(legend.position = "bottom")

Figura 6. Hipergeométrica (sin reposición) y binomial (con reposición) para los fallecimientos y las estancias largas entre 10 episodios auditados en cada establecimiento.

var_hip <- function(N, K, n) n * (K / N) * (1 - K / N) * (N - n) / (N - 1)
var_bin <- function(N, K, n) n * (K / N) * (1 - K / N)
fila_h <- function(h, K, ev, tabla) {
  data.table(Establecimiento = h$nombre_hospital, Desenlace = ev, N = fnum(h$N_h), K = fnum(K),
             `n / N` = fpct(n_aud / h$N_h, 1),
             `Mayor |diferencia|` = paste0(fnum(max(abs(tabla$hip - tabla$bin)) * 100, 2), " p.p. (k = ",
                                           tabla$k[which.max(abs(tabla$hip - tabla$bin))], ")"),
             `Var hipergeométrica` = fnum(var_hip(h$N_h, K, n_aud), 4),
             `Var binomial` = fnum(var_bin(h$N_h, K, n_aud), 4),
             `Factor (N − n) / (N − 1)` = fnum((h$N_h - n_aud) / (h$N_h - 1), 4),
             dif_num = max(abs(tabla$hip - tabla$bin)))
}
res_h <- rbind(fila_h(h_g, h_g$K_F, "Fallecimientos", hg_F), fila_h(h_c, h_c$K_F, "Fallecimientos", hc_F),
               fila_h(h_g, h_g$K_L, "Estancias largas", hg_L), fila_h(h_c, h_c$K_L, "Estancias largas", hc_L))
tabla_html(res_h[, !"dif_num"], caption = "Mayor diferencia entre la hipergeométrica y la binomial, y varianza de cada modelo",
           label = "hip_res")

Tabla 15. Mayor diferencia entre la hipergeométrica y la binomial, y varianza de cada modelo
Establecimiento	Desenlace	N	K	n / N	Mayor \|diferencia\|	Var hipergeométrica	Var binomial	Factor (N − n) / (N − 1)
Complejo Hospitalario Dr. Sótero del Río (Santiago, Puente Alto)	Fallecimientos	3.138	269	0,3%	0,07 p.p. (k = 1)	0,7815	0,7837	0,9971
Hospital Dr. Exequiel González Cortés (Santiago, San Miguel)	Fallecimientos	26	2	38,5%	11,80 p.p. (k = 1)	0,4544	0,7101	0,6400
Complejo Hospitalario Dr. Sótero del Río (Santiago, Puente Alto)	Estancias largas	3.138	1.041	0,3%	0,04 p.p. (k = 3)	2,2105	2,2169	0,9971
Hospital Dr. Exequiel González Cortés (Santiago, San Miguel)	Estancias largas	26	21	38,5%	8,17 p.p. (k = 8)	0,9941	1,5533	0,6400

P0_hc <- hc_F$hip[1]; P0_hc_b <- hc_F$bin[1]

Interpretación (Tabla 13, Tabla 14, Tabla 15 y Figura 6): en el establecimiento de mayor volumen, Complejo Hospitalario Dr. Sótero del Río (Santiago, Puente Alto), la auditoría revisa 10 de 3.138 episodios, una fracción muestral de 0,32%. Sacar un episodio casi no cambia la composición de los que quedan, así que revisar sin reposición equivale en la práctica a revisar con reposición: la mayor diferencia entre la hipergeométrica y la binomial es 0,069 puntos porcentuales y las varianzas casi coinciden (factor 0,9971).

En el establecimiento pequeño, Hospital Dr. Exequiel González Cortés (Santiago, San Miguel), la auditoría revisa 10 de 26 episodios (n / N = 38,5%), y entre ellos hay solo K = 2 fallecimientos. Cada episodio auditado cambia la composición de los que quedan, y H no puede superar 2. P(H = 0) = 36,9%, mientras que la binomial con p = 2 / 26 da 44,9%, y la mayor diferencia llega a 11,8 puntos porcentuales. La varianza hipergeométrica es la binomial multiplicada por (N − n) / (N − 1) = 0,640: sin reposición la auditoría es más precisa porque cada episodio revisado deja de estar disponible. Con estancias largas (K = 21 de 26) ocurre lo mismo: la mayor diferencia es 8,2 puntos porcentuales y la varianza hipergeométrica (0,9941) es menor que la binomial (1,5533), mientras que en el establecimiento grande la diferencia es 0,04 puntos. Las dos distribuciones se parecen más en el establecimiento grande porque su fracción muestral es casi cero. Usar la binomial para auditar un establecimiento pequeño sin reposición exageraría la variabilidad.

6 Geométrica y binomial negativa: revisar fichas hasta encontrar casos

RESULTADO 04

Un auditor revisa, una por una, fichas de episodios marcados por el marcador M de la Actividad 2: especialidad médica Medicina interna, Cardiología o Geriatría y edad al ingreso de 65 años o más. La probabilidad de que una ficha marcada tenga I50 es el VPP = P(E | M) = 26,60%. Si las fichas se eligen al azar (con reposición), el número de fichas revisadas hasta encontrar el primer caso es T ~ Geométrica(VPP). Sin marcador, revisando episodios al azar de Ω, la probabilidad de éxito es P(E) = 5,99% y el número de fichas es T₀ ~ Geométrica(P(E)). La convención es la declarada en la sección 2.5: T cuenta ensayos, y P(T = k) = dgeom(k - 1, p).

6.1 Función de probabilidad de T y T₀ para k = 1 a 50

kT <- 1:50
fpT <- data.table(k = kT, dT = dgeom(kT - 1, vpp), FT = pgeom(kT - 1, vpp),
                  dT0 = dgeom(kT - 1, pE), FT0 = pgeom(kT - 1, pE))
tabla_html(fpT[, .(k, `P(T = k)` = fprob(dT), `P(T ≤ k)` = fprob_cdf(FT),
                   `P(T₀ = k)` = fprob(dT0), `P(T₀ ≤ k)` = fprob_cdf(FT0))],
           caption = "Función de probabilidad (dgeom(k − 1, p)) y acumulada (pgeom(k − 1, p)) de T (con marcador) y T₀ (sin marcador)",
           label = "fpT") |>
  add_header_above(c(" " = 1, "T: con marcador, p = VPP" = 2, "T₀: sin marcador, p = P(E)" = 2)) |>
  scroll_box(height = "420px")

Tabla 16. Función de probabilidad (dgeom(k − 1, p)) y acumulada (pgeom(k − 1, p)) de T (con marcador) y T₀ (sin marcador)
	T: con marcador, p = VPP		T₀: sin marcador, p = P(E)
k	P(T = k)	P(T ≤ k)	P(T₀ = k)	P(T₀ ≤ k)
1	26,60%	26,60%	5,99%	5,99%
2	19,52%	46,12%	5,63%	11,62%
3	14,33%	60,45%	5,29%	16,91%
4	10,52%	70,97%	4,97%	21,88%
5	7,72%	78,69%	4,68%	26,56%
6	5,67%	84,36%	4,40%	30,96%
7	4,16%	88,52%	4,13%	35,09%
8	3,05%	91,57%	3,89%	38,98%
9	2,24%	93,81%	3,65%	42,63%
10	1,65%	95,46%	3,43%	46,07%
11	1,21%	96,67%	3,23%	49,29%
12	0,89%	97,55%	3,04%	52,33%
13	0,65%	98,20%	2,85%	55,18%
14	0,48%	98,68%	2,68%	57,87%
15	0,35%	99,03%	2,52%	60,39%
16	0,26%	99,29%	2,37%	62,76%
17	0,19%	99,48%	2,23%	64,99%
18	0,14%	99,62%	2,10%	67,09%
19	0,10%	99,72%	1,97%	69,06%
20	0,0747%	99,79%	1,85%	70,91%
21	0,0549%	99,85%	1,74%	72,65%
22	0,0403%	99,89%	1,64%	74,29%
23	0,0296%	99,92%	1,54%	75,83%
24	0,0217%	99,94%	1,45%	77,28%
25	0,0159%	99,96%	1,36%	78,64%
26	0,0117%	99,97%	1,28%	79,92%
27	0,00858%	99,98%	1,20%	81,12%
28	0,00630%	99,98%	1,13%	82,25%
29	0,00462%	99,99%	1,06%	83,31%
30	0,00339%	99,9906%	1,00%	84,31%
31	0,00249%	99,9931%	0,94%	85,25%
32	0,00183%	99,9949%	0,88%	86,13%
33	0,00134%	99,9962%	0,83%	86,96%
34	0,000986%	99,9972%	0,78%	87,74%
35	0,000723%	99,9980%	0,73%	88,48%
36	0,000531%	99,9985%	0,69%	89,17%
37	0,000390%	99,9989%	0,65%	89,82%
38	0,000286%	99,99921%	0,61%	90,43%
39	0,000210%	99,99942%	0,57%	91,00%
40	0,000154%	99,99957%	0,54%	91,54%
41	0,000113%	99,99968%	0,51%	92,04%
42	0,0000831%	99,99977%	0,48%	92,52%
43	0,0000610%	99,99983%	0,45%	92,97%
44	0,0000448%	99,99987%	0,42%	93,39%
45	0,0000329%	99,999909%	0,40%	93,79%
46	0,0000241%	99,999933%	0,37%	94,16%
47	0,0000177%	99,999951%	0,35%	94,51%
48	0,0000130%	99,999964%	0,33%	94,84%
49	0,00000954%	99,999973%	0,31%	95,15%
50	0,00000700%	99,999980%	0,29%	95,44%

fig_T <- rbind(data.table(v = paste0("T: con marcador (p = ", fpct(vpp, 1), ")"), k = kT, d = fpT$dT, E = 1 / vpp),
               data.table(v = paste0("T0: sin marcador (p = ", fpct(pE, 1), ")"), k = kT, d = fpT$dT0, E = 1 / pE))
ggplot(fig_T, aes(x = k, y = d, fill = v)) +
  geom_col(width = 0.8, show.legend = FALSE) +
  geom_vline(aes(xintercept = E), linetype = "dashed", color = paleta[6]) +
  facet_wrap(~ v, ncol = 1, scales = "free_y") +
  scale_fill_manual(values = c(paleta[1], paleta[2])) +
  scale_x_continuous(breaks = c(1, seq(5, 50, 5))) +
  scale_y_continuous(labels = lab_pct) +
  labs(title = "Fichas revisadas hasta encontrar el primer episodio con insuficiencia cardíaca",
       x = "k (fichas revisadas)", y = "P(· = k)")

Figura 7. Función de probabilidad de T (con marcador, p = VPP) y T₀ (sin marcador, p = P(E)) para k = 1 a 50. Las líneas discontinuas marcan la esperanza de cada una.

geo_res <- function(pp) c(E = 1 / pp, V = (1 - pp) / pp^2, med = qgeom(0.5, pp) + 1,
                         P50 = pgeom(49, pp, lower.tail = FALSE), P50_mano = (1 - pp)^50,
                         P10 = pgeom(9, pp))
gT <- geo_res(vpp); gT0 <- geo_res(pE)
ahorro <- gT0[["E"]] - gT[["E"]]
res_geo <- data.table(
  Medida = c("Probabilidad de éxito p", "E[T] = 1 / p", "Var[T] = (1 − p) / p²", "Desviación estándar",
             "Mediana: qgeom(0,5, p) + 1", "P(T ≤ 10) = pgeom(9, p)", "P(T > 50) = 1 − pgeom(49, p)",
             "P(T > 50) a mano: (1 − p)^50"),
  T = c(fpct(vpp), fnum(gT[["E"]], 2), fnum(gT[["V"]], 2), fnum(sqrt(gT[["V"]]), 2), gT[["med"]],
        fpct(gT[["P10"]]), fprob(gT[["P50"]]), fprob(gT[["P50_mano"]])),
  T0 = c(fpct(pE), fnum(gT0[["E"]], 2), fnum(gT0[["V"]], 2), fnum(sqrt(gT0[["V"]]), 2), gT0[["med"]],
         fpct(gT0[["P10"]]), fprob(gT0[["P50"]]), fprob(gT0[["P50_mano"]])))
setnames(res_geo, c("Medida", "T: con marcador M", "T₀: sin marcador"))
tabla_html(res_geo, caption = "Esperanza, varianza, mediana y probabilidades de cola de T y T₀", label = "resT")

Tabla 17. Esperanza, varianza, mediana y probabilidades de cola de T y T₀
Medida	T: con marcador M	T₀: sin marcador
Probabilidad de éxito p	26,60%	5,99%
E[T] = 1 / p	3,76	16,70
Var[T] = (1 − p) / p²	10,38	262,26
Desviación estándar	3,22	16,19
Mediana: qgeom(0,5, p) + 1	3	12
P(T ≤ 10) = pgeom(9, p)	95,46%	46,07%
P(T > 50) = 1 − pgeom(49, p)	0,0000193%	4,56%
P(T > 50) a mano: (1 − p)^50	0,0000193%	4,56%

Interpretación (Tabla 16, Tabla 17 y Figura 7): con el marcador, el número esperado de fichas revisadas hasta el primer caso es E[T] = 1 / 0,2660 = 3,76, la mediana es 3 y 95,5% de las búsquedas termina en 10 fichas o menos. P(T > 50) = (1 − VPP)^50 = 0,0000193%: es prácticamente imposible revisar 50 fichas marcadas sin encontrar un caso. Sin marcador, E[T₀] = 16,70 fichas, la mediana es 12 y P(T₀ > 50) = 4,6%. Gracias al marcador, el auditor revisa en promedio 12,9 fichas menos por caso encontrado (4,4 veces menos). La geométrica tiene su moda en k = 1 y decrece siempre, así que la media es mayor que la mediana; la desviación estándar de T₀ (16,2) muestra que sin marcador una búsqueda puede alargarse mucho.

6.2 T₅ ~ Binomial negativa(5, VPP): revisar hasta encontrar 5 casos

E5 <- 5 / vpp; V5 <- 5 * (1 - vpp) / vpp^2
P5_100 <- pnbinom(100 - 5, size = 5, prob = vpp)          # P(T5 ≤ 100): a lo más 95 fracasos
E5_0 <- 5 / pE; P5_100_0 <- pnbinom(100 - 5, size = 5, prob = pE)
med5 <- qnbinom(0.5, 5, vpp) + 5
nb_tab <- data.table(
  Medida = c("E[T₅] = 5 / p", "Var[T₅] = 5(1 − p) / p²", "Mediana: qnbinom(0,5, 5, p) + 5",
             "P(T₅ ≤ 100) = pnbinom(95, 5, p)", "P(T₅ > 100)"),
  `Con marcador (p = VPP)` = c(fnum(E5, 2), fnum(V5, 2), med5, fprob_alta(P5_100), fprob(1 - P5_100)),
  `Sin marcador (p = P(E))` = c(fnum(E5_0, 2), fnum(5 * (1 - pE) / pE^2, 2), qnbinom(0.5, 5, pE) + 5,
                                fprob_alta(P5_100_0), fprob(1 - P5_100_0)))
tabla_html(nb_tab, caption = "Fichas revisadas hasta encontrar 5 episodios con I50: binomial negativa con y sin marcador", label = "nb")

Tabla 18. Fichas revisadas hasta encontrar 5 episodios con I50: binomial negativa con y sin marcador
Medida	Con marcador (p = VPP)	Sin marcador (p = P(E))
E[T₅] = 5 / p	18,80	83,51
Var[T₅] = 5(1 − p) / p²	51,89	1.311,29
Mediana: qnbinom(0,5, 5, p) + 5	18	78
P(T₅ ≤ 100) = pnbinom(95, 5, p)	99,99999971%	72,15%
P(T₅ > 100)	0,000000284%	27,85%

Interpretación (Tabla 18): para reunir 5 casos con el marcador se esperan E[T₅] = 5 / 0,2660 = 18,80 fichas, cinco veces E[T], porque T₅ es la suma de cinco tiempos geométricos independientes. P(T₅ ≤ 100) = 99,99999971%: revisar 100 fichas marcadas basta casi con certeza (la probabilidad de no lograrlo es 0,000000284%). Sin marcador se esperarían 83,5 fichas, y 100 fichas alcanzarían solo en 72,1% de las búsquedas.

6.3 Verificación con los datos: revisar fichas marcadas en orden aleatorio

Se toman los 98.239 episodios marcados por M, se desordenan al azar y se anota la posición del primer episodio con I50 y la del quinto. Se repite 1.000 veces. En los datos la revisión es sin reposición: una ficha revisada no vuelve a aparecer.

set.seed(20261030)
eM <- omega$E[omega$M]                                  # TRUE si el episodio marcado tiene I50
R_rev <- 1000
pos <- t(replicate(R_rev, {
  orden <- sample.int(length(eM))                       # una permutación aleatoria de las fichas marcadas
  w <- which(eM[orden])
  c(primero = w[1], quinto = w[5])
}))
E_sinrep <- (nM + 1) / (nME + 1)                       # esperanza exacta del primer éxito sin reposición
# error estándar aproximado de la varianza muestral de una geométrica: Var · √((κ − 1) / n),
# con κ = 9 + p² / (1 − p) la curtosis de la geométrica
se_var_geo <- function(pp, n) (1 - pp) / pp^2 * sqrt((8 + pp^2 / (1 - pp)) / n)
verif <- data.table(
  Variable = c("T (primer caso)", "T₅ (quinto caso)"),
  `Teórico con reposición` = c(paste0("1 / VPP = ", fnum(1 / vpp, 4)), paste0("5 / VPP = ", fnum(5 / vpp, 4))),
  `Teórico sin reposición` = c(paste0("(n(M) + 1) / (n(M ∩ E) + 1) = ", fnum(E_sinrep, 4)),
                               paste0("5 (n(M) + 1) / (n(M ∩ E) + 1) = ", fnum(5 * E_sinrep, 4))),
  `Promedio en 1.000 revisiones desordenadas` = fnum(colMeans(pos), 4),
  `Varianza teórica (geom. / BN)` = fnum(c((1 - vpp) / vpp^2, 5 * (1 - vpp) / vpp^2), 3),
  `Varianza en las revisiones` = fnum(apply(pos, 2, var), 3),
  `Máximo observado` = apply(pos, 2, max))
tabla_html(verif, caption = paste0("Posición del primer y del quinto episodio con I50 al revisar en orden aleatorio los ", fnum(nM), " episodios marcados (1.000 repeticiones)"),
           label = "verifT")

Tabla 19. Posición del primer y del quinto episodio con I50 al revisar en orden aleatorio los 98.239 episodios marcados (1.000 repeticiones)
Variable	Teórico con reposición	Teórico sin reposición	Promedio en 1.000 revisiones desordenadas	Varianza teórica (geom. / BN)	Varianza en las revisiones	Máximo observado
T (primer caso)	1 / VPP = 3,7601	(n(M) + 1) / (n(M ∩ E) + 1) = 3,7600	3,9100	10,378	12,218	32
T₅ (quinto caso)	5 / VPP = 18,8003	5 (n(M) + 1) / (n(M ∩ E) + 1) = 18,7998	19,1460	51,890	55,991	54

ecdf_T <- data.table(k = 1:max(15, max(pos[, 1])))
ecdf_T[, obs := sapply(k, function(x) mean(pos[, 1] == x))][, teo := dgeom(k - 1, vpp)]

ggplot(ecdf_T, aes(x = k)) +
  geom_col(aes(y = obs), fill = paleta[4], width = 0.8) +
  geom_point(aes(y = teo), color = paleta[2], size = 2.6) +
  geom_line(aes(y = teo), color = paleta[2], linewidth = 0.6) +
  scale_x_continuous(breaks = ecdf_T$k) +
  scale_y_continuous(labels = lab_pct) +
  labs(title = "Primer caso al revisar fichas marcadas: datos contra geométrica",
       subtitle = "Barras: 1.000 revisiones desordenadas de los episodios marcados · Puntos: dgeom(k − 1, VPP)",
       x = "k (posición del primer episodio con I50)", y = "Frecuencia relativa / probabilidad")

Figura 8. Frecuencia relativa de la posición del primer episodio con I50 en las 1.000 revisiones desordenadas (barras) y probabilidad de la geométrica con p = VPP (puntos).

Interpretación (Tabla 19 y Figura 8): en las 1.000 revisiones desordenadas, el primer caso apareció en promedio en la ficha 3,91, frente a 1 / VPP = 3,76, y el quinto en la 19,15, frente a 18,80. La frecuencia de cada posición sigue de cerca a la geométrica (Figura 8). La diferencia entre revisar con o sin reposición es despreciable: sin reposición la esperanza exacta es (n(M) + 1) / (n(M ∩ E) + 1) = 3,7600, apenas 0,0001 fichas menos que 1 / VPP, porque una búsqueda revisa pocas fichas de un total de 98.239 y casi no cambia la proporción de casos entre las que quedan. La diferencia que se ve entre el promedio observado y 1 / VPP (0,150) es del orden del azar de 1.000 repeticiones: el error estándar del promedio es √(Var / 1.000) = 0,102. La varianza de las revisiones (12,22) también difiere de la teórica (10,38). Con solo 1.000 repeticiones la varianza de una distribución con cola larga fluctúa bastante: su error estándar aproximado es 0,93. Por eso el Resultado 06 repite la revisión desordenada 10.000 veces.

7 Poisson: egresos por día

RESULTADO 05

D cuenta los episodios de la cohorte con fecha de egreso (FECHAALTA) en cada uno de los 366 días de 2024. La serie se construye uniendo el calendario completo con los conteos y reemplazando por 0 los días sin egresos, de modo que tiene exactamente 366 días.

serie <- data.table(
  Medida = c("Días en la serie", "Días sin egresos de la cohorte (D = 0)", "Episodios de E con egreso fuera de 2024",
             "Suma de D en los 366 días", "λ estimado = promedio diario", "Mínimo y máximo de D",
             "Media de D", "Varianza de D (var())", "Índice de dispersión = varianza / media"),
  Valor = c(fnum(nrow(dias)), fnum(sum(dias$D == 0)), fnum(n_fuera_2024), fnum(sum(dias$D)),
            paste0(fnum(sum(dias$D)), " / 366 = ", fnum(lambda, 2), " egresos por día"),
            paste0(min(dias$D), " y ", max(dias$D)), fnum(mean(dias$D), 2), fnum(var(dias$D), 2),
            fnum(var(dias$D) / mean(dias$D), 2)))
tabla_html(serie, caption = "Serie diaria D de egresos de la cohorte en 2024: construcción y resumen", label = "serie")

Tabla 20. Serie diaria D de egresos de la cohorte en 2024: construcción y resumen
Medida	Valor
Días en la serie	366
Días sin egresos de la cohorte (D = 0)	0
Episodios de E con egreso fuera de 2024	0
Suma de D en los 366 días	52.312
λ estimado = promedio diario	52.312 / 366 = 142,93 egresos por día
Mínimo y máximo de D	29 y 288
Media de D	142,93
Varianza de D (var())	3.249,78
Índice de dispersión = varianza / media	22,74

ID_todo <- var(dias$D) / mean(dias$D)

Interpretación (Tabla 20): la serie tiene los 366 días de 2024, y todos los egresos de la cohorte caen en 2024 (0 fuera de rango). En este caso ningún día quedó con D = 0, porque la cohorte es grande: hay al menos 29 egresos cada día. Aun así, el código completa el calendario para que un día sin egresos no desaparezca de la serie. El promedio es λ = 142,93 egresos por día, muy por encima de un egreso diario, así que no hace falta pasar a conteos semanales.

7.1 Frecuencias observadas contra esperadas por la Poisson

D toma más de 200 valores distintos, y cada valor aparece en muy pocos días. Por eso las frecuencias se agrupan en intervalos de 20 egresos. La frecuencia esperada de cada intervalo [a, b] es 366 · P(a ≤ D ≤ b) = 366 · (ppois(b, λ) − ppois(a − 1, λ)), y el último intervalo queda abierto para que las esperadas sumen 366.

ancho <- 20
cortes <- seq(0, ceiling((max(dias$D) + 1) / ancho) * ancho, by = ancho)
tabla_frec <- function(x, m, size_nb = NULL) {
  a <- head(cortes, -1); b <- c(tail(a, -1) - 1, Inf)
  obs <- sapply(seq_along(a), function(i) sum(x >= a[i] & x <= b[i]))
  # sobre la media se restan colas superiores, para no perder precisión al restar números cercanos a 1
  pp <- ifelse(a > m, ppois(a - 1, m, lower.tail = FALSE) - ppois(b, m, lower.tail = FALSE),
               ppois(b, m) - ppois(a - 1, m))
  out <- data.table(int = ifelse(is.infinite(b), paste0(a, " o más"), paste0(a, "–", b)), a = a,
                    obs = obs, p_obs = obs / length(x), p_poi = pp)
  if (!is.null(size_nb)) {
    pn <- ifelse(a > m, pnbinom(a - 1, size = size_nb, mu = m, lower.tail = FALSE) -
                   pnbinom(b, size = size_nb, mu = m, lower.tail = FALSE),
                 pnbinom(b, size = size_nb, mu = m) - pnbinom(a - 1, size = size_nb, mu = m))
    out[, p_nb := pn]
  }
  out[]
}
m_todo <- mean(dias$D); v_todo <- var(dias$D)
size_todo <- m_todo^2 / (v_todo - m_todo)               # binomial negativa por momentos
fr_todo <- tabla_frec(dias$D, m_todo, size_todo)
tvd <- function(po, pm) 0.5 * sum(abs(po - pm))
tabla_html(fr_todo[, .(`Intervalo de D` = int, `Días observados` = obs,
                       `Esperados Poisson 366 · P` = fsig(366 * p_poi),
                       `Proporción observada` = fpct(p_obs, 1), `Probabilidad Poisson` = fprob(p_poi))],
           caption = paste0("Frecuencias observadas de D en los 366 días de 2024 y esperadas por Poisson(λ = ", fnum(lambda, 2), "), en intervalos de ", ancho, " egresos"),
           label = "frecD")

Tabla 21. Frecuencias observadas de D en los 366 días de 2024 y esperadas por Poisson(λ = 142,93), en intervalos de 20 egresos
Intervalo de D	Días observados	Esperados Poisson 366 · P	Proporción observada	Probabilidad Poisson
0–19	0	< 0,0000000001	0,0%	< 0,0000000001%
20–39	7	< 0,0000000001	1,9%	< 0,0000000001%
40–59	36	< 0,0000000001	9,8%	< 0,0000000001%
60–79	36	0,00000136	9,8%	0,000000373%
80–99	35	0,0235	9,6%	0,00641%
100–119	3	8,25	0,8%	2,25%
120–139	8	135,22	2,2%	36,95%
140–159	55	191,51	15,0%	52,33%
160–179	70	30,42	19,1%	8,31%
180–199	69	0,57	18,9%	0,16%
200–219	35	0,00140	9,6%	0,000381%
220–239	8	0,000000506	2,2%	0,000000138%
240–259	2	< 0,0000000001	0,5%	< 0,0000000001%
260–279	1	< 0,0000000001	0,3%	< 0,0000000001%
280 o más	1	< 0,0000000001	0,3%	< 0,0000000001%

fig_D <- melt(fr_todo[, .(a, Poisson = 366 * p_poi, `Binomial negativa` = 366 * p_nb)], id.vars = "a",
              variable.name = "modelo", value.name = "esp")
ggplot(fr_todo, aes(x = a + ancho / 2)) +
  geom_col(aes(y = obs), fill = paleta[4], width = ancho * 0.9) +
  geom_line(data = fig_D, aes(y = esp, color = modelo), linewidth = 0.9) +
  geom_point(data = fig_D, aes(y = esp, color = modelo), size = 2.4) +
  scale_color_manual(values = c(Poisson = paleta[2], `Binomial negativa` = paleta[3]), name = "Esperadas") +
  scale_x_continuous(breaks = cortes) +
  labs(title = "Egresos diarios de la cohorte: observado contra esperado",
       subtitle = paste0("λ = ", fnum(lambda, 1), " · varianza observada = ", fnum(v_todo, 0)),
       x = "D (egresos de la cohorte en el día)", y = "Número de días") +
  theme(legend.position = "bottom")

Figura 9. Frecuencias observadas de D (barras) y esperadas por la Poisson y por la binomial negativa ajustada por momentos (puntos), en los 366 días de 2024.

Interpretación (Tabla 21 y Figura 9): la Poisson no describe la serie. Con λ = 142,9 concentra casi todos los días entre 113 y 175 egresos (99% central), pero solo 29,8% de los días observados cae en ese rango. Hay muchos días con menos de 100 egresos y bastantes con más de 200, que para la Poisson serían prácticamente imposibles (por ejemplo, la Poisson espera 0,0235 días con menos de 100 egresos y se observan 114). Además, la distribución observada tiene dos grupos: uno alrededor de 70 y otro alrededor de 170 egresos, que se explican en la sección siguiente.

7.2 Media, varianza e índice de dispersión

En una Poisson la media y la varianza son iguales a λ, así que el índice de dispersión (varianza / media) vale 1. Si es cercano a 1, la variación diaria es la que espera la Poisson. Si es mayor que 1 hay sobredispersión: los días varían más de lo que permite un λ constante, típicamente porque λ cambia de un día a otro o porque los egresos no son independientes. Si es menor que 1 hay subdispersión: los conteos son más regulares que lo esperado.

En la serie, la media es 142,93 y la varianza 3.249,78, así que el índice de dispersión es 22,74 (Tabla 20): la varianza observada es 22,7 veces la que supone la Poisson. Hay una sobredispersión fuerte.

7.3 Promedio de D por día de la semana

dias_sem <- c("Lunes", "Martes", "Miércoles", "Jueves", "Viernes", "Sábado", "Domingo")
dias[, dow := factor(dias_sem[as.integer(format(fecha, "%u"))], levels = dias_sem)]
# feriados legales de Chile en 2024 (lunes a viernes); se usan para definir días hábiles
feriados <- as.Date(c("2024-01-01", "2024-03-29", "2024-05-01", "2024-05-21", "2024-06-20",
                      "2024-07-16", "2024-08-15", "2024-09-18", "2024-09-19", "2024-09-20",
                      "2024-10-31", "2024-11-01", "2024-12-25"))
dias[, feriado := fecha %in% feriados]
dias[, habil := !(dow %in% c("Sábado", "Domingo")) & !feriado]
res_dow <- dias[, .(n_dias = .N, media = mean(D), varianza = var(D)), by = dow][order(dow)]
res_dow[, ID := varianza / media]
tabla_html(res_dow[, .(`Día de la semana` = dow, `Días` = n_dias, `Promedio de D` = fnum(media, 1),
                       `Varianza de D` = fnum(varianza, 1), `Índice de dispersión` = fnum(ID, 2))],
           caption = "Promedio, varianza e índice de dispersión de D por día de la semana (cada fila sobre los días de ese tipo en 2024)",
           label = "dow")

Tabla 22. Promedio, varianza e índice de dispersión de D por día de la semana (cada fila sobre los días de ese tipo en 2024)
Día de la semana	Días	Promedio de D	Varianza de D	Índice de dispersión
Lunes	53	167,2	874,4	5,23
Martes	53	178,8	1.171,9	6,56
Miércoles	52	170,8	1.164,1	6,81
Jueves	52	163,5	868,4	5,31
Viernes	52	186,7	1.339,8	7,18
Sábado	52	80,7	115,3	1,43
Domingo	52	51,7	127,8	2,47

m_fer <- dias[feriado == TRUE, mean(D)]

ggplot(dias, aes(x = dow, y = D)) +
  geom_boxplot(fill = paleta[4], color = paleta[3], outlier.shape = NA, width = 0.6) +
  geom_jitter(data = dias[feriado == FALSE], width = 0.15, alpha = 0.35, size = 1.2, color = paleta[3]) +
  geom_point(data = dias[feriado == TRUE], color = paleta[2], size = 2.6) +
  geom_hline(yintercept = lambda, linetype = "dashed", color = paleta[6]) +
  labs(title = "Egresos diarios de la cohorte según el día de la semana",
       subtitle = paste0("Línea discontinua: λ = ", fnum(lambda, 1), " (promedio de los 366 días)"),
       x = NULL, y = "D (egresos en el día)")

Figura 10. Distribución de D por día de la semana. Los puntos naranjos son feriados que caen de lunes a viernes.

Interpretación (Tabla 22 y Figura 10): los fines de semana se comportan muy distinto. De lunes a viernes el promedio va de 163,5 a 186,7 egresos, mientras que el sábado baja a 80,7 y el domingo a 51,7. Los 13 feriados que caen de lunes a viernes se parecen a un fin de semana (promedio 72,7). Mezclar días con promedios tan distintos es lo que produce los dos grupos de Figura 9 y buena parte de la sobredispersión. Por eso se repite la comparación solo con los días hábiles (lunes a viernes que no son feriado: 249 días). Dentro de cada día de la semana el índice de dispersión sigue siendo mayor que 1 (entre 1,4 y 7,2), así que el día de la semana no explica toda la variación.

7.4 Días hábiles y binomial negativa ajustada por momentos

Como hay sobredispersión, se ajusta una binomial negativa por momentos: con media μ y varianza σ² observadas, size = μ² / (σ² − μ), y en R dnbinom(k, size = size, mu = μ). Su varianza es μ + μ² / size = σ², así que reproduce la media y la varianza observadas. Los dos modelos se comparan con la distancia de variación total, DVT = ½ Σ |proporción observada − probabilidad esperada|, calculada sobre los mismos intervalos de 20 egresos. La DVT va de 0 (ajuste perfecto) a 1 (sin ninguna coincidencia).

Dh <- dias[habil == TRUE, D]
m_h <- mean(Dh); v_h <- var(Dh); size_h <- m_h^2 / (v_h - m_h)
fr_h <- tabla_frec(Dh, m_h, size_h)
tabla_html(fr_h[obs > 0 | length(Dh) * p_poi >= 0.01 | length(Dh) * p_nb >= 0.01,
                .(`Intervalo de D` = int, `Días observados` = obs,
                  `Esperados Poisson` = fsig(length(Dh) * p_poi),
                  `Esperados binomial negativa` = fsig(length(Dh) * p_nb))],
           caption = paste0("Frecuencias observadas de D en los ", length(Dh), " días hábiles y esperadas por Poisson(", fnum(m_h, 2), ") y por la binomial negativa (size = ", fnum(size_h, 2), ")"),
           label = "frecDh")

Tabla 23. Frecuencias observadas de D en los 249 días hábiles y esperadas por Poisson(178,64) y por la binomial negativa (size = 73,29)
Intervalo de D	Días observados	Esperados Poisson	Esperados binomial negativa
80–99	0	0,0000000136	0,0288
100–119	0	0,000328	1,07
120–139	8	0,30	11,21
140–159	55	18,14	43,56
160–179	70	113,66	75,99
180–199	69	101,63	68,24
200–219	35	14,89	35,02
220–239	8	0,38	11,15
240–259	2	0,00179	2,35
260–279	1	0,00000173	0,35
280 o más	1	0,000000000378	0,0406

comp <- data.table(
  Serie = c("Todos los días (366)", paste0("Días hábiles (", length(Dh), ")")),
  Media = fnum(c(m_todo, m_h), 2), Varianza = fnum(c(v_todo, v_h), 2),
  `Índice de dispersión` = fnum(c(v_todo / m_todo, v_h / m_h), 2),
  `size = μ² / (σ² − μ)` = fnum(c(size_todo, size_h), 2),
  `DVT Poisson` = fnum(c(tvd(fr_todo$p_obs, fr_todo$p_poi), tvd(fr_h$p_obs, fr_h$p_poi)), 3),
  `DVT binomial negativa` = fnum(c(tvd(fr_todo$p_obs, fr_todo$p_nb), tvd(fr_h$p_obs, fr_h$p_nb)), 3))
tabla_html(comp, caption = "Poisson contra binomial negativa ajustada por momentos: índice de dispersión y distancia de variación total (DVT, sobre intervalos de 20 egresos)",
           label = "dvt")

Tabla 24. Poisson contra binomial negativa ajustada por momentos: índice de dispersión y distancia de variación total (DVT, sobre intervalos de 20 egresos)
Serie	Media	Varianza	Índice de dispersión	size = μ² / (σ² − μ)	DVT Poisson	DVT binomial negativa
Todos los días (366)	142,93	3.249,78	22,74	6,58	0,735	0,353
Días hábiles (249)	178,64	614,05	3,44	73,29	0,306	0,055

mes_h <- dias[habil == TRUE, .(media = mean(D)), by = .(mes = as.integer(format(fecha, "%m")))]
meses <- c("enero", "febrero", "marzo", "abril", "mayo", "junio", "julio", "agosto", "septiembre", "octubre", "noviembre", "diciembre")
dvt_poi_h <- tvd(fr_h$p_obs, fr_h$p_poi); dvt_nb_h <- tvd(fr_h$p_obs, fr_h$p_nb)
dvt_poi_t <- tvd(fr_todo$p_obs, fr_todo$p_poi); dvt_nb_t <- tvd(fr_todo$p_obs, fr_todo$p_nb)

fig_Dh <- melt(fr_h[, .(a, Poisson = length(Dh) * p_poi, `Binomial negativa` = length(Dh) * p_nb)], id.vars = "a",
               variable.name = "modelo", value.name = "esp")
ggplot(fr_h, aes(x = a + ancho / 2)) +
  geom_col(aes(y = obs), fill = paleta[4], width = ancho * 0.9) +
  geom_line(data = fig_Dh, aes(y = esp, color = modelo), linewidth = 0.9) +
  geom_point(data = fig_Dh, aes(y = esp, color = modelo), size = 2.4) +
  scale_color_manual(values = c(Poisson = paleta[2], `Binomial negativa` = paleta[3]), name = "Esperadas") +
  scale_x_continuous(breaks = cortes) +
  labs(title = "Egresos de la cohorte en días hábiles: observado contra esperado",
       subtitle = paste0("Media = ", fnum(m_h, 1), " · varianza = ", fnum(v_h, 0), " · índice de dispersión = ", fnum(v_h / m_h, 2)),
       x = "D (egresos de la cohorte en el día hábil)", y = "Número de días") +
  theme(legend.position = "bottom")

Figura 11. Días hábiles: frecuencias observadas de D (barras) y esperadas por la Poisson y por la binomial negativa ajustada por momentos.

Interpretación (Tabla 23, Tabla 24 y Figura 11): al dejar solo los días hábiles el índice de dispersión baja de 22,7 a 3,4, pero sigue lejos de 1: incluso entre días hábiles la varianza es 3,4 veces la media. Una razón visible en los datos es que el promedio de los días hábiles cambia a lo largo del año, desde 151,5 en febrero hasta 198,9 en julio, de modo que λ no es constante. La binomial negativa describe mucho mejor la serie: su DVT es 0,055 en días hábiles frente a 0,306 de la Poisson, y 0,353 frente a 0,735 con todos los días. Con todos los días la binomial negativa tampoco reproduce los dos grupos de Figura 9, porque es una distribución de una sola moda; en días hábiles, donde la distribución tiene una sola moda, el ajuste es bueno. Que la Poisson se ajuste mal muestra qué supuesto falla: los egresos de la cohorte no ocurren con una tasa constante a lo largo del año ni de la semana.

7.5 La Poisson como aproximación de la binomial

B es el número de episodios con I50 entre 1.000 episodios elegidos al azar de Ω: B ~ Binomial(1.000, P(E)), con P(E) = 5,99%. Como n es grande y P(E) pequeña, se puede aproximar con Poisson(1.000 · P(E)) = Poisson(59,87).

mB <- 1000 * pE
kB <- seq(floor(qbinom(0.0005, 1000, pE) / 5) * 5, ceiling(qbinom(0.9995, 1000, pE) / 5) * 5, by = 5)
apx <- data.table(k = kB, bin = dbinom(kB, 1000, pE), poi = dpois(kB, mB),
                  Fbin = pbinom(kB, 1000, pE), Fpoi = ppois(kB, mB))
tabla_html(apx[, .(k, `Binomial(1.000, P(E)): P(B = k)` = fprob(bin), `Poisson(1.000 · P(E)): P = k` = fprob(poi),
                   `Diferencia` = fpp_dif(bin - poi), `Binomial P(B ≤ k)` = fprob_cdf(Fbin), `Poisson P(≤ k)` = fprob_cdf(Fpoi))],
           caption = paste0("Binomial(1.000; ", f4(pE), ") y Poisson(", fnum(mB, 2), ") en una misma tabla (k cada 5 valores, cubriendo el 99,9% central)"),
           label = "apx")

Tabla 25. Binomial(1.000; 0,0599) y Poisson(59,87) en una misma tabla (k cada 5 valores, cubriendo el 99,9% central)
k	Binomial(1.000, P(E)): P(B = k)	Poisson(1.000 · P(E)): P = k	Diferencia	Binomial P(B ≤ k)	Poisson P(≤ k)
35	0,0113%	0,0154%	-0,004 p.p.	0,0251%	0,0353%
40	0,12%	0,15%	-0,026 p.p.	0,33%	0,42%
45	0,71%	0,79%	-0,070 p.p.	2,41%	2,76%
50	2,32%	2,38%	-0,061 p.p.	10,35%	11,07%
55	4,45%	4,38%	0,069 p.p.	28,44%	29,10%
60	5,30%	5,14%	0,162 p.p.	54,10%	54,08%
65	4,07%	3,99%	0,079 p.p.	77,60%	76,96%
70	2,08%	2,11%	-0,039 p.p.	91,91%	91,25%
75	0,72%	0,79%	-0,063 p.p.	97,86%	97,51%
80	0,18%	0,21%	-0,034 p.p.	99,58%	99,46%
85	0,0305%	0,0410%	-0,010 p.p.	99,94%	99,91%
90	0,00384%	0,00597%	-0,002 p.p.	99,9935%	99,99%

k_all <- 0:1000
dif_max <- max(abs(dbinom(k_all, 1000, pE) - dpois(k_all, mB)))
dvt_apx <- tvd(dbinom(k_all, 1000, pE), dpois(k_all, mB))

Interpretación (Tabla 25): las dos distribuciones tienen la misma esperanza, 59,87 episodios con I50 entre 1.000. La varianza de la binomial es 1.000 · P(E) · (1 − P(E)) = 56,29, algo menor que la de la Poisson (59,87), porque el factor (1 − P(E)) = 0,9401 no es despreciable. Por eso la Poisson pone un poco menos de probabilidad cerca de la media y un poco más en las colas. Aun así la aproximación es buena: la mayor diferencia entre ambas en un valor de k es 0,162 puntos porcentuales y la DVT sobre k = 0, ..., 1.000 es 0,0149. La aproximación mejora cuanto menor es P(E).

8 Teoría, simulación y datos

RESULTADO 06

Para S, T y D se comparan tres valores de la esperanza y de la varianza:

Teórico: la fórmula del modelo con los parámetros de la base.
Simulación del modelo: 10.000 realizaciones con rbinom(10000, 20, p), rgeom(10000, vpp) + 1 (convención de ensayos) y rpois(10000, λ).
Remuestreo de los datos reales: 10.000 auditorías de 20 episodios de la cohorte elegidos con reposición para S; la revisión desordenada de las fichas marcadas del Resultado 04 para T, repetida ahora 10.000 veces para que las tres variables usen el mismo número de repeticiones; y 10.000 días elegidos con reposición entre los 366 días observados para D.

Se agregan dos filas que ayudan a interpretar: T₅ (binomial negativa) y una auditoría de 20 episodios dentro de un establecimiento elegido al azar (con probabilidad proporcional a sus episodios), que es como se haría una auditoría real.

set.seed(20261023)                                   # semilla fija: mismos resultados en cada compilación
B_sim <- 10000
# ---- simulación de los modelos (funciones r) ----
sim_S  <- rbinom(B_sim, 20, p)
sim_T  <- rgeom(B_sim, vpp) + 1                      # rgeom cuenta fracasos: se suma 1 para contar fichas
sim_T5 <- rnbinom(B_sim, size = 5, prob = vpp) + 5
sim_D  <- rpois(B_sim, lambda)
sim_Sh <- rbinom(B_sim, 20, p)                       # el modelo binomial no distingue establecimientos
# ---- remuestreo de los datos reales ----
fal <- coh$Fal
rem_S <- replicate(B_sim, sum(fal[sample.int(nE, 20, replace = TRUE)]))
pos10 <- t(replicate(B_sim, {                       # misma revisión desordenada del Resultado 04
  w <- which(eM[sample.int(length(eM))])
  c(primero = w[1], quinto = w[5])
}))
rem_T <- pos10[, 1]; rem_T5 <- pos10[, 2]
rem_D <- sample(dias$D, B_sim, replace = TRUE)
fal_h <- split(coh$Fal, coh$COD_HOSPITAL)
h_b <- coh$COD_HOSPITAL[sample.int(nE, B_sim, replace = TRUE)]   # establecimiento ∝ número de episodios
rem_Sh <- vapply(h_b, function(h) { x <- fal_h[[h]]; sum(x[sample.int(length(x), 20, replace = TRUE)]) }, numeric(1))

fila_sim <- function(var, modelo, E_t, V_t, s, r) {
  data.table(Variable = var, Modelo = modelo, E_t = E_t, E_s = mean(s), E_r = mean(r),
             V_t = V_t, V_s = var(s), V_r = var(r), n_r = length(r))
}
res6 <- rbind(
  fila_sim("S", "Binomial(20, p)", 20 * p, 20 * p * (1 - p), sim_S, rem_S),
  fila_sim("T", "Geométrica(VPP)", 1 / vpp, (1 - vpp) / vpp^2, sim_T, rem_T),
  fila_sim("D", "Poisson(λ)", lambda, lambda, sim_D, rem_D),
  fila_sim("T₅", "Binomial negativa(5, VPP)", 5 / vpp, 5 * (1 - vpp) / vpp^2, sim_T5, rem_T5),
  fila_sim("S dentro de un establecimiento", "Binomial(20, p)", 20 * p, 20 * p * (1 - p), sim_Sh, rem_Sh))
res6[, `:=`(dE_s = (E_s - E_t) / E_t, dE_r = (E_r - E_t) / E_t, dV_s = (V_s - V_t) / V_t, dV_r = (V_r - V_t) / V_t)]
t6 <- res6[, .(Variable, Modelo, `E teórica` = fnum(E_t, 3), `E simulada` = fnum(E_s, 3),
               `E remuestreo` = fnum(E_r, 3), `Var teórica` = fnum(V_t, 3), `Var simulada` = fnum(V_s, 3),
               `Var remuestreo` = fnum(V_r, 3))]
tabla_html(t6, caption = "Esperanza y varianza: valor teórico, simulación de 10.000 realizaciones del modelo y remuestreo de los datos reales",
           label = "sim") |>
  add_header_above(c(" " = 2, "Esperanza" = 3, "Varianza" = 3)) |>
  row_spec(1:3, bold = TRUE)

Tabla 26. Esperanza y varianza: valor teórico, simulación de 10.000 realizaciones del modelo y remuestreo de los datos reales
		Esperanza			Varianza
Variable	Modelo	E teórica	E simulada	E remuestreo	Var teórica	Var simulada	Var remuestreo
S	Binomial(20, p)	2,000	1,981	1,998	1,800	1,802	1,814
T	Geométrica(VPP)	3,760	3,790	3,800	10,378	10,293	10,719
D	Poisson(λ)	142,929	142,987	143,067	142,929	142,025	3.240,634
T₅	Binomial negativa(5, VPP)	18,800	18,831	18,720	51,890	51,834	51,329
S dentro de un establecimiento	Binomial(20, p)	2,000	1,985	2,001	1,800	1,786	1,948

t6b <- res6[, .(Variable, `Remuestreos usados` = fnum(n_r),
                `E: simulación vs teoría` = paste0(fnum(100 * dE_s, 2), "%"),
                `E: remuestreo vs teoría` = paste0(fnum(100 * dE_r, 2), "%"),
                `Var: simulación vs teoría` = paste0(fnum(100 * dV_s, 2), "%"),
                `Var: remuestreo vs teoría` = paste0(fnum(100 * dV_r, 2), "%"))]
tabla_html(t6b, caption = "Diferencias relativas respecto del valor teórico, (valor − teórico) / teórico", label = "sim_dif")

Tabla 27. Diferencias relativas respecto del valor teórico, (valor − teórico) / teórico
Variable	Remuestreos usados	E: simulación vs teoría	E: remuestreo vs teoría	Var: simulación vs teoría	Var: remuestreo vs teoría
S	10.000	-0,95%	-0,07%	0,10%	0,81%
T	10.000	0,80%	1,05%	-0,82%	3,29%
D	10.000	0,04%	0,10%	-0,63%	2.167,30%
T₅	10.000	0,16%	-0,43%	-0,11%	-1,08%
S dentro de un establecimiento	10.000	-0,75%	0,07%	-0,79%	8,24%

# varianza entre establecimientos de la proporción de fallecimientos (ponderada por episodios)
p_h <- est$K_F / est$N_h; w_h <- est$N_h / sum(est$N_h)
var_ph <- sum(w_h * (p_h - sum(w_h * p_h))^2)
V_mezcla <- 20 * p * (1 - p) + 20 * 19 * var_ph
D_dow <- data.table(D = rem_D)

fig6 <- rbind(
  data.table(v = "S: fallecimientos entre 20", origen = "Simulación del modelo", x = sim_S),
  data.table(v = "S: fallecimientos entre 20", origen = "Remuestreo de los datos", x = rem_S),
  data.table(v = "T: fichas hasta el primer caso", origen = "Simulación del modelo", x = sim_T),
  data.table(v = "T: fichas hasta el primer caso", origen = "Remuestreo de los datos", x = rem_T),
  data.table(v = "D: egresos por día", origen = "Simulación del modelo", x = sim_D),
  data.table(v = "D: egresos por día", origen = "Remuestreo de los datos", x = rem_D))
fig6[, origen := factor(origen, levels = c("Simulación del modelo", "Remuestreo de los datos"))]
fig6[, v := factor(v, levels = unique(v))]
# S y T se muestran por valor (T hasta 30 fichas); D en intervalos de 10 egresos
fig6 <- fig6[!(v == "T: fichas hasta el primer caso" & x > 30)]
fig6[, xb := ifelse(v == "D: egresos por día", floor(x / 10) * 10 + 5, x)]
fr6 <- fig6[, .N, by = .(v, origen, xb)][, prop := N / sum(N), by = .(v, origen)]
ggplot(fr6, aes(x = xb, y = prop, fill = origen)) +
  geom_col(data = fr6[v != "D: egresos por día"], position = position_dodge(width = 0.8), width = 0.8) +
  geom_col(data = fr6[v == "D: egresos por día"], position = position_dodge(width = 9), width = 9) +
  facet_wrap(~ v, scales = "free", ncol = 3) +
  scale_fill_manual(values = c(paleta[2], paleta[3]), name = NULL) +
  scale_x_continuous(breaks = function(l) { b <- pretty(l); b[b == round(b)] }, labels = function(x) fnum(x)) +
  scale_y_continuous(labels = lab_pct) +
  labs(title = "Simulación de los modelos contra remuestreo de los datos reales",
       subtitle = "T se muestra hasta 30 fichas y D en intervalos de 10 egresos",
       x = "Valor de la variable", y = "Proporción de realizaciones") +
  theme(legend.position = "bottom")

Figura 12. Distribución de la simulación del modelo y del remuestreo de los datos reales para S, T y D.

Interpretación (Tabla 26, Tabla 27 y Figura 12):

Diferencias que se deben al azar de la simulación. Con 10.000 realizaciones, la media simulada se aleja de la teórica en unas pocas veces √(Var / 10.000): para S el error estándar es 0,013, para T 0,032 y para D 0,120. Las diferencias de la columna "simulación vs teoría" (Tabla 27) son de ese tamaño, menores a 1,0%, y cambiarían de signo con otra semilla. Muestran que los modelos están bien programados, no que describan bien los datos.
S: el modelo coincide con el remuestreo de toda la cohorte. Elegir 20 episodios con reposición de toda la cohorte cumple exactamente los supuestos de la binomial, así que el remuestreo da E = 1,998 y Var = 1,814, igual a la teoría salvo el azar. En cambio, dentro de un establecimiento la media se mantiene (2,001) pero la varianza sube a 1,948, 8,2% más que la binomial. Ese exceso no es azar: corresponde a la variación de p entre establecimientos. La varianza de una mezcla de binomiales es 20p(1 − p) + 20 · 19 · Var(p_h) = 1,955, cercana al valor remuestreado. Es el supuesto de p común que no se cumplía en el Resultado 02.
T: el modelo coincide con la revisión real sin reposición. En 10.000 revisiones desordenadas la media es 3,800 y la varianza 10,72, frente a 3,760 y 10,38 de la geométrica (diferencias de 1,05% y 3,29%). Con 10.000 repeticiones el error estándar aproximado de la varianza es 0,30, así que la diferencia en la varianza equivale a 1,2 errores estándar. Es una diferencia del orden del azar, como se esperaba: sin reposición la esperanza exacta casi no cambia (Resultado 04), porque una búsqueda revisa muy pocas de las fichas marcadas. La diferencia mayor que se vio con 1.000 revisiones en el Resultado 04 se debía al menor número de repeticiones. T₅ se comporta igual (-0,43% en la media y -1,08% en la varianza).
D: el modelo no describe los datos. La media remuestreada (143,07) coincide con λ, pero la varianza remuestreada es 3.240,6, frente a 142,9 de la Poisson y 142,0 de su simulación: 22,7 veces más. Esta diferencia no se reduce aumentando las simulaciones. Muestra que λ no es constante (Resultado 05) y confirma que comparar solo la media no basta para declarar un buen ajuste.

9 Síntesis para la jefatura

RESULTADO 07

9.1 Cinco hallazgos principales

k90 <- qgeom(0.9, vpp) + 1                          # menor k con P(T ≤ k) ≥ 0,90 (convención de ensayos)
k90_bucle <- which(pgeom((1:200) - 1, vpp) >= 0.9)[1]   # verificación con un recorrido directo
k90_0 <- qgeom(0.9, pE) + 1
datos_k90 <- mean(pos[, 1] <= k90)

En una auditoría de 20 episodios de insuficiencia cardíaca elegidos al azar, lo esperable son 2,0 fallecimientos (S ~ Binomial(20; 0,100)), y hay 12,2% de probabilidad de no encontrar ninguno. P(S ≥ 2) = 60,8% y el percentil 90 es 4. Verificable en Tabla 11 y Figura 4.
Los fallecimientos por establecimiento varían más de lo que permite una p común (binomial por establecimiento). 11 de 72 establecimientos tienen P(S_h ≥ observado) ≤ 5% y 12 tienen P(S_h ≤ observado) ≤ 5%, cuando con p común se esperarían unos 3,1 y 3,0; la razón de dispersión es 3,29. Una auditoría dentro de un establecimiento tiene varianza 1,95 en vez de 1,80. Esto describe diferencias en los casos atendidos, no la calidad de la atención. Verificable en Tabla 12, Figura 5 y Tabla 26.
Auditar sin reposición un establecimiento pequeño cambia las probabilidades (hipergeométrica). En el establecimiento pequeño (Hospital Dr. Exequiel González Cortés (Santiago, San Miguel); N = 26, K = 2, n / N = 38,5%), P(H = 0) = 36,9% frente a 44,9% con la binomial, y la varianza se reduce por el factor 0,640. En el de mayor volumen (Complejo Hospitalario Dr. Sótero del Río (Santiago, Puente Alto); n / N = 0,32%) la diferencia es menor a 0,1 puntos porcentuales. Verificable en Tabla 14, Tabla 15 y Figura 6.
El marcador M reduce de 16,7 a 3,8 el número esperado de fichas revisadas hasta encontrar un caso (geométrica). Con VPP = 26,6%, P(T > 50) = 0,0000193%, y para 5 casos se esperan 18,8 fichas (binomial negativa). La revisión real de las fichas marcadas en orden aleatorio dio un promedio de 3,91 en 1.000 repeticiones y de 3,80 en 10.000. Verificable en Tabla 17, Tabla 18, Tabla 19 y Tabla 26.
Los egresos diarios de la cohorte promedian λ = 142,9, pero varían mucho más que una Poisson. El índice de dispersión es 22,7 con todos los días y 3,4 en días hábiles; el sábado y el domingo promedian 81 y 52 egresos, frente a 173 de lunes a viernes. La binomial negativa reduce la DVT de 0,306 a 0,055 en días hábiles. Verificable en Tabla 20, Tabla 22, Tabla 24 y Figura 11.

9.2 Limitaciones

Los episodios no son ensayos independientes. 39,5% de los episodios de la cohorte pertenece a un beneficiario con más de un episodio, y los episodios de un mismo establecimiento comparten forma de atención y de registro. La binomial, la geométrica y la Poisson suponen independencia, así que sus probabilidades son una referencia para muestras al azar de toda la cohorte y no describen bien grupos de episodios relacionados.
Los parámetros no son constantes, y los datos no cumplen ese supuesto. p cambia entre establecimientos (razón de dispersión 3,29), y λ cambia según el día de la semana y el mes (índice de dispersión 3,4 aun en días hábiles). Por eso la Poisson no sirve para planificar egresos diarios, y las probabilidades de cola por establecimiento solo ordenan y describen: no evalúan la atención.
Los parámetros y las variables dependen del registro y se estimaron con los mismos datos que se comparan. X está truncada en 35 campos y depende de cuán exhaustivo es cada establecimiento al codificar diagnósticos secundarios. El marcador M se eligió en la Actividad 2 con estos mismos datos, así que su VPP de 26,6% es optimista para otro año o para otra base. Finalmente, la DVT depende de los intervalos de 20 egresos usados para agrupar.

9.3 ¿Cuántas fichas marcadas por M debe revisar un auditor?

Debe revisar 8 fichas marcadas por M para tener al menos 90% de probabilidad de encontrar un episodio con insuficiencia cardíaca. Es el menor k con P(T ≤ k) ≥ 0,90, con T ~ Geométrica(VPP = 26,6%): P(T ≤ k) = 1 − (1 − VPP)^k, así que P(T ≤ 7) = 88,5% todavía no alcanza el 90% y P(T ≤ 8) = 91,6% sí. En R es qgeom(0.9, vpp) + 1 = 8 (un recorrido directo de k = 1, 2, ... da el mismo valor, 8). En las 1.000 revisiones reales desordenadas, el primer caso apareció dentro de las primeras 8 fichas en 90,0% de las veces. Sin el marcador harían falta 38 fichas para el mismo 90%.

10 Referencias

Fondo Nacional de Salud (FONASA), Chile. Base de datos pública de egresos hospitalarios y cirugía mayor ambulatoria bajo el mecanismo GRD, año 2024. Portal de datos abiertos de FONASA. Archivo GRD_PUBLICO_2024.txt, descargado el 21 de agosto de 2026.
FONASA. Tablas maestras GRD (hospitales), archivo Tablas_maestras_GRD.xlsx.
Organización Mundial de la Salud. Clasificación Internacional de Enfermedades, décima revisión (CIE-10), archivo CIE-10.xlsx.
Feriados legales de Chile 2024, según la legislación vigente, usados para definir los días hábiles.
R Core Team. R: A Language and Environment for Statistical Computing, funciones dbinom, dhyper, dgeom, dnbinom y dpois (paquete stats).
Carvajal, A. (2026). Caracterización estadística de la insuficiencia cardíaca en la base GRD 2024. Actividad 1, Estadística Computacional I, USACH.
Carvajal, A. (2026). Probabilidades y teorema de Bayes en los episodios GRD 2024 con insuficiencia cardíaca. Actividad 2, Estadística Computacional I, USACH.

Variables aleatorias y distribuciones discretas en los episodios GRD 2024 con insuficiencia cardíaca

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

Alejandro Carvajal

15 de octubre de 2026

1 Introducción

1.1 Contexto y objetivo

1.2 Declaración de uso de inteligencia artificial

2 Datos, correcciones y variables

2.1 Fuente de datos y tablas auxiliares

2.2 Correcciones a partir de la retroalimentación de las Actividades 1 y 2

2.3 Tabla de variables aleatorias

2.4 Parámetros estimados desde la base

2.5 Convención de R para la geométrica y la binomial negativa

3 Variable aleatoria empírica: número de diagnósticos

3.1 Función de probabilidad y función de distribución acumulada

3.2 Esperanza, varianza y medidas de posición

3.3 Comparación según el rol diagnóstico: Eₚ contra Eₛ

4 Bernoulli y binomial: fallecimientos en una auditoría de 20 episodios

4.1 La variable Bernoulli Y

4.2 S ~ Binomial(20, p) y S_L ~ Binomial(20, P(L | E))

4.3 Fallecimientos observados y esperados por establecimiento

4.4 Supuestos de la binomial en estos datos

5 Hipergeométrica: auditar sin reposición en un establecimiento

6 Geométrica y binomial negativa: revisar fichas hasta encontrar casos

6.1 Función de probabilidad de T y T₀ para k = 1 a 50

6.2 T₅ ~ Binomial negativa(5, VPP): revisar hasta encontrar 5 casos

6.3 Verificación con los datos: revisar fichas marcadas en orden aleatorio

7 Poisson: egresos por día

7.1 Frecuencias observadas contra esperadas por la Poisson

7.2 Media, varianza e índice de dispersión

7.3 Promedio de D por día de la semana

7.4 Días hábiles y binomial negativa ajustada por momentos

7.5 La Poisson como aproximación de la binomial

8 Teoría, simulación y datos

9 Síntesis para la jefatura

9.1 Cinco hallazgos principales

9.2 Limitaciones

9.3 ¿Cuántas fichas marcadas por M debe revisar un auditor?

10 Referencias

Variables aleatorias y distribuciones discretas en los episodios GRD 2024 con insuficiencia cardíaca

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

Alejandro Carvajal

15 de octubre de 2026

1 Introducción

1.1 Contexto y objetivo

1.2 Declaración de uso de inteligencia artificial

2 Datos, correcciones y variables

2.1 Fuente de datos y tablas auxiliares

2.2 Correcciones a partir de la retroalimentación de las Actividades 1 y 2

2.3 Tabla de variables aleatorias

2.4 Parámetros estimados desde la base

2.5 Convención de R para la geométrica y la binomial negativa

3 Variable aleatoria empírica: número de diagnósticos

3.1 Función de probabilidad y función de distribución acumulada

3.2 Esperanza, varianza y medidas de posición

3.3 Comparación según el rol diagnóstico: Eₚ contra Eₛ

4 Bernoulli y binomial: fallecimientos en una auditoría de 20 episodios

4.1 La variable Bernoulli Y

4.2 S ~ Binomial(20, p) y SL ~ Binomial(20, P(L | E))

4.3 Fallecimientos observados y esperados por establecimiento

4.4 Supuestos de la binomial en estos datos

5 Hipergeométrica: auditar sin reposición en un establecimiento

6 Geométrica y binomial negativa: revisar fichas hasta encontrar casos

6.1 Función de probabilidad de T y T₀ para k = 1 a 50

6.2 T₅ ~ Binomial negativa(5, VPP): revisar hasta encontrar 5 casos

6.3 Verificación con los datos: revisar fichas marcadas en orden aleatorio

7 Poisson: egresos por día

7.1 Frecuencias observadas contra esperadas por la Poisson

7.2 Media, varianza e índice de dispersión

7.3 Promedio de D por día de la semana

7.4 Días hábiles y binomial negativa ajustada por momentos

7.5 La Poisson como aproximación de la binomial

8 Teoría, simulación y datos

9 Síntesis para la jefatura

9.1 Cinco hallazgos principales

9.2 Limitaciones

9.3 ¿Cuántas fichas marcadas por M debe revisar un auditor?

10 Referencias

4.2 S ~ Binomial(20, p) y S_L ~ Binomial(20, P(L | E))