Modelo de VIH en HSH
Tu nombre
2025-11-27
1 Introducción
El objetivo de este documento es:
- Representar la dinámica de transmisión del VIH en HSH mediante un sistema de ecuaciones diferenciales ordinarias (EDO) estructurado en 12 compartimentos epidemiológicos para dos grupos de riesgo (24 estados).
- Estratificar la población HSH en bajo y alto riesgo (LRMSM y HRMSM).
- Incluir demografía, uso de condón, pruebas de VIH, tratamiento antirretroviral (TAR) y profilaxis preexposición (PrEP).
- Explorar escenarios de intervención (status quo, TasP, PrEP, combinados) y calcular indicadores epidemiológicos (prevalencia, incidencia, infecciones evitadas e indicadores de eliminación).
- Incorporar un análisis de sensibilidad e incertidumbre sobre la fuerza de infección y parámetros clave.
El horizonte temporal de la simulación se fija en 15 años con un paso
de integración de 0.01 años, consistente con el protocolo. El modelo es
determinístico y se resuelve numéricamente con el paquete
deSolve.
2 1. Paquetes
3 2. Compartimentos y nombres de estado
Se definen 12 compartimentos epidemiológicos:
S1 : Susceptible sin tamizaje (sin prueba reciente)
S2 : Susceptible prueba VIH (-) reciente (identificado en riesgo)
S3 : Susceptible en PrEP
A1 : Infección aguda sin identificar
A2 : Infección aguda identificada, sin TAR
A3 : Infección aguda en TAR
C1 : Infección crónica sin identificar
C2 : Infección crónica identificada, sin TAR
C3 : Infección crónica en TAR
D1 : SIDA sin identificar
D2 : SIDA identificado, sin TAR
D3 : SIDA en TAR
Cada compartimento se duplica para dos grupos de riesgo: LR (bajo) y HR (alto), para un total de 24 estados.
# 12 compartimentos epidemiológicos base:
comp_base <- c(
"S1","S2","S3",
"A1","A2","A3",
"C1","C2","C3",
"D1","D2","D3"
)
# Grupos de riesgo
risk_groups <- c("LR", "HR")
# Nombres completos de estados (24 en total)
state_names <- c(
paste0(comp_base, "_LR"),
paste0(comp_base, "_HR")
)
length(state_names)## [1] 24state_names## [1] "S1_LR" "S2_LR" "S3_LR" "A1_LR" "A2_LR" "A3_LR" "C1_LR" "C2_LR" "C3_LR"
## [10] "D1_LR" "D2_LR" "D3_LR" "S1_HR" "S2_HR" "S3_HR" "A1_HR" "A2_HR" "A3_HR"
## [19] "C1_HR" "C2_HR" "C3_HR" "D1_HR" "D2_HR" "D3_HR"# Índices
idx <- setNames(seq_along(state_names), state_names)4 3. Condiciones iniciales
Las condiciones iniciales distribuyen la población total de HSH entre susceptibles e infectados, respetando:
Tamaño total de hombres 18–59 años en Colombia.
Porcentaje de HSH en la población masculina.
Fracción de bajo y alto riesgo (LRMSM y HRMSM).
Prevalencia inicial en cada grupo de riesgo.
Proporción de infectados en infección aguda, crónica y SIDA, así como su estado de diagnóstico y TAR.
make_init <- function(
men_18_59 = 15140556, # hombres 18–59 (DANE)
msm_frac = 0.035, # % HSH (dentro del rango 0,0296–0,044)
lr_frac = 0.80, # % LRMSM dentro de HSH
prev_LR0 = 0.18, # prevalencia inicial LR (Asume teniendo en cuenta la prevalencia global de alrededor del 20%)
prev_HR0 = 0.30, # prevalencia inicial HR(Asume teniendo en cuenta la prevalencia global de alrededor del 20%)
#distribución de infectados por estadio
frac_stage = c(acute = 0.05, chronic = 0.75, aids = 0.20),
#dentro de cada estadio: sin diagnóstico / sin TAR / en TAR
frac_diag = c(undiagn = 0.20, noART = 0.15, ART = 0.65),#Se estima que 80,84% de las personas que viven con VIH conocen su diagnóstico.En TAR (del total de PVV): 65,27%;Entonces, diagnosticados sin TAR ≈ 80,84% – 65,27% ≈ 15,6%
#distribución de susceptibles en S1, S2, S3
frac_S = c(S1 = 0.79, S2 = 0.20, S3 = 0.01)# susceptibles sin prueba reciente;# prueba VIH negativa en últimos 12 meses-ASUME;# en PrEP (muy baja cobertura actual)
) {
#Población total HSH y por riesgo
N_hsh <- men_18_59 * msm_frac #total de HSH en Colombia (modelo)
N_LR <- N_hsh * lr_frac #HSH de bajo riesgo (LR)
N_HR <- N_hsh - N_LR #HSH de alto riesgo (HR), calculado como el complemento.
#Vector de estados en cero
init <- setNames(numeric(length(state_names)), state_names) #crea un vector de ceros, de longitud igual al número de compartimentos (24 si tienes 12×2).y le asigna el nombre.
#Susceptibles e infectados totales por grupo INICIALES
#misma distribución S1/S2/S3 se aplica a LR y HR
I_LR <- N_LR * prev_LR0
I_HR <- N_HR * prev_HR0
S_LR <- N_LR - I_LR
S_HR <- N_HR - I_HR
#Susceptibles LR
init["S1_LR"] <- S_LR * frac_S["S1"]
init["S2_LR"] <- S_LR * frac_S["S2"]
init["S3_LR"] <- S_LR * frac_S["S3"]
#Susceptibles HR
init["S1_HR"] <- S_HR * frac_S["S1"]
init["S2_HR"] <- S_HR * frac_S["S2"]
init["S3_HR"] <- S_HR * frac_S["S3"]
#Infectados: función auxiliar que reparte según estadio y diagnóstico
#dentro de cada estadio clínico (agudo, crónico, SIDA), las proporciones de diagnóstico/TAR (frac_diag) son iguales.
dist_infected <- function(init, I_total, prefix,
frac_stage, frac_diag) {
I_acute <- I_total * frac_stage["acute"]
I_chronic <- I_total * frac_stage["chronic"]
I_aids <- I_total * frac_stage["aids"]
fU <- frac_diag["undiagn"]
fN <- frac_diag["noART"]
fA <- frac_diag["ART"]
# Agudo A1–A3
init[paste0("A1_", prefix)] <- init[paste0("A1_", prefix)] + I_acute * fU
init[paste0("A2_", prefix)] <- init[paste0("A2_", prefix)] + I_acute * fN
init[paste0("A3_", prefix)] <- init[paste0("A3_", prefix)] + I_acute * fA
# Crónico C1–C3
init[paste0("C1_", prefix)] <- init[paste0("C1_", prefix)] + I_chronic * fU
init[paste0("C2_", prefix)] <- init[paste0("C2_", prefix)] + I_chronic * fN
init[paste0("C3_", prefix)] <- init[paste0("C3_", prefix)] + I_chronic * fA
# SIDA D1–D3
init[paste0("D1_", prefix)] <- init[paste0("D1_", prefix)] + I_aids * fU
init[paste0("D2_", prefix)] <- init[paste0("D2_", prefix)] + I_aids * fN
init[paste0("D3_", prefix)] <- init[paste0("D3_", prefix)] + I_aids * fA
init
}
#Aplicar dist_infected a LR y HR
init <- dist_infected(init, I_LR, "LR", frac_stage, frac_diag)
init <- dist_infected(init, I_HR, "HR", frac_stage, frac_diag)
stopifnot(all(names(init) == state_names)) #Aquí comprueba que el orden y los nombres de init coinciden exactamente con state_names
init
}
init <- make_init()
sum(init) # tamaño total de la población HSH modelada## [1] 529919.55 4. Parámetros base del modelo
Se definen parámetros demográficos, biológicos, conductuales y de intervención. En el modelo, los parámetros rho_LR y rho_HR son tasas de entrada demográfica (1/año) de nuevos HSH al sistema, específicas para los grupos de bajo y alto riesgo. Estos parámetros aparecen en las ecuaciones diferenciales de los susceptibles S1_LR y S1_HR como:
\(dS1_LR/dt = rho_LR * N_LR + ... - mu * S1_LR\) \(dS1_HR/dt = rho_HR * N_HR + ... - mu * S1_HR\) donde:
N_LR = suma de todos los compartimentos del grupo LR
N_HR = suma de todos los compartimentos del grupo HR
mu = tasa de mortalidad de población general (1/año)
Si agregamos todos los compartimentos de un grupo de riesgo g (LR o HR), y denotamos \(N_g(t) = población total del grupo g\), de forma esquemática podemos escribir (ignorando por un momento los flujos debidos al VIH): \(dN_g/dt ≈ rho_g * N_g - mu * N_g\) es decir: \(dN_g/dt ≈ (rho_g - mu) * N_g\)
A partir de esta expresión se deduce:
Si \(rho_g > mu\) -> \((rho_g - mu) > 0\) -> la población \(N_g\) crece en el tiempo.
Si \(rho_g = mu\) -> \((rho_g - mu) = 0\) -> N_g se mantiene aproximadamente constante.
Si \(rho_g < mu\) -> \((rho_g - mu) < 0\) -> N_g decrece en el tiempo.
Por lo tanto, rho_LR y rho_HR controlan el balance entre entradas demográficas (jóvenes que alcanzan 18 años, migración, etc.) y salidas por mortalidad dela población general en cada grupo de riesgo.
En el escenario base se fija: \(rho_LR = 0\) \(rho_HR = 0\) lo cual implica: \(dN_g/dt ≈ - mu * N_g < 0\) Es decir, se modela una COHORTE CERRADA de HSH sin entradas demográficas; sólo hay pérdidas por mortalidad (mu + mortalidad adicional por VIH).
Esta decisión tiene dos ventajas:
Simplifica la dinámica demográfica en la versión inicial del modelo y permite interpretar la evolución de la prevalencia e incidencia dentro de una cohorte fija.
Deja explícito que rho_LR y rho_HR son parámetros “calibrables”: en una fase posterior se pueden ajustar para representar una población aproximadamente estacionaria o con crecimiento neto.
Por ejemplo, si en ausencia de VIH se quisiera que el tamaño de cada grupo permanezca aproximadamente estable, una elección natural sería: \(rho_LR ≈ mu\) \(rho_HR ≈ mu\) lo que daría: \(dN_g/dt ≈ (mu - mu) * N_g = 0\).
En ese caso, el flujo de entrada demográfica compensaría las muertes de fondo. Adicionalmente, si se quiere compensar también la mortalidad extra debida al VIH (b4, …, b12), se pueden calibrar rho_LR y rho_HR ligeramente por encima de mu usando datos demográficos observados.
Conversión de propición observada de realización de pruebas en los últimos 12 meses en una tasa continua anual de prueba , bajo el supuesto que las pruebas ocurren como un proceso de poisson con tasa constante.
La probabilidad de NO hacerse ninguna prueba en 1 año es :
\[ P( 0 \mathrm{pruebas .en . 12. meses})=e^{\psi} \]
Entoces, la probabilidad de hacerse al menos una prueba en 12 meses es:
\[ P(\geqslant1 prueba) =1- e^{\psi} \]
\[-ln(1-p) =\psi\]
#Función auxiliar: convierte proporción tamizada en 12 meses
#Dada la proporción observada de personas que se hicieron al menos una prueba en 12 meses, asumo que las pruebas ocurren como un proceso de Poisson con tasa constante, y calculo la tasa continua anual (hazard) que genera esa proporción.
#en una tasa anual constante de pruebas (hazard)
derive_psi_from_prop <- function(p_test_12m) {
-log(1 - p_test_12m)
}
#Convierte una proporción observada de personas que se hicieron al menos una prueba en los últimos 12 meses (p_test_12m) en una tasa continua anual de prueba (psi), bajo el supuesto de que las pruebas ocurren como un proceso de Poisson con tasa constante.
parms_base <- list(
#Demografía y riesgo
men_18_59 = 15140556, # hombres 18–59
msm_frac = 0.035, # % HSH (dentro del rango 0,0296–0,044)
lr_frac = 0.80,
hr_frac = 0.20,
#Mortalidad general
mu = 0.0052, # 5,2 defunciones / 1000 hab-año
#Mortalidad específica por estadio (exceso sobre mu) (tasas anuales por individuo)
b4 = 0.0030, # Agudo sin TAR (A1)
b5 = 0.0030, # Agudo sin TAR (A2)
b6 = 0.0104, # Agudo con TAR (A3)
b7 = 0.0052, # Crónico sin TAR (C1)
b8 = 0.0052, # Crónico sin TAR (C2)
b9 = 0.0104, # Crónico con TAR (C3)
b10 = 1.68, # SIDA sin TAR (D1)
b11 = 1.68, # SIDA sin TAR (D2)
b12 = 0.05, # SIDA con TAR (D3)
#Tasa de entrada demográfica por grupo (calibrable)
rho_LR = 0.0,
rho_HR = 0.0,
#Progresión VIH (duraciones en años, protocolo)
dur_acute_y = 1.19, # Agudo -> Crónico
dur_chronic_y = 7.93, # Crónico -> SIDA
dur_aids_y = 2.0, # SIDA -> muerte (supuesto)
theta_acute = 1 / 1.19,
theta_chronic = 1 / 7.93,
theta_aids = 1 / 2.0,
#Probabilidades de transmisión por pareja
#Estos sigma_* representan la probabilidad de transmisión por pareja (no por acto individual), por año o por relación según cómo lo definas en tu fuerza de infección.
sigma_acute = 0.0082,
sigma_chronic = 0.0007,
sigma_aids = 0.0028,
#Reducción multiplicativa de infectividad por TAR
art_inf_red = 0.99,
#Conducta sexual y mezcla
partners_year = c(LR = 4.39, HR = 15.0), # nº parejas/año por grupo
condom_use = 0.17, # proporción de actos con condón
condom_eff = 0.95, # eficacia del condón
mixing_m = 0.3, # mezcla asortativa parcial LR/HR
#Reducción de conducta tras diagnóstico, TAR y SIDA
beh_red_dx = 0.20, # reducción tras diagnóstico VIH
beh_red_art = 0.50, # reducción adicional en TAR (supuesto)
beh_red_aids = 0.90, # reducción en SIDA (pacientes muy sintomáticos)
#Pruebas de VIH proporción de población que se realiza prueba en últimos 12 meses
tested_12m_prop = c(LR = 0.20, HR = 0.20),#proporción de cada grupo que reporta al menos una prueba de VIH en los últimos 12 meses (aquí 20% tanto en LR como en HR).
dur_recent_test_y = 1.0, #tiempo promedio que un susceptible permanece en el estado de “test negativo reciente” (S2)
omega2 = 1 / 1.0, #tasa de salida de S2 hacia S1 (pierde el “estatus de prueba reciente”).
#Tasas de detección por estadio (búsqueda pasiva)
psi4_LR = 0.10, # detección en agudo LR
psi4_HR = 0.10, # detección en agudo HR
psi7_LR = 0.10, # detección en crónico LR
psi7_HR = 0.10, # detección en crónico HR
psi10_LR = 1.00, # detección en SIDA LR
psi10_HR = 1.00, # detección en SIDA HR
#Tasas de prueba en susceptibles (se derivan de tested_12m_prop)
psi1_LR = NA_real_,
psi1_HR = NA_real_,
#TAR (tasas de inicio, independiente del estadio clínico)
alpha5 = 0.9506,
alpha8 = 0.9506,
alpha11 = 0.9506,
#Ahora bien: Es decir, en promedio ~1 año para entrar a TAR después de estar en A2, C2 o D2.
art_cov_target = 0.95, # objetivo de cobertura (parámetro auxiliar)
#PrEP
phi2 = 0.0, # entrada anual a PrEP (se modifica por escenarios)
phi3 = 0.0, # abandono de PrEP (0 -> 100% adherencia)
prep_cov_target = 0.0,#cobertura poblacional objetivo de PrEP (proporción de HSH susceptibles en PrEP) en cada escenario.
prep_uptake_m = 0.0,#parámetro auxiliar para modelar la velocidad / forma con la que se alcanza esa cobertura (por ejemplo, una curva logística). En la versión actual casi no lo usamos aún, pero queda declarado para poder implementar una trayectoria suave hacia la meta.
#En el escenario base se pondrán en 0 porque el protocolo dice que la cobertura de PrEP en Colombia es todavía muy baja / piloto, así que el status quo puede asumirse con PrEP casi inexistente.
prep_eff = 0.90, # eficacia biológica de PrEP
prep_adherence = 1.00, # adherencia promedio (1 = 100%)
#En el escenario base se asume que, dado que la PrEP casi no está implementada, quienes están en PrEP (escenarios futuros) tienen adherencia “ideal” (sirve como caso de referencia).En los escenarios de intervención se puede reducir la adherencia (por ejemplo, prep_adherence <- 0.8) o introducir abandono (phi3 > 0) para reflejar comportamientos más realistas, tal como se comenta en el protocolo cuando muestra cómo la eficacia cae con adherencia <80 %.
#Horizonte de simulación (40 años, paso 0.1)
t_start = 0,
t_end = 40,
dt = 0.1
)
#Completar tasas de pruebas basadas en la proporción tamizada
# Convertir proporción tamizada en 12 meses (p) a tasa constante psi_1
# usando: p = 1 - exp(-psi_1) => psi_1 = -log(1 - p)
parms_base$psi1_LR <- derive_psi_from_prop(parms_base$tested_12m_prop["LR"])
parms_base$psi1_HR <- derive_psi_from_prop(parms_base$tested_12m_prop["HR"])
#Completar tasas de progresión internas (mapeo a compartimentos 4–12)
# Se asume que el tiempo medio en el estadio no depende del diagnóstico ni del TAR:
# theta4 = theta5 = theta6 = theta_acute
# theta7 = theta8 = theta9 = theta_chronic
# theta10 = theta11 = theta12 = theta_aids
parms_base <- within(parms_base, {
theta4 <- theta_acute
theta5 <- theta_acute
theta6 <- theta_acute
theta7 <- theta_chronic
theta8 <- theta_chronic
theta9 <- theta_chronic
theta10 <- theta_aids
theta11 <- theta_aids
theta12 <- theta_aids
})6 5. Funciones auxiliares: tamaños, mezcla, infectividad y fuerza de infección:
Se implementan funciones que:
Calculan tamaño de los grupos de riesgo (LR y HR).
Construyen la matriz de mezcla sexual H, que respeta la reciprocidad de contactos.
Calculan la infectividad promedio por pareja en cada grupo, ponderando compartimentos (agudo, crónico, SIDA) y efectos de diagnóstico, TAR, SIDA y condón.
Derivan la fuerza de infección (λ) para susceptibles S1, S2 y S3, incluyendo la protección por PrEP.
#Tamaños de grupo LR / HR:
get_group_sizes <- function(state) {
# Seleccionar y sumar todos los compartimentos que terminan en "_LR"
N_LR <- sum(state[grepl("_LR$", names(state))], na.rm = TRUE)
# Seleccionar y sumar todos los compartimentos que terminan en "_HR"
N_HR <- sum(state[grepl("_HR$", names(state))], na.rm = TRUE)
# Devolver vector con tamaños por grupo de riesgo
c(LR = N_LR, HR = N_HR)
}
#Matriz de mezcla H entre LR y HR (mezcla asortativa parcial)
# N_LR, N_HR: tamaños de los grupos de bajo y alto riesgo
# m: índice de mezcla asortativa (0 = proporcional-Es decir, todos eligen parejas según la composición de la población-, 1 = totalmente asortativa-LR sólo con LR y HR sólo con HR)
mixing_matrix <- function(N_LR, N_HR, m) {
N_tot <- N_LR + N_HR
# Si el total no es válido, devolver mezcla 50/50 como valor de seguridad
if (!is.finite(N_tot) || N_tot <= 0) {
return(matrix(c(0.5, 0.5,
0.5, 0.5),
nrow = 2, byrow = TRUE,
dimnames = list(c("LR","HR"), c("LR","HR"))))
}
# Proporciones poblacionales de cada grupo
p_LR <- N_LR / N_tot
p_HR <- N_HR / N_tot
# Entradas de la matriz de mezcla:
# H11: P(pareja LR | individuo LR)
# H12: P(pareja HR | individuo LR)
# H21: P(pareja LR | individuo HR)
# H22: P(pareja HR | individuo HR)
H11 <- m + (1 - m) * p_LR
H12 <- (1 - m) * p_HR
H21 <- (1 - m) * p_LR
H22 <- m + (1 - m) * p_HR
# Construir la matriz 2x2 con nombres de filas y columnas
matrix(c(H11, H21,
H12, H22),
nrow = 2, byrow = FALSE,
dimnames = list(c("LR","HR"), c("LR","HR")))
}
#Infectividad promedio por pareja en cada grupo (LR / HR)
infectivity_by_group <- function(state, parms) {
with(as.list(c(parms, as.list(state))), {
# factor de protección por condón
condom_factor <- 1 - condom_use * condom_eff
if (!is.finite(condom_factor) || condom_factor < 0) condom_factor <- 0
# modificadores de conducta
f_dx <- 1 - beh_red_dx#reducción proporcional de conducta tras diagnóstico.
f_art <- 1 - beh_red_art#reducción adicional en TAR.
f_aids <- 1 - beh_red_aids#reducción en SIDA (por enfermedad avanzada).
#Nº total de infectados en cada grupo
I_LR <- sum(c(A1_LR, A2_LR, A3_LR,
C1_LR, C2_LR, C3_LR,
D1_LR, D2_LR, D3_LR), na.rm = TRUE)
I_HR <- sum(c(A1_HR, A2_HR, A3_HR,
C1_HR, C2_HR, C3_HR,
D1_HR, D2_HR, D3_HR), na.rm = TRUE)
#sigma_*: probabilidad básica de transmisión,
#condom_factor: efecto del condón,
#f_dx, f_art, f_aids: reducciones de conducta,
#(1 - art_inf_red) para los que están en TAR, porque el TAR reduce además la infectividad biológica.
# Infectividad promedio en LR
if (is.finite(I_LR) && I_LR > 0) {
LR_A1 <- sigma_acute * condom_factor * A1_LR
LR_A2 <- sigma_acute * condom_factor * f_dx * A2_LR
LR_A3 <- sigma_acute * condom_factor * f_dx * f_art *
(1 - art_inf_red) * A3_LR
LR_C1 <- sigma_chronic * condom_factor * C1_LR
LR_C2 <- sigma_chronic * condom_factor * f_dx * C2_LR
LR_C3 <- sigma_chronic * condom_factor * f_dx * f_art *
(1 - art_inf_red) * C3_LR
LR_D1 <- sigma_aids * condom_factor * f_aids * D1_LR
LR_D2 <- sigma_aids * condom_factor * f_dx * f_aids * D2_LR
LR_D3 <- sigma_aids * condom_factor * f_dx * f_art * f_aids *
(1 - art_inf_red) * D3_LR
# Promedio por individuo infectado LR
beta_LR <- (LR_A1 + LR_A2 + LR_A3 +
LR_C1 + LR_C2 + LR_C3 +
LR_D1 + LR_D2 + LR_D3) / I_LR
} else {
beta_LR <- 0
}
# Infectividad promedio en HR (análoga)
if (is.finite(I_HR) && I_HR > 0) {
HR_A1 <- sigma_acute * condom_factor * A1_HR
HR_A2 <- sigma_acute * condom_factor * f_dx * A2_HR
HR_A3 <- sigma_acute * condom_factor * f_dx * f_art *
(1 - art_inf_red) * A3_HR
HR_C1 <- sigma_chronic * condom_factor * C1_HR
HR_C2 <- sigma_chronic * condom_factor * f_dx * C2_HR
HR_C3 <- sigma_chronic * condom_factor * f_dx * f_art *
(1 - art_inf_red) * C3_HR
HR_D1 <- sigma_aids * condom_factor * f_aids * D1_HR
HR_D2 <- sigma_aids * condom_factor * f_dx * f_aids * D2_HR
HR_D3 <- sigma_aids * condom_factor * f_dx * f_art * f_aids *
(1 - art_inf_red) * D3_HR
# Promedio por individuo infectado HR
beta_HR <- (HR_A1 + HR_A2 + HR_A3 +
HR_C1 + HR_C2 + HR_C3 +
HR_D1 + HR_D2 + HR_D3) / I_HR
} else {
beta_HR <- 0
}
#Limpieza defensiva
if (!is.finite(beta_LR) || beta_LR < 0) beta_LR <- 0
if (!is.finite(beta_HR) || beta_HR < 0) beta_HR <- 0
c(LR = beta_LR, HR = beta_HR)
})
}
#Fuerza de infección (λ) para susceptibles S1, S2, S3
foi_susceptibles <- function(t, state, parms) {
with(as.list(c(parms, as.list(state))), {
# 1) Tamaño de los grupos de riesgo y matriz de mezcla sexual H
N_groups <- get_group_sizes(state)
H <- mixing_matrix(N_groups["LR"], N_groups["HR"], mixing_m)
# 2) Infectividad promedio por pareja en cada grupo (β_LR, β_HR)
beta <- infectivity_by_group(state, parms)
# 3) Número de parejas por año en cada grupo
c_LR <- partners_year["LR"]
c_HR <- partners_year["HR"]
# Fuerza de infección base en cada grupo(sin distinguir S1/S2/S3)
# λ_g = c_g * (H_{g,LR} * β_LR + H_{g,HR} * β_HR)
lambda_LR_base <- H["LR","LR"] * c_LR * beta["LR"] +
H["LR","HR"] * c_LR * beta["HR"]
lambda_HR_base <- H["HR","LR"] * c_HR * beta["LR"] +
H["HR","HR"] * c_HR * beta["HR"]
if (!is.finite(lambda_LR_base) || lambda_LR_base < 0) lambda_LR_base <- 0
if (!is.finite(lambda_HR_base) || lambda_HR_base < 0) lambda_HR_base <- 0
# 5) Factor de protección por PrEP:
# λ_PrEP = (1 - prep_eff * prep_adherence) * λ_base
prep_factor <- 1 - (prep_eff * prep_adherence)
if (!is.finite(prep_factor) || prep_factor < 0) prep_factor <- 0
if (prep_factor > 1) prep_factor <- 1
# Susceptibles sin PrEP (S1, S2): misma λ base
#λ para cada tipo de susceptible
# S1, S2: sin PrEP -> λ = λ_base
lambda1_LR <- lambda_LR_base
lambda2_LR <- lambda_LR_base
lambda1_HR <- lambda_HR_base
lambda2_HR <- lambda_HR_base
# Susceptibles en PrEP (S3): λ reducida
# S3: con PrEP -> λ = prep_factor * λ_base
lambda3_LR <- lambda_LR_base * prep_factor
lambda3_HR <- lambda_HR_base * prep_factor
lambdas <- list(
lambda1_LR = lambda1_LR,
lambda2_LR = lambda2_LR,
lambda3_LR = lambda3_LR,
lambda1_HR = lambda1_HR,
lambda2_HR = lambda2_HR,
lambda3_HR = lambda3_HR
)
# Limpieza defensiva de valores no finitos o negativos
lambdas <- lapply(lambdas, function(x) {
x[!is.finite(x) | x < 0] <- 0
x
})
lambdas
})
}7 6. Sistema de ecuaciones diferenciales (EDO):
#“El sistema de EDO que usamos es el límite determinístico de un proceso de Markov en tiempo continuo, donde cada transición (infección, progresión, diagnóstico, TAR, muerte) ocurre con una tasa constante. A nivel individual estos tiempos son exponenciales; a nivel poblacional, la ley de los grandes números nos permite describir la dinámica mediante ecuaciones diferenciales del tipo ‘entradas menos salidas’, donde las salidas son la suma de hazards multiplicada por el tamaño del compartimento.”
ode_vi_hsh <- function(t, state, parms) {
#Esta es la función que recibe deSolve::ode()
#t: tiempo actual (lo pasa ode()).
#state: vector con los 24 estados (S1_LR, S2_HR, A1_LR, etc.).
#parms: lista de parámetros (parms_base).
#Fuerza de infección en susceptibles
lambdas <- foi_susceptibles(t, state, parms)
#Tamaños de grupo LR y HR
N_LR <- sum(state[grepl("_LR$", names(state))], na.rm = TRUE)
N_HR <- sum(state[grepl("_HR$", names(state))], na.rm = TRUE)
with(as.list(c(parms,
as.list(state),
lambdas,
list(N_LR = N_LR, N_HR = N_HR))), {
d_state <- setNames(numeric(length(state)), names(state))
#d_state es el vector donde vas a guardar dX/dt para cada compartimento.
## Bloque 1 – Susceptibles
# LR
d_state["S1_LR"] <-
rho_LR * N_LR +
omega2 * S2_LR +
phi3 * S3_LR -
lambda1_LR * S1_LR -
psi1_LR * S1_LR -
mu * S1_LR
d_state["S2_LR"] <-
psi1_LR * S1_LR -
lambda2_LR * S2_LR -
(omega2 + phi2 + mu) * S2_LR
d_state["S3_LR"] <-
phi2 * S2_LR -
lambda3_LR * S3_LR -
(phi3 + mu) * S3_LR
# HR
d_state["S1_HR"] <-
rho_HR * N_HR +
omega2 * S2_HR +
phi3 * S3_HR -
lambda1_HR * S1_HR -
psi1_HR * S1_HR -
mu * S1_HR
d_state["S2_HR"] <-
psi1_HR * S1_HR -
lambda2_HR * S2_HR -
(omega2 + phi2 + mu) * S2_HR
d_state["S3_HR"] <-
phi2 * S2_HR -
lambda3_HR * S3_HR -
(phi3 + mu) * S3_HR
## Bloque 2 – Infección aguda (A1, A2, A3)
# LR
d_state["A1_LR"] <-
lambda1_LR * S1_LR +
lambda2_LR * S2_LR -
(psi4_LR + mu + theta4 + b4) * A1_LR
d_state["A2_LR"] <-
lambda3_LR * S3_LR +
psi4_LR * A1_LR -
(mu + alpha5 + theta5 + b5) * A2_LR
d_state["A3_LR"] <-
alpha5 * A2_LR -
(mu + theta6 + b6) * A3_LR
# HR
d_state["A1_HR"] <-
lambda1_HR * S1_HR +
lambda2_HR * S2_HR -
(psi4_HR + mu + theta4 + b4) * A1_HR
d_state["A2_HR"] <-
lambda3_HR * S3_HR +
psi4_HR * A1_HR -
(mu + alpha5 + theta5 + b5) * A2_HR
d_state["A3_HR"] <-
alpha5 * A2_HR -
(mu + theta6 + b6) * A3_HR
## Bloque 3 – Infección crónica (C1, C2, C3)
# LR
d_state["C1_LR"] <-
theta4 * A1_LR -
(psi7_LR + mu + theta7 + b7) * C1_LR
d_state["C2_LR"] <-
theta5 * A2_LR +
psi7_LR * C1_LR -
(mu + alpha8 + theta8 + b8) * C2_LR
d_state["C3_LR"] <-
theta6 * A3_LR +
alpha8 * C2_LR -
(mu + theta9 + b9) * C3_LR
# HR
d_state["C1_HR"] <-
theta4 * A1_HR -
(psi7_HR + mu + theta7 + b7) * C1_HR
d_state["C2_HR"] <-
theta5 * A2_HR +
psi7_HR * C1_HR -
(mu + alpha8 + theta8 + b8) * C2_HR
d_state["C3_HR"] <-
theta6 * A3_HR +
alpha8 * C2_HR -
(mu + theta9 + b9) * C3_HR
## Bloque 4 – SIDA (D1, D2, D3)
# LR
d_state["D1_LR"] <-
theta7 * C1_LR -
(psi10_LR + mu + theta10 + b10) * D1_LR
d_state["D2_LR"] <-
theta8 * C2_LR +
psi10_LR * D1_LR -
alpha11 * D2_LR -
(mu + b11) * D2_LR
d_state["D3_LR"] <-
theta9 * C3_LR +
alpha11 * D2_LR -
(mu + theta12 + b12) * D3_LR
# HR
d_state["D1_HR"] <-
theta7 * C1_HR -
(psi10_HR + mu + theta10 + b10) * D1_HR
d_state["D2_HR"] <-
theta8 * C2_HR +
psi10_HR * D1_HR -
alpha11 * D2_HR -
(mu + b11) * D2_HR
d_state["D3_HR"] <-
theta9 * C3_HR +
alpha11 * D2_HR -
(mu + theta12 + b12) * D3_HR
list(d_state)
})
}8 7. Función de simulación:
sim_vi_hsh <- function(parms = parms_base,
init = init,
t_start = parms$t_start,
t_end = parms$t_end,
dt = parms$dt) {
times <- seq(t_start, t_end, by = dt)
out <- deSolve::ode(
y = init,
times = times,
func = ode_vi_hsh,
parms = parms
)
out <- as.data.frame(out)
lr_cols <- grepl("_LR$", names(out))
hr_cols <- grepl("_HR$", names(out))
out$N_LR <- rowSums(out[, lr_cols, drop = FALSE], na.rm = TRUE)
out$N_HR <- rowSums(out[, hr_cols, drop = FALSE], na.rm = TRUE)
out$N_tot <- out$N_LR + out$N_HR
inf_LR_idx <- grepl("^(A|C|D)[123]_LR$", names(out))
inf_HR_idx <- grepl("^(A|C|D)[123]_HR$", names(out))
out$I_LR <- rowSums(out[, inf_LR_idx, drop = FALSE], na.rm = TRUE)
out$I_HR <- rowSums(out[, inf_HR_idx, drop = FALSE], na.rm = TRUE)
out$I_tot <- out$I_LR + out$I_HR
out$prev_LR <- ifelse(out$N_LR > 0, out$I_LR / out$N_LR, NA_real_)
out$prev_HR <- ifelse(out$N_HR > 0, out$I_HR / out$N_HR, NA_real_)
out$prev_tot <- ifelse(out$N_tot > 0, out$I_tot / out$N_tot, NA_real_)
tibble::as_tibble(out)
}
#Ejecución base con parámetros iniciales
res_base <- sim_vi_hsh(parms = parms_base, init = init)
head(res_base)## # A tibble: 6 × 34
## time S1_LR S2_LR S3_LR A1_LR A2_LR A3_LR C1_LR C2_LR C3_LR D1_LR D2_LR
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 274625. 69525. 3476. 763. 572. 2480. 11446. 8585. 37200. 3052. 2289.
## 2 0.1 275227. 68706. 3474. 732. 485. 2325. 11241. 7850. 37652. 2342. 2084.
## 3 0.2 275739. 67977. 3473. 701. 412. 2175. 11038. 7182. 38019. 1824. 1865.
## 4 0.3 276173. 67328. 3471. 672. 350. 2031. 10837. 6576. 38308. 1444. 1651.
## 5 0.4 276538. 66750. 3469. 645. 299. 1894. 10638. 6026. 38527. 1166. 1452.
## 6 0.5 276842. 66234. 3467. 618. 255. 1763. 10442. 5527. 38682. 962. 1272.
## # ℹ 22 more variables: D3_LR <dbl>, S1_HR <dbl>, S2_HR <dbl>, S3_HR <dbl>,
## # A1_HR <dbl>, A2_HR <dbl>, A3_HR <dbl>, C1_HR <dbl>, C2_HR <dbl>,
## # C3_HR <dbl>, D1_HR <dbl>, D2_HR <dbl>, D3_HR <dbl>, N_LR <dbl>, N_HR <dbl>,
## # N_tot <dbl>, I_LR <dbl>, I_HR <dbl>, I_tot <dbl>, prev_LR <dbl>,
## # prev_HR <dbl>, prev_tot <dbl>ggplot(res_base, aes(x = time, y = prev_tot)) +
geom_line() +
labs(x = "Tiempo (años)", y = "Prevalencia total VIH HSH (status quo)") +
theme_minimal()
library(dplyr)
library(ggplot2)
library(purrr)
library(tidyr)
# 1. Definir los niveles de rho a comparar
mu_val <- parms_base$mu
rho_levels <- c(
0, # cohorte cerrada
mu_val, # rho = mu (aprox. estacionario sin VIH)
1.5 * mu_val, # rho > mu
2.0 * mu_val # rho aún mayor
)
# 2. Simular para cada valor de rho, con rho_LR = rho_HR = mismo valor
sim_rho <- purrr::map_dfr(rho_levels, function(rho) {
# Copia de parámetros para no alterar parms_base global
parms_rho <- parms_base
parms_rho$rho_LR <- rho
parms_rho$rho_HR <- rho
res <- sim_vi_hsh(parms = parms_rho, init = init)
res |>
mutate(rho = rho)
})
# 3. (Opcional) Normalizar N_tot por el valor inicial para ver proporción
N0 <- sim_rho |> filter(time == min(time)) |> slice(1) |> pull(N_tot)
sim_rho <- sim_rho |>
mutate(N_tot_rel = N_tot / N0)
# 4. Graficar tamaño total de la población HSH en el tiempo por escenario de rho
ggplot(sim_rho, aes(x = time, y = N_tot_rel,
colour = factor(round(rho, 5)))) +
geom_line(linewidth = 1) +
labs(
x = "Tiempo (años)",
y = "Tamaño relativo de la cohorte HSH (N_tot / N0)",
colour = expression(rho),
title = "Evolución del tamaño de la cohorte HSH según tasa de entrada demográfica (rho)",
subtitle = "rho_LR = rho_HR; valores comparados: 0, mu, 1.5*mu, 2*mu"
) +
theme_minimal() +
theme(
legend.position = "bottom"
)
library(dplyr)
library(ggplot2)
library(purrr)
library(tidyr)
## 1. Definir niveles de rho
mu_val <- parms_base$mu
rho_levels <- c(
0, # cohorte cerrada
mu_val, # rho = mu
1.5 * mu_val, # rho > mu
2.0 * mu_val # rho aún mayor
)
## 2. Simular para cada rho
sim_rho <- purrr::map_dfr(rho_levels, function(rho) {
parms_rho <- parms_base
parms_rho$rho_LR <- rho
parms_rho$rho_HR <- rho
res <- sim_vi_hsh(parms = parms_rho, init = init)
res |>
mutate(rho = rho)
})
## 3. Normalizar N_tot y dejar prevalencia
N0 <- sim_rho |>
filter(time == min(time)) |>
slice(1) |>
pull(N_tot)
sim_rho <- sim_rho |>
mutate(
N_tot_rel = N_tot / N0, # tamaño relativo cohorte
prev_tot = prev_tot # ya viene de sim_vi_hsh
)
## 4. Pasar a formato largo para graficar N_tot_rel y prev_tot en facetas
sim_long <- sim_rho |>
select(time, rho, N_tot_rel, prev_tot) |>
tidyr::pivot_longer(
cols = c(N_tot_rel, prev_tot),
names_to = "metric",
values_to = "value"
) |>
mutate(
metric = dplyr::recode(
metric,
N_tot_rel = "Tamaño relativo cohorte (N_tot / N0)",
prev_tot = "Prevalencia VIH (prev_tot)"
)
)
## 5. Graficar: dos paneles, una línea por rho
ggplot(sim_long,
aes(x = time,
y = value,
colour = factor(round(rho, 5)))) +
geom_line(linewidth = 1) +
facet_wrap(~ metric, scales = "free_y", ncol = 1) +
labs(
x = "Tiempo (años)",
y = NULL,
colour = expression(rho),
title = "Dinámica de la cohorte HSH y prevalencia de VIH según tasa de entrada demográfica (rho)",
subtitle = "rho_LR = rho_HR; se comparan valores 0, mu, 1.5*mu y 2*mu"
) +
theme_minimal() +
theme(
legend.position = "bottom"
)
9 8. Funciones de indicadores (incidencia sin calibrar)
Se cuantifica:
Nuevas infecciones por unidad de tiempo (new_inf_year).
Incidencia instantánea por persona-año (inc_rate_inst) y por 1.000 persona-año.
# Añade indicadores de incidencia a una simulación (sim_vi_hsh)
#La función NO vuelve a resolver el modelo.
#Toma una simulación ya hecha y le añade:
#número de nuevas infecciones por unidad de tiempo (new_inf_year),
#tasa de incidencia instantánea (inc_rate_inst),
#incidencia por 1000 persona-año (inc_per1000_inst).
add_indicators_single <- function(sim, parms) {
times <- sim$time
n <- nrow(sim)
# dt entre puntos de tiempo (se repite el último para igualar longitud)
if (n > 1) {
dt_vec <- c(diff(times), tail(diff(times), 1))
} else {
dt_vec <- 0
}
# Vectores para nuevas infecciones e incidencia instantánea
new_inf_year <- numeric(n)
inc_rate_inst <- numeric(n)
for (i in seq_len(n)) {
# Vector de estados en el tiempo i
state_vec <- as.numeric(sim[i, state_names])
names(state_vec) <- state_names
# Fuerzas de infección λ1, λ2, λ3 en LR/HR
lambdas <- foi_susceptibles(
t = times[i],
state = state_vec,
parms = parms
)
# Calcular nuevas infecciones y tasa de incidencia instantánea
tmp <- with(as.list(c(parms, as.list(state_vec), lambdas)), {
# Nuevas infecciones en el instante t_i:
# suma de λ_k,g * S_k,g sobre k = 1,2,3; g = LR, HR
new_inf <- lambda1_LR * S1_LR + lambda2_LR * S2_LR + lambda3_LR * S3_LR +
lambda1_HR * S1_HR + lambda2_HR * S2_HR + lambda3_HR * S3_HR
# Incidencia instantánea por persona-año
inc_rate <- ifelse(sim$N_tot[i] > 0, new_inf / sim$N_tot[i], NA_real_)
list(new_inf = new_inf, inc_rate = inc_rate)
})
new_inf_year[i] <- tmp$new_inf
inc_rate_inst[i] <- tmp$inc_rate
}
# Devolver la simulación con columnas adicionales
sim %>%
dplyr::mutate(
dt = dt_vec,
new_inf_year = new_inf_year,
inc_rate_inst = inc_rate_inst,
inc_per1000_inst = inc_rate_inst * 1000
)
}10 9. Calibración de la incidencia inicial:
El objetivo de esta sección es ajustar un subconjunto de parámetros de transmisión de forma que la incidencia media del primer año en el escenario Status quo sea consistente con la evidencia disponible para HSH en Colombia.
Según el estudio de Zea et al. (2017), la incidencia anual en HSH en Bogotá se estimó en aproximadamente 0,25 casos por 100 HSH-año, es decir:
\[ I_{\text{obs}} = 0.25 \; \text{casos / 100 HSH-año} = 2.5 \; \text{casos / 1000 HSH-año} = 0.0025 \; \text{casos / HSH-año}. \]
En el modelo, la incidencia instantánea se define como
\[ \text{inc}(t) = \frac{\text{new\_inf}(t)}{N_{\text{tot}}(t)}, \]
y la incidencia media durante el primer año (año 1 del horizonte de simulación) para un conjunto de parámetros \(\theta\) se aproxima por
\[ I_1(\theta) \approx \frac{ \displaystyle \sum_{t_i \in [0,1)} \text{new\_inf}(t_i;\theta)\,\Delta t_i }{ \displaystyle \sum_{t_i \in [0,1)} N_{\text{tot}}(t_i;\theta)\,\Delta t_i }. \]
Buscamos un factor de escala \(\kappa\) que multiplique las probabilidades de transmisión por pareja
\[ \sigma_{\text{agudo}},\; \sigma_{\text{crónico}},\; \sigma_{\text{SIDA}} \]
de forma que la incidencia media en el primer año se acerque al valor observado \(I_{\text{obs}} = 0.0025\). Es decir, definimos
\[ \sigma_{\text{agudo}}^{(\kappa)} = \kappa \,\sigma_{\text{agudo}}^{(0)},\quad \sigma_{\text{crónico}}^{(\kappa)} = \kappa \,\sigma_{\text{crónico}}^{(0)},\quad \sigma_{\text{SIDA}}^{(\kappa)} = \kappa \,\sigma_{\text{SIDA}}^{(0)}, \]
y buscamos
\[ \kappa^\ast = \arg\min_{\kappa > 0} \left\{ \left[ I_1\bigl(\theta(\kappa)\bigr) - I_{\text{obs}} \right]^2 \right\}. \]
En la práctica, este problema se resuelve numéricamente mediante una búsqueda unidimensional sobre \(\kappa\).
# Función que calcula la incidencia media del primer año
# para una lista de parámetros 'parms' y condiciones iniciales 'init'.
compute_inc_year1 <- function(parms, init) {
sim <- sim_vi_hsh(parms = parms, init = init)
sim_i <- add_indicators_single(sim, parms)
year1 <- sim_i$time >= 0 & sim_i$time < 1
cases_year1 <- sum(
sim_i$new_inf_year[year1] * sim_i$dt[year1],
na.rm = TRUE
)
pya_year1 <- sum(
sim_i$N_tot[year1] * sim_i$dt[year1],
na.rm = TRUE
)
inc1 <- cases_year1 / pya_year1 # casos por persona-año
inc1
}
# Incidencia objetivo según Zea et al. (2017):
# 0.25 casos por 100 HSH-año = 2.5 por 1000 HSH-año = 0.0025 por HSH-año
inc_target <- 0.0025 # casos por persona-año
# --- 9.1 Función objetivo para buscar kappa (escala común de sigma_*) ----
obj_fun <- function(kappa, parms_base, init, inc_target) {
# 1) Clonar parámetros base
p <- parms_base
# 2) Escalar probabilidades de transmisión por pareja
# (se aplica el mismo factor kappa a agudo, crónico y SIDA)
p$sigma_acute <- kappa * parms_base$sigma_acute
p$sigma_chronic <- kappa * parms_base$sigma_chronic
p$sigma_aids <- kappa * parms_base$sigma_aids
# 3) Incidencia media en el primer año con estos parámetros
inc1_k <- compute_inc_year1(p, init)
# 4) Diferencia respecto al objetivo (esto es lo que uniroot hace = 0)
inc1_k - inc_target
}
# --- 9.2 Búsqueda numérica de kappa* usando uniroot -----------------------
# Intervalo razonable para kappa (ajustable si es necesario)
kappa_lower <- 0.1
kappa_upper <- 5
# (Opcional) Diagnóstico: incidencia con los extremos del intervalo
inc_lower <- compute_inc_year1(
parms = within(parms_base, {
sigma_acute <- kappa_lower * sigma_acute
sigma_chronic <- kappa_lower * sigma_chronic
sigma_aids <- kappa_lower * sigma_aids
}),
init = init
)
inc_upper <- compute_inc_year1(
parms = within(parms_base, {
sigma_acute <- kappa_upper * sigma_acute
sigma_chronic <- kappa_upper * sigma_chronic
sigma_aids <- kappa_upper * sigma_aids
}),
init = init
)
cat("Incidencia año 1 con kappa =", kappa_lower, ":",
inc_lower * 1000, "por 1000 p-año\n")## Incidencia año 1 con kappa = 0.1 : 0.106057 por 1000 p-añocat("Incidencia año 1 con kappa =", kappa_upper, ":",
inc_upper * 1000, "por 1000 p-año\n")## Incidencia año 1 con kappa = 5 : 7.398172 por 1000 p-año# Resolver kappa* tal que inc1(kappa) = inc_target
kappa_hat <- tryCatch(
{
uniroot(
f = obj_fun,
interval = c(kappa_lower, kappa_upper),
parms_base = parms_base,
init = init,
inc_target = inc_target
)$root
},
error = function(e) {
warning(
"No se pudo encontrar kappa en el intervalo dado: ",
conditionMessage(e)
)
NA_real_
}
)
kappa_hat## [1] 2.069015# --- 9.3 Actualizar parámetros base con el kappa calibrado ---------------
parms_calibrated <- parms_base
if (!is.na(kappa_hat)) {
parms_calibrated$sigma_acute <- kappa_hat * parms_base$sigma_acute
parms_calibrated$sigma_chronic <- kappa_hat * parms_base$sigma_chronic
parms_calibrated$sigma_aids <- kappa_hat * parms_base$sigma_aids
}
# Incidencia resultante con los parámetros calibrados
inc1_cal <- compute_inc_year1(parms_calibrated, init)
cat("Incidencia objetivo (Zea et al.):",
inc_target * 1000, "casos / 1000 p-año\n")## Incidencia objetivo (Zea et al.): 2.5 casos / 1000 p-añocat("Incidencia calibrada año 1:",
inc1_cal * 1000, "casos / 1000 p-año\n")## Incidencia calibrada año 1: 2.5 casos / 1000 p-añocat("Factor de escala kappa_hat:",
kappa_hat, "\n")## Factor de escala kappa_hat: 2.069015# OJO: a partir de aquí usamos siempre parms_base calibrado
parms_base <- parms_calibrated10.1 Volver a correr todo lo anterior: Ya calibrado
sim_vi_hsh <- function(parms = parms_base,
init = init,
t_start = parms$t_start,
t_end = parms$t_end,
dt = parms$dt) {
times <- seq(t_start, t_end, by = dt)
out <- deSolve::ode(
y = init,
times = times,
func = ode_vi_hsh,
parms = parms
)
out <- as.data.frame(out)
lr_cols <- grepl("_LR$", names(out))
hr_cols <- grepl("_HR$", names(out))
out$N_LR <- rowSums(out[, lr_cols, drop = FALSE], na.rm = TRUE)
out$N_HR <- rowSums(out[, hr_cols, drop = FALSE], na.rm = TRUE)
out$N_tot <- out$N_LR + out$N_HR
inf_LR_idx <- grepl("^(A|C|D)[123]_LR$", names(out))
inf_HR_idx <- grepl("^(A|C|D)[123]_HR$", names(out))
out$I_LR <- rowSums(out[, inf_LR_idx, drop = FALSE], na.rm = TRUE)
out$I_HR <- rowSums(out[, inf_HR_idx, drop = FALSE], na.rm = TRUE)
out$I_tot <- out$I_LR + out$I_HR
out$prev_LR <- ifelse(out$N_LR > 0, out$I_LR / out$N_LR, NA_real_)
out$prev_HR <- ifelse(out$N_HR > 0, out$I_HR / out$N_HR, NA_real_)
out$prev_tot <- ifelse(out$N_tot > 0, out$I_tot / out$N_tot, NA_real_)
tibble::as_tibble(out)
}
#Ejecución base con parámetros iniciales
res_base <- sim_vi_hsh(parms = parms_base, init = init)
head(res_base)## # A tibble: 6 × 34
## time S1_LR S2_LR S3_LR A1_LR A2_LR A3_LR C1_LR C2_LR C3_LR D1_LR D2_LR
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 274625. 69525. 3476. 763. 572. 2480. 11446. 8585. 37200. 3052. 2289.
## 2 0.1 275192. 68697. 3474. 773. 485. 2325. 11243. 7850. 37652. 2342. 2084.
## 3 0.2 275671. 67960. 3472. 779. 412. 2175. 11044. 7182. 38019. 1824. 1865.
## 4 0.3 276072. 67303. 3471. 782. 352. 2031. 10851. 6577. 38308. 1444. 1651.
## 5 0.4 276404. 66717. 3469. 784. 301. 1894. 10662. 6027. 38527. 1166. 1452.
## 6 0.5 276674. 66194. 3467. 783. 259. 1764. 10478. 5528. 38683. 962. 1272.
## # ℹ 22 more variables: D3_LR <dbl>, S1_HR <dbl>, S2_HR <dbl>, S3_HR <dbl>,
## # A1_HR <dbl>, A2_HR <dbl>, A3_HR <dbl>, C1_HR <dbl>, C2_HR <dbl>,
## # C3_HR <dbl>, D1_HR <dbl>, D2_HR <dbl>, D3_HR <dbl>, N_LR <dbl>, N_HR <dbl>,
## # N_tot <dbl>, I_LR <dbl>, I_HR <dbl>, I_tot <dbl>, prev_LR <dbl>,
## # prev_HR <dbl>, prev_tot <dbl>ggplot(res_base, aes(x = time, y = prev_tot)) +
geom_line() +
labs(x = "Tiempo (años)", y = "Prevalencia total VIH HSH (status quo)") +
theme_minimal()
10.2 Efecto del tamaño ya calibrado
library(dplyr)
library(ggplot2)
library(purrr)
library(tidyr)
# 1. Definir los niveles de rho a comparar
mu_val <- parms_base$mu
rho_levels <- c(
0, # cohorte cerrada
mu_val, # rho = mu (aprox. estacionario sin VIH)
1.5 * mu_val, # rho > mu
2.0 * mu_val # rho aún mayor
)
# 2. Simular para cada valor de rho, con rho_LR = rho_HR = mismo valor
sim_rho <- purrr::map_dfr(rho_levels, function(rho) {
# Copia de parámetros para no alterar parms_base global
parms_rho <- parms_base
parms_rho$rho_LR <- rho
parms_rho$rho_HR <- rho
res <- sim_vi_hsh(parms = parms_rho, init = init)
res |>
mutate(rho = rho)
})
# 3. (Opcional) Normalizar N_tot por el valor inicial para ver proporción
N0 <- sim_rho |> filter(time == min(time)) |> slice(1) |> pull(N_tot)
sim_rho <- sim_rho |>
mutate(N_tot_rel = N_tot / N0)
# 4. Graficar tamaño total de la población HSH en el tiempo por escenario de rho
ggplot(sim_rho, aes(x = time, y = N_tot_rel,
colour = factor(round(rho, 5)))) +
geom_line(linewidth = 1) +
labs(
x = "Tiempo (años)",
y = "Tamaño relativo de la cohorte HSH (N_tot / N0)",
colour = expression(rho),
title = "Evolución del tamaño de la cohorte HSH según tasa de entrada demográfica (rho)",
subtitle = "rho_LR = rho_HR; valores comparados: 0, mu, 1.5*mu, 2*mu"
) +
theme_minimal() +
theme(
legend.position = "bottom"
)
10.3 Efecto en la prevalencia calibrado:
library(dplyr)
library(purrr)
library(ggplot2)
# 1. Definir niveles de rho: iguales o mayores que mu
mu_val <- parms_base$mu
rho_levels <- c(
mu_val, # rho = mu
1.5 * mu_val, # rho = 1.5 * mu
2.0 * mu_val # rho = 2 * mu
)
# 2. Simular para cada valor de rho (rho_LR = rho_HR = rho)
sim_rho_prev <- map_dfr(rho_levels, function(rho) {
parms_rho <- parms_base
parms_rho$rho_LR <- rho
parms_rho$rho_HR <- rho
sim_vi_hsh(parms = parms_rho, init = init) |>
mutate(rho = rho)
})
# 3. Graficar sólo la prevalencia total VIH HSH
ggplot(sim_rho_prev,
aes(x = time,
y = prev_tot,
colour = factor(round(rho, 5)))) +
geom_line(linewidth = 1) +
labs(
x = "Tiempo (años)",
y = "Prevalencia total VIH HSH",
colour = expression(rho),
title = "Impacto de distintos valores de rho en la prevalencia de VIH en HSH",
subtitle = "rho_LR = rho_HR; se muestran escenarios con rho = mu, 1.5*mu y 2*mu"
) +
theme_minimal() +
theme(
legend.position = "bottom"
)
10.4 Efecto en la tasa de incidencia calibrado:
library(dplyr)
library(purrr)
library(ggplot2)
library(tidyr)
library(knitr)
## 1. Definir niveles de rho (incluyo 0 como referencia + >= mu)
mu_val <- parms_base$mu
rho_levels <- c(
0, # cohorte cerrada
mu_val, # rho = mu
1.5 * mu_val, # rho = 1.5 * mu
2.0 * mu_val # rho = 2 * mu
)
rho_levels## [1] 0.0000 0.0052 0.0078 0.0104# sólo para ver los valores numéricos
## 2. Simular para cada rho y añadir indicadores de incidencia
sim_rho_inc <- map_dfr(rho_levels, function(rho) {
# copiar parámetros y modificar rho_LR / rho_HR
parms_rho <- parms_base
parms_rho$rho_LR <- rho
parms_rho$rho_HR <- rho
# correr modelo
sim <- sim_vi_hsh(parms = parms_rho, init = init)
# añadir incidencia
sim_ind <- add_indicators_single(sim, parms_rho)
# añadir columna con el rho usado
sim_ind |>
mutate(rho = rho)
})
## 3. Gráfico de incidencia instantánea por 1000 persona-año
ggplot(sim_rho_inc,
aes(x = time,
y = inc_per1000_inst,
colour = factor(round(rho, 5)))) +
geom_line(linewidth = 1) +
labs(
x = "Tiempo (años)",
y = "Incidencia instantánea (casos por 1.000 persona-año)",
colour = expression(rho),
title = "Impacto de distintos valores de rho en la incidencia de VIH en HSH",
subtitle = "rho_LR = rho_HR; se comparan 0, mu, 1.5*mu y 2*mu"
) +
theme_minimal() +
theme(
legend.position = "bottom"
)
## 4. Tabla de resumen / “tabla de convención” para evaluar diferencias
# función auxiliar para tomar el valor más cercano a un tiempo objetivo
closest_at <- function(x, time, target_time) {
x[which.min(abs(time - target_time))]
}
tabla_inc <- sim_rho_inc |>
group_by(rho) |>
summarise(
inc_ini = closest_at(inc_per1000_inst, time, 0), # incidencia en t ≈ 0
inc_5 = closest_at(inc_per1000_inst, time, 5), # a 5 años
inc_20 = closest_at(inc_per1000_inst, time, 20), # a 20 años
inc_40 = closest_at(inc_per1000_inst, time, 40), # a 40 años
inc_media = mean(inc_per1000_inst, na.rm = TRUE) # media en todo el horizonte
) |>
ungroup() |>
mutate(rho_label = round(rho, 5)) |>
select(rho_label, everything(), -rho)
# Mostrar la tabla bonita
kable(
tabla_inc,
digits = 3,
caption = "Incidencia de VIH (casos por 1.000 persona-año) según diferentes valores de rho (rho_LR = rho_HR)"
)## Warning in attr(x, "align"): 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")## Warning in attr(x, "format"): 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")| rho_label | inc_ini | inc_5 | inc_20 | inc_40 | inc_media |
|---|---|---|---|---|---|
| 0.000 | 2.737 | 2.516 | 12.040 | 9.513 | 8.529 |
| 0.005 | 2.737 | 2.598 | 12.193 | 9.934 | 8.906 |
| 0.008 | 2.737 | 2.640 | 12.266 | 10.138 | 9.094 |
| 0.010 | 2.737 | 2.683 | 12.338 | 10.338 | 9.280 |
11 10. Escenarios de intervención:
Se codifican escenarios de intervención que se comparan con el Status quo. Cada escenario ajusta parámetros clave:
Tasas de prueba y TAR (TasP).
Cobertura y adherencia a PrEP.
Uso de condón, número de parejas y conducta sexual.
library(tibble)
library(dplyr)
tabla_parametros_escenarios <- tibble(
parametro = c(
"tested_12m_prop",
"psi1_LR/psi1_HR",
"psi4/psi7/psi10",
"alpha5/alpha8/alpha11",
"art_cov_target",
"prep_cov_target",
"phi2",
"phi3",
"prep_adherence",
"condom_use",
"partners_year_HR",
"beh_red_dx",
"beh_red_art",
"beh_red_aids"
),
interpretacion = c(
"Proporción tamizada en 12 meses (LR, HR)",
"Tasa anual de test en S1 (LR, HR)",
"Tasas de diagnóstico en infectados (agudo, crónico, SIDA)",
"Tasas de inicio de TAR (agudo, crónico, SIDA)",
"Meta de cobertura de TAR",
"Cobertura objetivo de PrEP en susceptibles",
"Tasa de entrada a PrEP",
"Tasa de abandono de PrEP",
"Adherencia media a PrEP",
"Uso de condón en actos sexuales",
"Número de parejas/año en HR",
"Reducción de conducta tras diagnóstico",
"Reducción de conducta en TAR",
"Reducción de conducta en SIDA"
),
status_quo = c(
"LR = 0.50, HR = 0.60 (base_parms)",
"Derivada de tested_12m_prop base",
"Valores base (psi4_LR/HR, psi7_LR/HR, psi10_LR/HR)",
"Valores base (alpha5, alpha8, alpha11)",
"Valor base (p.ej. ≈ 0.90)",
"Muy baja / ~0 según base_parms",
"Valor base (p.ej. 0.10)",
"Valor base (p.ej. 0.20)",
"Valor base (p.ej. 0.60 o 1.00)",
"Valor base (p.ej. 0.35)",
"Valor base (p.ej. 14.1)",
"Valor base (p.ej. 0.30)",
"Valor base (p.ej. 0.50)",
"Valor base (definido en parms_base)"
),
tasp_80 = c(
"LR = 0.80, HR = 0.80",
"Derivada de 0.80 / 0.80",
"psi4/7/10 × 1.5",
"alpha5/8/11 × 1.5",
"0.95",
"= Status quo",
"= Status quo",
"= Status quo",
"= Status quo",
"= Status quo",
"= Status quo",
"= Status quo",
"= Status quo",
"= Status quo"
),
prep_20 = c(
"= Status quo",
"= Status quo",
"= Status quo",
"= Status quo",
"= Status quo",
"0.20",
"0.25",
"0.0 (abandono casi nulo)",
"min(0.85, base × 0.95)",
"= Status quo",
"= Status quo",
"= Status quo",
"= Status quo",
"= Status quo"
),
combinado_moderado = c(
"LR = 0.70, HR = 0.75",
"Derivada de 0.70 / 0.75",
"psi4/7/10 × 1.3",
"alpha5/8/11 × 1.3",
"0.93",
"0.20",
"0.20",
"0.0",
"min(0.80, base × 0.90)",
"base + 0.15 (tope 0.80)",
"= Status quo",
"= Status quo",
"= Status quo",
"= Status quo"
),
combinado_fuerte = c(
"LR = 0.90, HR = 0.90",
"Derivada de 0.90 / 0.90",
"psi4/7/10 × 2.0",
"alpha5/8/11 × 1.8",
"0.98",
"0.30",
"0.35",
"0.0",
"min(0.90, base × 0.95)",
"base + 0.25 (tope 0.95)",
"base × 0.7 (−30 % parejas)",
"base + 0.20 (tope 0.50)",
"base + 0.20 (tope 0.70)",
"base + 0.05 (tope 0.95)"
)
)
tabla_parametros_escenarios## # A tibble: 14 × 7
## parametro interpretacion status_quo tasp_80 prep_20 combinado_moderado
## <chr> <chr> <chr> <chr> <chr> <chr>
## 1 tested_12m_prop Proporción ta… LR = 0.50… LR = 0… = Stat… LR = 0.70, HR = 0…
## 2 psi1_LR/psi1_HR Tasa anual de… Derivada … Deriva… = Stat… Derivada de 0.70 …
## 3 psi4/psi7/psi10 Tasas de diag… Valores b… psi4/7… = Stat… psi4/7/10 × 1.3
## 4 alpha5/alpha8/a… Tasas de inic… Valores b… alpha5… = Stat… alpha5/8/11 × 1.3
## 5 art_cov_target Meta de cober… Valor bas… 0.95 = Stat… 0.93
## 6 prep_cov_target Cobertura obj… Muy baja … = Stat… 0.20 0.20
## 7 phi2 Tasa de entra… Valor bas… = Stat… 0.25 0.20
## 8 phi3 Tasa de aband… Valor bas… = Stat… 0.0 (a… 0.0
## 9 prep_adherence Adherencia me… Valor bas… = Stat… min(0.… min(0.80, base × …
## 10 condom_use Uso de condón… Valor bas… = Stat… = Stat… base + 0.15 (tope…
## 11 partners_year_HR Número de par… Valor bas… = Stat… = Stat… = Status quo
## 12 beh_red_dx Reducción de … Valor bas… = Stat… = Stat… = Status quo
## 13 beh_red_art Reducción de … Valor bas… = Stat… = Stat… = Status quo
## 14 beh_red_aids Reducción de … Valor bas… = Stat… = Stat… = Status quo
## # ℹ 1 more variable: combinado_fuerte <chr># Construir parámetros de intervención según escenario
make_parms_scenario <- function(base_parms, scenario) {
p <- base_parms
p$scenario <- scenario
# Escenario base: parámetros sin cambios
if (scenario == "Status quo") {
return(p)
}
# Escenario TasP 80% (test and treat)
if (scenario == "TasP 80%") {
# Aumentar proporción de tamizados a 80% y derivar tasa anual
# 80% testeados en 12 meses -> derivar tasa anual psi1
p$tested_12m_prop <- c(LR = 0.80, HR = 0.80)
p$psi1_LR <- derive_psi_from_prop(p$tested_12m_prop["LR"])
p$psi1_HR <- derive_psi_from_prop(p$tested_12m_prop["HR"])
# Aumentar tasa de diagnóstico en infectados (agudo, crónico, SIDA)
p$psi4_LR <- p$psi4_LR * 1.5
p$psi4_HR <- p$psi4_HR * 1.5
p$psi7_LR <- p$psi7_LR * 1.5
p$psi7_HR <- p$psi7_HR * 1.5
p$psi10_LR <- p$psi10_LR * 1.5
p$psi10_HR <- p$psi10_HR * 1.5
# Aumentar inicio de TAR
p$alpha5 <- p$alpha5 * 1.5
p$alpha8 <- p$alpha8 * 1.5
p$alpha11 <- p$alpha11 * 1.5
# Meta de cobertura de TAR
p$art_cov_target <- 0.95
return(p)
}
# Escenario PrEP 20%
if (scenario == "PrEP 20%") {# 20% de susceptibles en PrEP
p$prep_cov_target <- 0.20
# Mayor entrada a PrEP y casi sin abandono
p$phi2 <- 0.25
p$phi3 <- 0.0 # abandono casi nulo
# Adherencia alta pero < 100%
p$prep_adherence <- min(0.85, p$prep_adherence * 0.95)
return(p)
}
# Escenario Combinado moderado
if (scenario == "Combinado moderado") {
# Testeo moderadamente alto
p$tested_12m_prop <- c(LR = 0.70, HR = 0.75)
p$psi1_LR <- derive_psi_from_prop(p$tested_12m_prop["LR"])
p$psi1_HR <- derive_psi_from_prop(p$tested_12m_prop["HR"])
# Mayor diagnóstico
p$psi4_LR <- p$psi4_LR * 1.3
p$psi4_HR <- p$psi4_HR * 1.3
p$psi7_LR <- p$psi7_LR * 1.3
p$psi7_HR <- p$psi7_HR * 1.3
p$psi10_LR <- p$psi10_LR * 1.3
p$psi10_HR <- p$psi10_HR * 1.3
# Mayor inicio de TAR
p$alpha5 <- p$alpha5 * 1.3
p$alpha8 <- p$alpha8 * 1.3
p$alpha11 <- p$alpha11 * 1.3
p$art_cov_target <- 0.93
# PrEP moderada
p$prep_cov_target <- 0.20
p$phi2 <- 0.20
p$phi3 <- 0.0
p$prep_adherence <- min(0.80, p$prep_adherence * 0.90)
# Aumento uso de condón
p$condom_use <- min(0.80, p$condom_use + 0.15)
return(p)
}
# Escenario Combinado fuerte
if (scenario == "Combinado fuerte") {
# Testeo muy intensivo
p$tested_12m_prop <- c(LR = 0.90, HR = 0.90)
p$psi1_LR <- derive_psi_from_prop(p$tested_12m_prop["LR"])
p$psi1_HR <- derive_psi_from_prop(p$tested_12m_prop["HR"])
# Detección duplicada
p$psi4_LR <- p$psi4_LR * 2.0
p$psi4_HR <- p$psi4_HR * 2.0
p$psi7_LR <- p$psi7_LR * 2.0
p$psi7_HR <- p$psi7_HR * 2.0
p$psi10_LR <- p$psi10_LR * 2.0
p$psi10_HR <- p$psi10_HR * 2.0
# Inicio de TAR muy alto
p$alpha5 <- p$alpha5 * 1.8
p$alpha8 <- p$alpha8 * 1.8
p$alpha11 <- p$alpha11 * 1.8
p$art_cov_target <- 0.98
# PrEP alta
p$prep_cov_target <- 0.30
p$phi2 <- 0.35
p$phi3 <- 0.0
p$prep_adherence <- min(0.90, p$prep_adherence * 0.95)
# Uso de condón muy alto
p$condom_use <- min(0.95, p$condom_use + 0.25)
# Reducción del número de parejas en HR
p$partners_year["HR"] <- p$partners_year["HR"] * 0.7
# Reducciones de conducta
p$beh_red_dx <- min(0.50, p$beh_red_dx + 0.20)
p$beh_red_art <- min(0.70, p$beh_red_art + 0.20)
p$beh_red_aids <- min(0.95, p$beh_red_aids + 0.05)
return(p)
}
stop("Escenario no reconocido: ", scenario)
}
#Correr escenarios y añadir indicadores
run_scenarios_with_indicators <- function(base_parms,
init,
scenario_names = c("Status quo",
"TasP 80%",
"PrEP 20%",
"Combinado moderado",
"Combinado fuerte")) {
purrr::map_dfr(scenario_names, function(sc) {
p_sc <- make_parms_scenario(base_parms, sc)
sim <- sim_vi_hsh(parms = p_sc, init = init)
sim_i <- add_indicators_single(sim, p_sc)
sim_i$escenario <- sc
sim_i
})
}
#Resumen anual por escenario
summarise_yearly_scenarios <- function(sim_all, base_year = 2025) {
sim_all %>%
dplyr::mutate(year = floor(time) + base_year) %>%
dplyr::group_by(escenario, year) %>%
dplyr::summarise(
N_tot_py = sum(N_tot * dt, na.rm = TRUE),
new_inf_tot = sum(new_inf_year * dt, na.rm = TRUE),
inc_per1000 = ifelse(N_tot_py > 0,
new_inf_tot / N_tot_py * 1000,
NA_real_),
prev_tot_end = dplyr::last(prev_tot),
.groups = "drop"
)
}
#Infecciones evitadas vs. Status quo
calc_infections_averted <- function(yearly_all,
base_name = "Status quo") {
base <- yearly_all %>%
dplyr::filter(escenario == base_name) %>%
dplyr::select(year, new_inf_base = new_inf_tot)
yearly_all %>%
dplyr::left_join(base, by = "year") %>%
dplyr::mutate(
inf_averted = new_inf_base - new_inf_tot,
inf_averted_pct = ifelse(new_inf_base > 0,
100 * inf_averted / new_inf_base,
NA_real_)
)
}12 11.Resultados de escenarios
Se ejecutan los escenarios, se resume la dinámica anual y se calculan infecciones evitadas frente al status quo.
# Base de parámetros para escenarios
base_parms_for_scenarios <- parms_base
#Correr todos los escenarios
res_all <- run_scenarios_with_indicators(
base_parms = base_parms_for_scenarios,
init = init
)
# Resumen anual
yearly_all <- summarise_yearly_scenarios(res_all, base_year = 2025)
# Infecciones evitadas vs Status quo
yearly_all2 <- calc_infections_averted(yearly_all, base_name = "Status quo")
#Diagnóstico rápido
yearly_all2 %>%
dplyr::group_by(escenario) %>%
dplyr::summarise(
min_inf = min(inf_averted, na.rm = TRUE),
max_inf = max(inf_averted, na.rm = TRUE),
.groups = "drop"
)## # A tibble: 5 × 3
## escenario min_inf max_inf
## <chr> <dbl> <dbl>
## 1 Combinado fuerte 224. 5121.
## 2 Combinado moderado 223. 4988.
## 3 PrEP 20% 20.2 3251.
## 4 Status quo 0 0
## 5 TasP 80% 9.47 701.12.1 11.1 Prevalencia total por escenario
ggplot(yearly_all2,
aes(x = year, y = prev_tot_end * 100, color = escenario)) +
geom_line(linewidth = 1) +
labs(x = "Año calendario",
y = "Prevalencia VIH HSH (%)",
color = "Escenario") +
theme_minimal()
12.1.1 Comparación con diferentes niveles de rho: Prevalencia
library(dplyr)
library(purrr)
library(tidyr)
library(ggplot2)
# niveles de rho (como antes)
mu_val <- parms_base$mu
rho_levels <- c(
0,
mu_val,
1.5 * mu_val,
2.0 * mu_val
)
# escenarios a comparar: usar EXACTAMENTE los nombres de make_parms_scenario
escenarios <- c(
"Status quo",
"TasP 80%",
"PrEP 20%",
"Combinado moderado",
"Combinado fuerte"
)
# función auxiliar para pasar de "tiempo modelo" a "año calendario"
make_yearly_prev <- function(sim, year0 = 2023L) {
sim %>%
mutate(year = year0 + floor(time)) %>%
group_by(year) %>%
summarise(
prev_tot_end = dplyr::last(prev_tot),
.groups = "drop"
)
}
# 1. Simular todas las combinaciones escenario × rho
yearly_all_rho <- purrr::map_dfr(escenarios, function(esc) {
purrr::map_dfr(rho_levels, function(rho) {
# 👇 aquí va el cambio importante: orden correcto de argumentos
parms_sc <- make_parms_scenario(parms_base, esc)
parms_sc$rho_LR <- rho
parms_sc$rho_HR <- rho
sim <- sim_vi_hsh(parms = parms_sc, init = init)
yearly <- make_yearly_prev(sim, year0 = 2023L)
yearly %>%
mutate(
escenario = esc,
rho = rho
)
})
})
# 2. Gráfico de prevalencia con escenarios y rho
ggplot(yearly_all_rho,
aes(x = year,
y = prev_tot_end * 100,
color = escenario,
linetype = factor(round(rho, 5)))) +
geom_line(linewidth = 1) +
labs(
x = "Año calendario",
y = "Prevalencia VIH HSH (%)",
color = "Escenario",
linetype = expression(rho),
title = "Prevalencia VIH HSH según escenario e intensidad de reposición demográfica (rho)"
) +
theme_minimal() +
theme(
legend.position = "bottom"
)
12.1.2 Comparación con diferentes niveles de rho: Prevalencia SEPRARDO
library(dplyr)
library(ggplot2)
library(scales)##
## Attaching package: 'scales'## The following object is masked from 'package:purrr':
##
## discardyearly_all_rho |>
mutate(
escenario = factor(
escenario,
levels = c("Status quo", "TasP 80%", "PrEP 20%",
"Combinado moderado", "Combinado fuerte")
),
rho_lab = factor(
round(rho, 5),
levels = round(sort(unique(rho)), 5),
labels = paste0("rho = ", round(sort(unique(rho)), 5))
)
) |>
ggplot(
aes(
x = year,
y = prev_tot_end * 100,
colour = escenario,
linetype = rho_lab
)
) +
geom_line(linewidth = 1) +
facet_wrap(~ rho_lab, ncol = 2) +
scale_y_continuous(labels = label_number(accuracy = 0.1)) +
scale_colour_brewer(palette = "Set1") +
labs(
x = "Año calendario",
y = "Prevalencia VIH HSH (%)",
colour = "Escenario",
linetype = "Intensidad de reposición (rho)",
title = "Prevalencia VIH HSH según escenario e intensidad de reposición demográfica",
subtitle = "Cada panel muestra un escenario; dentro de cada uno se comparan distintos valores de rho"
) +
theme_minimal(base_size = 12) +
theme(
panel.grid.minor = element_blank(),
legend.position = "bottom",
strip.text = element_text(face = "bold"),
plot.title = element_text(face = "bold", size = 14)
)
12.2 11.2 Incidencia anual por 1.000 persona-año
ggplot(yearly_all2,
aes(x = year, y = inc_per1000, color = escenario)) +
geom_line(linewidth = 1) +
labs(x = "Año calendario",
y = "Incidencia (casos / 1.000 p-año)",
color = "Escenario") +
theme_minimal()
12.2.1 Comparación de la tasa de incidencia por diferentes rho
library(dplyr)
library(purrr)
library(tidyr)
library(ggplot2)
## Niveles de rho (mismo para LR y HR)
mu_val <- parms_base$mu
rho_levels <- c(
0,
mu_val,
1.5 * mu_val,
2.0 * mu_val
)
## Escenarios EXACTAMENTE como en make_parms_scenario()
escenarios <- c(
"Status quo",
"TasP 80%",
"PrEP 20%",
"Combinado moderado",
"Combinado fuerte"
)
## 1. Simulación completa (tiempo continuo) + indicadores, para cada escenario × rho
sim_all_rho <- purrr::map_dfr(escenarios, function(sc) {
purrr::map_dfr(rho_levels, function(rho) {
# parámetros del escenario (función original)
parms_sc <- make_parms_scenario(parms_base, sc)
# ajustar rho en LR y HR
parms_sc$rho_LR <- rho
parms_sc$rho_HR <- rho
# simular EDO
sim <- sim_vi_hsh(parms = parms_sc, init = init)
# añadir indicadores de incidencia
sim_i <- add_indicators_single(sim, parms_sc)
# añadir columnas de escenario y rho
sim_i %>%
mutate(
escenario = sc,
rho = rho
)
})
})
summarise_yearly_scenarios_rho <- function(sim_all, base_year = 2025) {
sim_all %>%
mutate(year = floor(time) + base_year) %>%
group_by(escenario, rho, year) %>%
summarise(
N_tot_py = sum(N_tot * dt, na.rm = TRUE),
new_inf_tot = sum(new_inf_year * dt, na.rm = TRUE),
inc_per1000 = ifelse(
N_tot_py > 0,
new_inf_tot / N_tot_py * 1000,
NA_real_
),
prev_tot_end = dplyr::last(prev_tot),
.groups = "drop"
)
}
# Aplicar el resumen anual
yearly_all_rho <- summarise_yearly_scenarios_rho(sim_all_rho, base_year = 2025)
yearly_all_rho_plot <- yearly_all_rho %>%
mutate(
escenario = factor(
escenario,
levels = c("Status quo", "TasP 80%", "PrEP 20%",
"Combinado moderado", "Combinado fuerte")
),
rho_lab = factor(
round(rho, 5),
levels = round(sort(unique(rho)), 5),
labels = paste0("rho = ", round(sort(unique(rho)), 5))
)
)
ggplot(yearly_all_rho_plot,
aes(x = year,
y = inc_per1000,
colour = escenario,
linetype = rho_lab)) +
geom_line(linewidth = 0.9, alpha = 0.85) +
labs(
x = "Año calendario",
y = "Incidencia VIH HSH (casos / 1.000 persona-año)",
colour = "Escenario",
linetype = "Intensidad demográfica (rho)",
title = "Incidencia VIH HSH: comparación de escenarios y niveles de reposición demográfica",
subtitle = "Todas las curvas en un solo panel"
) +
scale_colour_brewer(palette = "Set1") +
theme_minimal(base_size = 12) +
theme(
legend.position = "bottom",
panel.grid.minor = element_blank(),
plot.title = element_text(face = "bold", size = 14),
plot.subtitle = element_text(size = 11)
)
12.3 11.3 Infecciones evitadas (absolutas y en porcentaje)
ggplot(
yearly_all2 %>%
dplyr::filter(escenario != "Status quo",
!is.na(inf_averted)),
aes(x = year, y = inf_averted, fill = escenario)
) +
geom_col(position = "dodge") +
labs(x = "Año calendario",
y = "Infecciones evitadas (vs. Status quo)",
fill = "Escenario") +
theme_minimal()
calc_infections_averted_rho <- function(yearly_all, base_name = "Status quo") {
# Filtrar solo el Status quo, pero uno por cada rho
base <- yearly_all %>%
filter(escenario == base_name) %>%
select(year, rho, new_inf_base = new_inf_tot)
# Unir por año Y rho para comparar correctamente
yearly_all %>%
left_join(base, by = c("year", "rho")) %>%
mutate(
inf_averted = new_inf_base - new_inf_tot,
inf_averted_pct = ifelse(
new_inf_base > 0,
100 * inf_averted / new_inf_base,
NA_real_
)
)
}
yearly_all_rho2 <- summarise_yearly_scenarios_rho(sim_all_rho, base_year = 2025)
yearly_all_rho_averted <- calc_infections_averted_rho(yearly_all_rho2)
ggplot(
yearly_all_rho_averted %>%
filter(
escenario != "Status quo",
!is.na(inf_averted)
) %>%
mutate(
rho_lab = paste0("rho = ", round(rho, 5)),
escenario = factor(
escenario,
levels = c("TasP 80%", "PrEP 20%", "Combinado moderado", "Combinado fuerte")
)
),
aes(
x = year,
y = inf_averted,
fill = rho_lab,
colour = escenario
)
) +
geom_col(
position = position_dodge(width = 0.9),
alpha = 0.8
) +
labs(
x = "Año calendario",
y = "Infecciones evitadas (vs Status quo con mismo rho)",
fill = "Nivel de reposición (rho)",
colour = "Escenario",
title = "Infecciones evitadas por escenarios y niveles de reposición demográfica (rho)",
subtitle = "Comparación anual vs Status quo con igual rho"
) +
theme_minimal(base_size = 12) +
theme(
legend.position = "bottom",
panel.grid.minor = element_blank(),
plot.title = element_text(face = "bold"),
strip.text = element_text(face = "bold")
)
library(dplyr)
library(ggplot2)
plot_data <- yearly_all_rho_averted %>%
filter(
escenario != "Status quo",
!is.na(inf_averted)
) %>%
mutate(
rho_lab = factor(
round(rho, 5),
levels = round(sort(unique(rho)), 5),
labels = paste0("rho = ", round(sort(unique(rho)), 5))
),
escenario = factor(
escenario,
levels = c("TasP 80%", "PrEP 20%", "Combinado moderado", "Combinado fuerte")
)
)
ggplot(
plot_data,
aes(x = year, y = inf_averted, fill = rho_lab)
) +
geom_col(position = "dodge", width = 0.8) +
facet_wrap(~ escenario, ncol = 2, scales = "free_y") +
labs(
x = "Año calendario",
y = "Infecciones evitadas (vs Status quo con mismo rho)",
fill = "Nivel de reposición (rho)",
title = "Infecciones evitadas por escenario y nivel de reposición demográfica",
subtitle = "Cada panel muestra un escenario; barras agrupadas por rho"
) +
theme_minimal(base_size = 12) +
theme(
legend.position = "bottom",
panel.grid.minor = element_blank(),
strip.text = element_text(face = "bold"),
plot.title = element_text(face = "bold", size = 14)
)
plot_cum <- plot_data %>%
group_by(escenario, rho_lab) %>%
arrange(year) %>%
mutate(cum_inf_averted = cumsum(inf_averted)) %>%
ungroup()
ggplot(
plot_cum,
aes(x = year, y = cum_inf_averted,
colour = escenario, linetype = rho_lab)
) +
geom_line(linewidth = 1) +
labs(
x = "Año calendario",
y = "Infecciones evitadas acumuladas (vs Status quo)",
colour = "Escenario",
linetype = "Nivel de reposición (rho)",
title = "Infecciones evitadas acumuladas por escenario y nivel de reposición demográfica"
) +
theme_minimal(base_size = 12) +
theme(
legend.position = "bottom",
panel.grid.minor = element_blank(),
plot.title = element_text(face = "bold", size = 14)
)
12.3.0.1 Ahora en porcentaje:
ggplot(
yearly_all2 %>%
dplyr::filter(escenario != "Status quo",
!is.na(inf_averted_pct)),
aes(x = year, y = inf_averted_pct, color = escenario)
) +
geom_line(linewidth = 1) +
labs(x = "Año calendario",
y = "% infecciones evitadas",
color = "Escenario") +
theme_minimal()
library(dplyr)
library(ggplot2)
plot_data <- yearly_all_rho_averted %>%
filter(
escenario != "Status quo",
!is.na(inf_averted_pct)
) %>%
mutate(
escenario = factor(
escenario,
levels = c("TasP 80%", "PrEP 20%", "Combinado moderado", "Combinado fuerte")
),
rho_lab = factor(
round(rho, 5),
levels = round(sort(unique(rho)), 5),
labels = paste0("rho = ", round(sort(unique(rho)), 5))
)
)
ggplot(
plot_data,
aes(
x = year,
y = inf_averted_pct,
color = escenario,
linetype = rho_lab
)
) +
geom_line(linewidth = 1) +
labs(
x = "Año calendario",
y = "% infecciones evitadas respecto al Status quo con mismo rho",
color = "Escenario",
linetype = "Nivel de reposición (rho)",
title = "% infecciones evitadas según escenario y niveles de reposición demográfica"
) +
theme_minimal(base_size = 12) +
theme(
legend.position = "bottom",
panel.grid.minor = element_blank(),
plot.title = element_text(face = "bold", size = 14)
)
12.4 11.4 Tabla resumen de escenarios
#Total de infecciones en el escenario base (Status quo)
base_total_inf <- yearly_all2 %>%
dplyr::filter(escenario == "Status quo") %>%
dplyr::summarise(total_inf = sum(new_inf_tot, na.rm = TRUE)) %>%
dplyr::pull(total_inf)
# Último año disponible por escenario
last_year <- yearly_all2 %>%
dplyr::group_by(escenario) %>%
dplyr::slice_max(year, n = 1, with_ties = FALSE) %>%
dplyr::ungroup() %>%
dplyr::select(
escenario,
anio_fin = year,
prev_fin = prev_tot_end,
inc_fin_1000 = inc_per1000
)
#Infecciones totales en el horizonte por escenario
inf_tot <- yearly_all2 %>%
dplyr::group_by(escenario) %>%
dplyr::summarise(
inf_totales = sum(new_inf_tot, na.rm = TRUE),
.groups = "drop"
)
# Unir todo y calcular infecciones evitadas
tabla_escenarios <- last_year %>%
dplyr::left_join(inf_tot, by = "escenario") %>%
dplyr::mutate(
prev_fin_pct = prev_fin * 100,
inf_totales_1e4 = inf_totales / 1e4,
inf_ev_total = dplyr::if_else(
escenario == "Status quo",
0,
base_total_inf - inf_totales
),
inf_ev_1e4 = inf_ev_total / 1e4,
inf_ev_pct = dplyr::if_else(
escenario == "Status quo",
NA_real_,
100 * inf_ev_total / base_total_inf
)
) %>%
dplyr::arrange(match(
escenario,
c("Status quo",
"TasP 80%",
"PrEP 20%",
"Combinado moderado",
"Combinado fuerte")
)) %>%
dplyr::transmute(
Escenario = escenario,
`Año` = anio_fin,
`Prevalencia (%)` = round(prev_fin_pct, 2),
`Incidencia (/1000 p-año)` = round(inc_fin_1000, 3),
`Inf. totales (×10.000)` = round(inf_totales_1e4, 2),
`Inf. evitadas (×10.000)` = round(inf_ev_1e4, 2),
`Inf. evitadas (%)` = round(inf_ev_pct, 1)
)
tabla_escenarios## # A tibble: 5 × 7
## Escenario Año `Prevalencia (%)` `Incidencia (/1000 p-año)`
## <chr> <dbl> <dbl> <dbl>
## 1 Status quo 2065 10.4 9.51
## 2 TasP 80% 2065 9.66 8.88
## 3 PrEP 20% 2065 3.28 2.34
## 4 Combinado moderado 2065 1.01 0.77
## 5 Combinado fuerte 2065 0.22 0.704
## # ℹ 3 more variables: `Inf. totales (×10.000)` <dbl>,
## # `Inf. evitadas (×10.000)` <dbl>, `Inf. evitadas (%)` <dbl>library(dplyr)
# 1) Total de infecciones en el escenario base (Status quo) POR rho
base_total_inf_rho <- yearly_all_rho2 %>%
filter(escenario == "Status quo") %>%
group_by(rho) %>%
summarise(
base_total_inf = sum(new_inf_tot, na.rm = TRUE),
.groups = "drop"
)
# 2) Último año disponible por escenario Y rho
last_year_rho <- yearly_all_rho2 %>%
group_by(escenario, rho) %>%
slice_max(year, n = 1, with_ties = FALSE) %>%
ungroup() %>%
select(
escenario,
rho,
anio_fin = year,
prev_fin = prev_tot_end,
inc_fin_1000 = inc_per1000
)
# 3) Infecciones totales en el horizonte por escenario Y rho
inf_tot_rho <- yearly_all_rho2 %>%
group_by(escenario, rho) %>%
summarise(
inf_totales = sum(new_inf_tot, na.rm = TRUE),
.groups = "drop"
)
# 4) Unir todo y calcular infecciones evitadas POR escenario Y rho
tabla_escenarios_rho <- last_year_rho %>%
left_join(inf_tot_rho, by = c("escenario", "rho")) %>%
left_join(base_total_inf_rho, by = "rho") %>%
mutate(
prev_fin_pct = prev_fin * 100,
inf_totales_1e4 = inf_totales / 1e4,
inf_ev_total = if_else(
escenario == "Status quo",
0,
base_total_inf - inf_totales
),
inf_ev_1e4 = inf_ev_total / 1e4,
inf_ev_pct = if_else(
escenario == "Status quo",
NA_real_,
100 * inf_ev_total / base_total_inf
)
) %>%
arrange(
rho,
match(
escenario,
c("Status quo",
"TasP 80%",
"PrEP 20%",
"Combinado moderado",
"Combinado fuerte")
)
) %>%
transmute(
rho = round(rho, 5),
Escenario = escenario,
`Año` = anio_fin,
`Prevalencia (%)` = round(prev_fin_pct, 2),
`Incidencia (/1000 p-año)`= round(inc_fin_1000, 3),
`Inf. totales (×10.000)` = round(inf_totales_1e4, 2),
`Inf. evitadas (×10.000)` = round(inf_ev_1e4, 2),
`Inf. evitadas (%)` = round(inf_ev_pct, 1)
)
tabla_escenarios_rho## # A tibble: 20 × 8
## rho Escenario Año `Prevalencia (%)` `Incidencia (/1000 p-año)`
## <dbl> <chr> <dbl> <dbl> <dbl>
## 1 0 Status quo 2065 10.4 9.51
## 2 0 TasP 80% 2065 9.66 8.88
## 3 0 PrEP 20% 2065 3.28 2.34
## 4 0 Combinado moderado 2065 1.01 0.77
## 5 0 Combinado fuerte 2065 0.22 0.704
## 6 0.0052 Status quo 2065 10.4 9.93
## 7 0.0052 TasP 80% 2065 9.62 9.27
## 8 0.0052 PrEP 20% 2065 3.8 3.18
## 9 0.0052 Combinado moderado 2065 1.39 1.24
## 10 0.0052 Combinado fuerte 2065 0.48 0.51
## 11 0.0078 Status quo 2065 10.3 10.1
## 12 0.0078 TasP 80% 2065 9.6 9.46
## 13 0.0078 PrEP 20% 2065 4.04 3.58
## 14 0.0078 Combinado moderado 2065 1.56 1.47
## 15 0.0078 Combinado fuerte 2065 0.56 0.488
## 16 0.0104 Status quo 2065 10.3 10.3
## 17 0.0104 TasP 80% 2065 9.58 9.64
## 18 0.0104 PrEP 20% 2065 4.27 3.97
## 19 0.0104 Combinado moderado 2065 1.72 1.71
## 20 0.0104 Combinado fuerte 2065 0.62 0.508
## # ℹ 3 more variables: `Inf. totales (×10.000)` <dbl>,
## # `Inf. evitadas (×10.000)` <dbl>, `Inf. evitadas (%)` <dbl>tabla_escenarios_rho_panel <- tabla_escenarios_rho %>%
mutate(
`Prev (%)` = round(`Prevalencia (%)`, 1),
`Inc (/1000)` = round(`Incidencia (/1000 p-año)`, 2),
`Inf tot (10k)` = round(`Inf. totales (×10.000)`, 1),
`Inf ev (10k)` = round(`Inf. evitadas (×10.000)`, 1),
`% ev` = round(`Inf. evitadas (%)`, 1)
) %>%
select(rho, Escenario, Año,
`Prev (%)`, `Inc (/1000)`,
`Inf tot (10k)`, `Inf ev (10k)`, `% ev`)
tabla_escenarios_rho_panel## # A tibble: 20 × 8
## rho Escenario Año `Prev (%)` `Inc (/1000)` `Inf tot (10k)`
## <dbl> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 0 Status quo 2065 10.4 9.51 12
## 2 0 TasP 80% 2065 9.7 8.88 11
## 3 0 PrEP 20% 2065 3.3 2.34 5.1
## 4 0 Combinado moderado 2065 1 0.77 1.2
## 5 0 Combinado fuerte 2065 0.2 0.7 0.2
## 6 0.0052 Status quo 2065 10.4 9.93 14.1
## 7 0.0052 TasP 80% 2065 9.6 9.27 12.9
## 8 0.0052 PrEP 20% 2065 3.8 3.18 6.6
## 9 0.0052 Combinado moderado 2065 1.4 1.24 1.9
## 10 0.0052 Combinado fuerte 2065 0.5 0.51 0.4
## 11 0.0078 Status quo 2065 10.3 10.1 15.3
## 12 0.0078 TasP 80% 2065 9.6 9.46 14.1
## 13 0.0078 PrEP 20% 2065 4 3.58 7.4
## 14 0.0078 Combinado moderado 2065 1.6 1.47 2.3
## 15 0.0078 Combinado fuerte 2065 0.6 0.49 0.5
## 16 0.0104 Status quo 2065 10.3 10.3 16.6
## 17 0.0104 TasP 80% 2065 9.6 9.64 15.3
## 18 0.0104 PrEP 20% 2065 4.3 3.97 8.4
## 19 0.0104 Combinado moderado 2065 1.7 1.71 2.7
## 20 0.0104 Combinado fuerte 2065 0.6 0.51 0.6
## # ℹ 2 more variables: `Inf ev (10k)` <dbl>, `% ev` <dbl>13 12. Estimación de un proxy de (R_{})
En modelos compartimentales determinísticos, el número reproductivo efectivo \(R_{\text{efectivo}}(t)\) representa el número esperado de infecciones secundarias producidas por un individuo infeccioso típico en el estado actual de la población (fracción de susceptibles, coberturas de TAR, PrEP, uso de condón, etc.).
En este modelo, las nuevas infecciones instantáneas se calculan como:
\[ \text{new\_inf}(t) = \lambda_{1,LR}(t)\,S_{1,LR}(t) + \lambda_{2,LR}(t)\,S_{2,LR}(t) + \lambda_{3,LR}(t)\,S_{3,LR}(t) + \lambda_{1,HR}(t)\,S_{1,HR}(t) + \lambda_{2,HR}(t)\,S_{2,HR}(t) + \lambda_{3,HR}(t)\,S_{3,HR}(t), \]
y el tamaño total de la población es
\[ N_{\text{tot}}(t) = N_{LR}(t) + N_{HR}(t). \]
La incidencia instantánea (casos por persona-año) viene dada por
\[ \text{inc}(t) = \frac{\text{new\_inf}(t)}{N_{\text{tot}}(t)}. \]
Dado que las intervenciones (TasP, PrEP, etc.) modifican sobre todo la dinámica a lo largo del tiempo y no únicamente el estado inicial (\(t = 0\)), se define un proxy de \(R_{\text{efectivo},0}\) basado en la incidencia promedio durante el primer año de simulación.
Sea \(\theta\) un conjunto de parámetros (por ejemplo, un escenario de intervención) y \(\theta_{\text{Status quo}}\) los parámetros del escenario de referencia. Definimos la incidencia media en el intervalo \([0,1]\) como:
\[ I_1(\theta) = \frac{ \displaystyle \int_0^1 \text{new\_inf}(t;\theta)\,dt }{ \displaystyle \int_0^1 N_{\text{tot}}(t;\theta)\,dt }. \]
El proxy relativo del número reproductivo efectivo en el año inicial se define entonces como
\[ R_{\text{rel},0}(\theta) = \frac{I_1(\theta)} {I_1(\theta_{\text{Status quo}})}. \]
Bajo el supuesto de proporcionalidad entre incidencia y número reproductivo efectivo, este cociente se interpreta como
\[ R_{\text{rel},0}(\theta) \approx \frac{ R_{\text{efectivo},0}(\theta) }{ R_{\text{efectivo},0}(\text{Status quo}) }. \]
Así:
- \(R_{\text{rel},0}(\theta) = 1\): mismo potencial de transmisión que el Status quo durante el primer año.
- \(R_{\text{rel},0}(\theta) = 0{,}7\): reducción del 30 % en la incidencia media del primer año.
- \(R_{\text{rel},0}(\theta) = 1{,}4\): incremento del 40 %.
De forma complementaria, se define el tiempo hasta eliminación como el primer instante en que la incidencia instantánea cae por debajo de un umbral \(\tau\) (casos por persona-año). En lugar de fijar un valor absoluto, se usa un umbral relativo respecto a la incidencia media de referencia:
\[ \tau = \alpha\,I_1(\theta_{\text{Status quo}}), \qquad 0 < \alpha < 1. \]
En este trabajo se toma \(\alpha = 0{,}1\), es decir, el umbral se fija en el 10 % de la incidencia media del primer año en el escenario de referencia. Entonces:
\[ T_{\text{elim}}(\theta) = \min\left\{ t \ge 0 : \text{inc}\bigl(t;\theta\bigr) < \tau \right\}. \]
Si la incidencia no alcanza dicho umbral en el horizonte simulado, se considera \(T_{\text{elim}}(\theta) = \text{NA}\).
| Aspecto | Enfoque: ABSOLUTO (Basado en dinámica de infectados) | Enfoque: RELATIVO (Basado en la incidencia del primer año) |
|---|---|---|
| Fórmula base | $$$$Refectivo(t)=G(t)F(t)$$$$ | \[ Rrel,0(θ)=I1(Status quo)I1(θ) \] |
| Interpreta | Derivar \[dI/dt\], definir flujos completos | Solo incidencia del primer año |
| Qué mide | Mide transmisión real del sistema dinámico | Mide reducción relativa de transmisión bajo un escenario |
| Más sensible a | Cambios en progresión natural, mortalidad, entrada | Cambios en medidas de prevención, testeo, TAR |
| Dificultad | Más complejo (requiere definir \[G(t)\]) | Muy simple y estable |
| Interpretación | Si \[R_{ef} < 1\] → epidemia decrece | Si $$ R_{rel,0} < 1$$ → menor incidencia que status quo |
13.1 12.1 Incidencia de referencias Status quo
#Simulación con parámetros base (Status quo)
sim_ref <- sim_vi_hsh(parms = parms_base, init = init)
sim_ref_i <- add_indicators_single(sim_ref, parms_base)
#Ventana de tiempo [0,1) años
year1_ref <- sim_ref_i$time >= 0 & sim_ref_i$time < 1
#Casos incidentes acumulados en el primer año
cases_ref_year1 <- sum(
sim_ref_i$new_inf_year[year1_ref] * sim_ref_i$dt[year1_ref],
na.rm = TRUE
)
#Persona-años acumulados en el primer año
pya_ref_year1 <- sum(
sim_ref_i$N_tot[year1_ref] * sim_ref_i$dt[year1_ref],
na.rm = TRUE
)
#Incidencia media en el primer año (Status quo)
inc1_ref <- cases_ref_year1 / pya_ref_year1
inc1_ref # para inspeccionar el valor## [1] 0.002513.2 12.2 Función para calcular R_rel,0 y tiempo hasta eliminación:
La función compute_R_rel_elim() recibe:
parms: lista de parámetros (para un escenario concreto)
init : condiciones iniciales del sistema
inc0_ref: incidencia inicial de referencia (Status quo)
inc_threshold: umbral de eliminación (casos / persona-año)
Devuelve una lista con:
R_rel_0: proxy relativo de R_efectivo,0
R_rel_0 = inc(0; parms) / inc_ref(0)
inc0: incidencia instantánea inicial inc(0; parms)
years_to_elim: tiempo hasta eliminación (años)
compute_R_rel_elim <- function(parms,
init,
inc1_ref,
alpha = 0.5) {
out <- tryCatch(
{
# 1) Simular el modelo con estos parámetros
sim <- sim_vi_hsh(parms = parms, init = init)
sim_i <- add_indicators_single(sim, parms)
# 2) Incidencia media en el primer año para este escenario
year1 <- sim_i$time >= 0 & sim_i$time < 1
cases_year1 <- sum(
sim_i$new_inf_year[year1] * sim_i$dt[year1],
na.rm = TRUE
)
pya_year1 <- sum(
sim_i$N_tot[year1] * sim_i$dt[year1],
na.rm = TRUE
)
inc1 <- cases_year1 / pya_year1
# 3) Proxy relativo de R_efectivo,0
R_rel_0 <- if (is.finite(inc1_ref) && inc1_ref > 0) {
inc1 / inc1_ref
} else {
NA_real_
}
# 4) Umbral de eliminación (fracción alpha de inc1_ref)
inc_threshold <- alpha * inc1_ref
# 5) Tiempo hasta eliminación:
# T_elim = min{ t : inc_rate_inst(t) < inc_threshold }
idx_elim <- which(sim_i$inc_rate_inst < inc_threshold)
years_to_elim <- if (length(idx_elim) == 0) {
NA_real_ # no se alcanza eliminación en el horizonte simulado
} else {
sim_i$time[min(idx_elim)]
}
list(
R_rel_0 = R_rel_0,
inc1 = inc1,
years_to_elim = years_to_elim
)
},
error = function(e) {
warning(
"compute_R_rel_elim: fallo de deSolve para este set de parámetros: ",
conditionMessage(e)
)
list(
R_rel_0 = NA_real_,
inc1 = NA_real_,
years_to_elim = NA_real_
)
}
)
out
}13.3 12.3 Comparaciones de escenarios
## 11.3 Comparaciones de escenarios para distintos niveles de α
# Escenarios que se van a comparar
escenarios <- c(
"Status quo",
"TasP 80%",
"PrEP 20%",
"Combinado moderado",
"Combinado fuerte"
)
# Valores de alpha (fracción del I1_ref que define el umbral de eliminación)
alpha_vec <- c(0.10, 0.15, 0.20, 0.25, 0.30, 0.40, 0.50)
# Grid escenario x alpha y cálculo de R_rel_0 e indicador de eliminación
tabla_Reff_alpha <- tidyr::expand_grid(
escenario = escenarios,
alpha = alpha_vec
) %>%
dplyr::rowwise() %>%
dplyr::mutate(
# 1) Parámetros del escenario concreto
parms_sc = list(make_parms_scenario(parms_base, escenario)),
# 2) Cálculo de R_rel_0, inc1 y years_to_elim para este alpha
stats = list(
compute_R_rel_elim(
parms = parms_sc,
init = init,
inc1_ref = inc1_ref,
alpha = alpha
)
),
R_rel_0 = stats$R_rel_0,
inc1 = stats$inc1,
years_to_elim = stats$years_to_elim
) %>%
dplyr::ungroup() %>%
dplyr::mutate(
alpha_pct = alpha * 100 # para mostrar el nivel de alpha en %
) %>%
dplyr::select(
Escenario = escenario,
`α (%)` = alpha_pct,
R_rel_0,
inc1,
years_to_elim
)
tabla_Reff_alpha## # A tibble: 35 × 5
## Escenario `α (%)` R_rel_0 inc1 years_to_elim
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Status quo 10 1 0.00250 NA
## 2 Status quo 15 1 0.00250 NA
## 3 Status quo 20 1 0.00250 NA
## 4 Status quo 25 1 0.00250 NA
## 5 Status quo 30 1 0.00250 NA
## 6 Status quo 40 1 0.00250 NA
## 7 Status quo 50 1 0.00250 NA
## 8 TasP 80% 10 0.946 0.00237 NA
## 9 TasP 80% 15 0.946 0.00237 NA
## 10 TasP 80% 20 0.946 0.00237 NA
## # ℹ 25 more rows## 11.4 Tabla de comparación por escenario (filas = indicadores, columnas = escenarios)
# Asegurarnos de tener dplyr / tidyr cargados
# library(dplyr)
# library(tidyr)
# 1) Tabla wide para la incidencia media del primer año (no depende de α)
tabla_inc1_wide <- tabla_Reff_alpha %>%
dplyr::group_by(Escenario) %>%
dplyr::summarise(
inc1 = dplyr::first(inc1), # es el mismo valor para todos los α
.groups = "drop"
) %>%
dplyr::mutate(
indicador = "Incidencia media año 1 (I1)"
) %>%
dplyr::select(indicador, Escenario, valor = inc1) %>%
tidyr::pivot_wider(
names_from = Escenario,
values_from = valor
)
# 2) Tabla wide para R_rel_0, con una fila por cada α
tabla_Rrel0_wide <- tabla_Reff_alpha %>%
dplyr::mutate(
indicador = paste0("R_rel_0 (α = ", `α (%)`, "%)")
) %>%
dplyr::select(indicador, Escenario, valor = R_rel_0) %>%
dplyr::distinct() %>%
tidyr::pivot_wider(
names_from = Escenario,
values_from = valor
)
# 3) Tabla wide para T_elim (years_to_elim), una fila por cada α
tabla_Telim_wide <- tabla_Reff_alpha %>%
dplyr::mutate(
indicador = paste0("T_elim (α = ", `α (%)`, "%)")
) %>%
dplyr::select(indicador, Escenario, valor = years_to_elim) %>%
dplyr::distinct() %>%
tidyr::pivot_wider(
names_from = Escenario,
values_from = valor
)
# 4) Unir todo en una sola tabla:
# - Primero la incidencia media
# - Luego todas las filas de R_rel_0 por α
# - Luego todas las filas de T_elim por α
tabla_indicadores_escenarios <- dplyr::bind_rows(
tabla_inc1_wide,
tabla_Rrel0_wide,
tabla_Telim_wide
)
tabla_indicadores_escenarios## # A tibble: 15 × 6
## indicador `Combinado fuerte` `Combinado moderado` `PrEP 20%` `Status quo`
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Incidencia m… 0.00117 0.00192 0.00246 0.00250
## 2 R_rel_0 (α =… 0.466 0.767 0.984 1
## 3 R_rel_0 (α =… 0.466 0.767 0.984 1
## 4 R_rel_0 (α =… 0.466 0.767 0.984 1
## 5 R_rel_0 (α =… 0.466 0.767 0.984 1
## 6 R_rel_0 (α =… 0.466 0.767 0.984 1
## 7 R_rel_0 (α =… 0.466 0.767 0.984 1
## 8 R_rel_0 (α =… 0.466 0.767 0.984 1
## 9 T_elim (α = … 4.1 NA NA NA
## 10 T_elim (α = … 3 13.7 NA NA
## 11 T_elim (α = … 2.3 9.5 NA NA
## 12 T_elim (α = … 1.7 7.2 NA NA
## 13 T_elim (α = … 1.3 5.6 NA NA
## 14 T_elim (α = … 0.7 3.5 NA NA
## 15 T_elim (α = … 0.3 2.1 NA NA
## # ℹ 1 more variable: `TasP 80%` <dbl>13.3.1 Ordenar escenarios y preparar etiquetas
## Ordenar escenarios y crear etiqueta de alpha
tabla_Reff_alpha_plot <- tabla_Reff_alpha %>%
dplyr::mutate(
Escenario = factor(
Escenario,
levels = c("Status quo",
"TasP 80%",
"PrEP 20%",
"Combinado moderado",
"Combinado fuerte")
),
alpha_lab = paste0(`α (%)`, "%")
)Heatmap elegante de 𝑅rel,0Rrel,0 (proxy relativo de R efectivo)
## Heatmap de R_rel_0 por escenario y nivel de α
ggplot(tabla_Reff_alpha_plot,
aes(x = alpha_lab, y = Escenario, fill = R_rel_0)) +
geom_tile(color = "white", linewidth = 0.4) +
geom_text(aes(label = sprintf("%.2f", R_rel_0)),
color = "white", size = 3) +
scale_fill_viridis_c(
option = "plasma",
name = expression(R[rel,0]),
limits = c(min(tabla_Reff_alpha_plot$R_rel_0, na.rm = TRUE),
max(tabla_Reff_alpha_plot$R_rel_0, na.rm = TRUE))
) +
labs(
x = expression(alpha~~"(umbral como % de"~~I[1]^ref*")"),
y = "Escenario",
title = expression("Proxy relativo de"~~R[efectivo,0]),
subtitle = "Valores < 1 indican menor potencial de transmisión que el Status quo"
) +
theme_minimal(base_size = 12) +
theme(
panel.grid = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
plot.title = element_text(face = "bold"),
legend.position = "right",
legend.title = element_text(face = "bold")
)
Heatmap de tiempo hasta eliminación 𝑇elimTelim
## Heatmap de años hasta eliminación por escenario y α
ggplot(tabla_Reff_alpha_plot,
aes(x = alpha_lab, y = Escenario, fill = years_to_elim)) +
geom_tile(color = "white", linewidth = 0.4) +
geom_text(
aes(label = ifelse(is.na(years_to_elim),
"NA",
sprintf("%.1f", years_to_elim))),
color = "white", size = 3
) +
scale_fill_viridis_c(
option = "magma",
name = "Años hasta eliminación",
na.value = "grey70"
) +
labs(
x = expression(alpha~~"(umbral como % de"~~I[1]^ref*")"),
y = "Escenario",
title = expression("Tiempo hasta eliminación"~~T[elim]),
subtitle = "NA indica que no se alcanza el umbral en el horizonte simulado"
) +
theme_minimal(base_size = 12) +
theme(
panel.grid = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
plot.title = element_text(face = "bold"),
legend.position = "right",
legend.title = element_text(face = "bold")
)
14 13. Estimación de un \(R_{\text{efectivo}}\) absoluto y tiempo hasta eliminación
En este modelo, agregamos todos los compartimentos infectados (A, C, D en LR y HR) en una sola cantidad:
\[ I(t) = I_{LR}(t) + I_{HR}(t), \]
donde \(I_{LR}(t)\) y \(I_{HR}(t)\) son las sumas de los compartimentos infectados en cada grupo de riesgo.
La dinámica agregada de los infectados puede escribirse de forma genérica como:
\[ \frac{dI(t)}{dt} = F(t) - G(t), \]
donde:
- \(F(t)\) = flujo de nuevas infecciones que entran al conjunto de compartimentos infectados en el tiempo \(t\).
- \(G(t)\) = flujo de salida de infectados en el tiempo \(t\), debido a progresión, muerte u otras salidas que los sacan del conjunto de estados infectantes.
Por construcción del modelo:
- \(F(t)\) es exactamente \(\text{new\_inf}(t)\), el número de nuevas infecciones instantáneas (tasas por año).
- Entonces, a partir de la identidad anterior:
\[ G(t) = F(t) - \frac{dI(t)}{dt} = \text{new\_inf}(t) - \frac{dI(t)}{dt}. \]
Definimos el número reproductivo efectivo instantáneo como el cociente:
\[ R_{\text{efectivo}}(t) = \frac{F(t)}{G(t)} = \frac{\text{new\_inf}(t)} {\text{new\_inf}(t) - \dfrac{dI(t)}{dt}}. \]
Intuitivamente:
- Si la epidemia está en equilibrio (\(dI/dt \approx 0\)), entonces \(\text{new\_inf}(t) \approx G(t)\) y \(R_{\text{efectivo}}(t) \approx 1\).
- Si la epidemia está creciendo (\(dI/dt > 0\)), entonces \(\text{new\_inf}(t) > G(t)\) y \(R_{\text{efectivo}}(t) > 1\).
- Si la epidemia está decreciendo (\(dI/dt < 0\)), entonces \(\text{new\_inf}(t) < G(t)\) y \(R_{\text{efectivo}}(t) < 1\).
Para obtener un resumen del comportamiento inicial de cada escenario, se define un \(R_{\text{efectivo},0}\) promedio como el promedio ponderado de \(R_{\text{efectivo}}(t)\) durante el primer año:
\[ R_{\text{efectivo},0}(\theta) = \frac{ \displaystyle \int_0^1 R_{\text{efectivo}}(t;\theta)\,dt }{ \displaystyle \int_0^1 dt } \approx \frac{ \sum_{t_i \in [0,1)} R_{\text{efectivo}}(t_i;\theta)\,\Delta t_i }{ \sum_{t_i \in [0,1)} \Delta t_i }. \]
Como antes, la incidencia instantánea se define como
\[ \text{inc}(t) = \frac{\text{new\_inf}(t)}{N_{\text{tot}}(t)}, \qquad N_{\text{tot}}(t) = N_{LR}(t) + N_{HR}(t). \]
y la incidencia media en el primer año para un escenario con parámetros \(\theta\) se define como:
\[ I_1(\theta) = \frac{ \displaystyle \int_0^1 \text{new\_inf}(t;\theta)\,dt }{ \displaystyle \int_0^1 N_{\text{tot}}(t;\theta)\,dt }. \]
Sea \(\theta_{\text{Status quo}}\) el escenario de referencia; su incidencia media en el primer año es \(I_1^{\text{ref}} = I_1(\theta_{\text{Status quo}})\).
Con estos elementos, el tiempo hasta eliminación se define mediante un umbral relativo de incidencia:
\[ \tau = \alpha\,I_1^{\text{ref}}, \qquad 0 < \alpha < 1, \]
y
\[ T_{\text{elim}}(\theta) = \min\left\{ t \ge 0 : \text{inc}(t;\theta) < \tau \right\}. \]
Si la incidencia no desciende por debajo de \(\tau\) en el horizonte de simulación, se toma \(T_{\text{elim}}(\theta) = \text{NA}\).
14.1 13.1 Incidencia media de referencia en Status quo
#Incidencia media del primer año en Status quo
sim_ref <- sim_vi_hsh(parms = parms_base, init = init)
sim_ref_i <- add_indicators_single(sim_ref, parms_base)
year1_ref <- sim_ref_i$time >= 0 & sim_ref_i$time < 1
cases_ref_year1 <- sum(
sim_ref_i$new_inf_year[year1_ref] * sim_ref_i$dt[year1_ref],
na.rm = TRUE
)
pya_ref_year1 <- sum(
sim_ref_i$N_tot[year1_ref] * sim_ref_i$dt[year1_ref],
na.rm = TRUE
)
inc1_ref <- cases_ref_year1 / pya_ref_year1
inc1_ref # incidencia media año 1 en Status quo## [1] 0.002514.2 13.2 Función para 𝑅efectivo,0 Refectivo,0 absoluto y 𝑇elimTeli:
Esta función:
Simula el modelo con un conjunto de parámetros ‘parms’.
Calcula R_efectivo(t) a partir de: R_efectivo(t) = new_inf(t) / (new_inf(t) - dI/dt) donde I(t) es el total de infectados.
Calcula R_efectivo,0 como promedio en [0, 1) años.
Calcula la incidencia media en el primer año (inc1).
Calcula el tiempo hasta eliminación según: inc(t) < alpha * inc1_ref
compute_R_eff_abs_elim <- function(parms,
init,
inc1_ref,
alpha = 0.1) {
out <- tryCatch(
{
# 1) Simular el modelo con estos parámetros
sim <- sim_vi_hsh(parms = parms, init = init)
sim_i <- add_indicators_single(sim, parms)
# 2) Extraer cantidades necesarias
time <- sim_i$time
dt_vec <- sim_i$dt
I_tot <- sim_i$I_tot # total de infectados
new_inf <- sim_i$new_inf_year # flujo instantáneo de nuevas infecciones (tasa/año)
# 3) Aproximar dI/dt numéricamente
dI_dt <- c(
diff(I_tot) / diff(time),
tail(diff(I_tot) / diff(time), 1)
)
# 4) Flujo de salida de infectados: G(t) = new_inf - dI/dt
removal <- new_inf - dI_dt
# 5) Número reproductivo efectivo instantáneo
R_eff_t <- ifelse(removal > 0,
new_inf / removal,
NA_real_)
# 6) Promedio de R_efectivo en el primer año [0,1)
year1 <- time >= 0 & time < 1
if (any(year1 & is.finite(R_eff_t))) {
R_eff_0 <- sum(
R_eff_t[year1] * dt_vec[year1],
na.rm = TRUE
) / sum(
dt_vec[year1][is.finite(R_eff_t[year1])],
na.rm = TRUE
)
} else {
R_eff_0 <- NA_real_
}
# 7) Incidencia media en el primer año para este escenario
cases_year1 <- sum(
new_inf[year1] * dt_vec[year1],
na.rm = TRUE
)
pya_year1 <- sum(
sim_i$N_tot[year1] * dt_vec[year1],
na.rm = TRUE
)
inc1 <- cases_year1 / pya_year1
# 8) Umbral de eliminación
inc_threshold <- alpha * inc1_ref
# 9) Tiempo hasta eliminación
idx_elim <- which(sim_i$inc_rate_inst < inc_threshold)
years_to_elim <- if (length(idx_elim) == 0) {
NA_real_ # no alcanza eliminación en el horizonte simulado
} else {
sim_i$time[min(idx_elim)]
}
list(
R_eff_0 = R_eff_0,
inc1 = inc1,
years_to_elim = years_to_elim
)
},
error = function(e) {
warning(
"compute_R_eff_abs_elim: fallo de deSolve para este set de parámetros: ",
conditionMessage(e)
)
list(
R_eff_0 = NA_real_,
inc1 = NA_real_,
years_to_elim = NA_real_
)
}
)
out
}14.3 13.3 Comparación de escenarios y niveles de α:
## Escenarios y alphas
escenarios <- c(
"Status quo",
"TasP 80%",
"PrEP 20%",
"Combinado moderado",
"Combinado fuerte"
)
alpha_vec <- c(0.10, 0.15, 0.20, 0.25, 0.30, 0.40, 0.50)
## Tabla escenario x alpha con R_eff_0, inc1 y T_elim
tabla_Reff_abs_alpha <- tidyr::expand_grid(
escenario = escenarios,
alpha = alpha_vec
) %>%
dplyr::rowwise() %>%
dplyr::mutate(
# parámetros del escenario
parms_sc = list(make_parms_scenario(parms_base, escenario)),
# cálculo de R_efectivo,0, inc1 y T_elim para este alpha
stats = list(
compute_R_eff_abs_elim(
parms = parms_sc,
init = init,
inc1_ref = inc1_ref,
alpha = alpha
)
),
R_eff_0 = stats$R_eff_0,
inc1 = stats$inc1,
years_to_elim = stats$years_to_elim
) %>%
dplyr::ungroup() %>%
dplyr::mutate(
alpha_pct = alpha * 100 # nivel de alpha en porcentaje
) %>%
dplyr::select(
Escenario = escenario,
alpha_pct,
R_eff_0,
inc1,
years_to_elim
)
tabla_Reff_abs_alpha## # A tibble: 35 × 5
## Escenario alpha_pct R_eff_0 inc1 years_to_elim
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Status quo 10 0.0822 0.00250 NA
## 2 Status quo 15 0.0822 0.00250 NA
## 3 Status quo 20 0.0822 0.00250 NA
## 4 Status quo 25 0.0822 0.00250 NA
## 5 Status quo 30 0.0822 0.00250 NA
## 6 Status quo 40 0.0822 0.00250 NA
## 7 Status quo 50 0.0822 0.00250 NA
## 8 TasP 80% 10 0.0791 0.00237 NA
## 9 TasP 80% 15 0.0791 0.00237 NA
## 10 TasP 80% 20 0.0791 0.00237 NA
## # ℹ 25 more rows## 1) Incidencia media del primer año (no depende de alpha)
tabla_inc1_wide <- tabla_Reff_abs_alpha %>%
dplyr::group_by(Escenario) %>%
dplyr::summarise(
inc1 = dplyr::first(inc1), # es el mismo valor para todos los alpha
.groups = "drop"
) %>%
dplyr::mutate(
indicador = "Incidencia media año 1 (I1)"
) %>%
dplyr::select(indicador, Escenario, valor = inc1) %>%
tidyr::pivot_wider(
names_from = Escenario,
values_from = valor
)
## 2) R_efectivo,0 por alpha
tabla_Reff0_wide <- tabla_Reff_abs_alpha %>%
dplyr::mutate(
indicador = paste0("R_efectivo,0 (α = ", alpha_pct, "%)")
) %>%
dplyr::select(indicador, Escenario, valor = R_eff_0) %>%
dplyr::distinct() %>%
tidyr::pivot_wider(
names_from = Escenario,
values_from = valor
)
## 3) T_elim por alpha
tabla_Telim_wide <- tabla_Reff_abs_alpha %>%
dplyr::mutate(
indicador = paste0("T_elim (α = ", alpha_pct, "%)")
) %>%
dplyr::select(indicador, Escenario, valor = years_to_elim) %>%
dplyr::distinct() %>%
tidyr::pivot_wider(
names_from = Escenario,
values_from = valor
)
## 4) Unir todo en una sola tabla
tabla_indicadores_escenarios <- dplyr::bind_rows(
tabla_inc1_wide,
tabla_Reff0_wide,
tabla_Telim_wide
)
tabla_indicadores_escenarios## # A tibble: 15 × 6
## indicador `Combinado fuerte` `Combinado moderado` `PrEP 20%` `Status quo`
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Incidencia m… 0.00117 0.00192 0.00246 0.00250
## 2 R_efectivo,0… 0.0390 0.0633 0.0807 0.0822
## 3 R_efectivo,0… 0.0390 0.0633 0.0807 0.0822
## 4 R_efectivo,0… 0.0390 0.0633 0.0807 0.0822
## 5 R_efectivo,0… 0.0390 0.0633 0.0807 0.0822
## 6 R_efectivo,0… 0.0390 0.0633 0.0807 0.0822
## 7 R_efectivo,0… 0.0390 0.0633 0.0807 0.0822
## 8 R_efectivo,0… 0.0390 0.0633 0.0807 0.0822
## 9 T_elim (α = … 4.1 NA NA NA
## 10 T_elim (α = … 3 13.7 NA NA
## 11 T_elim (α = … 2.3 9.5 NA NA
## 12 T_elim (α = … 1.7 7.2 NA NA
## 13 T_elim (α = … 1.3 5.6 NA NA
## 14 T_elim (α = … 0.7 3.5 NA NA
## 15 T_elim (α = … 0.3 2.1 NA NA
## # ℹ 1 more variable: `TasP 80%` <dbl># Tabla redondeada
tabla_indicadores_escenarios_red <- tabla_indicadores_escenarios %>%
dplyr::mutate(
dplyr::across(
-indicador,
~ ifelse(is.na(.x), NA_real_, round(.x, 3))
)
)
tabla_indicadores_escenarios_red## # A tibble: 15 × 6
## indicador `Combinado fuerte` `Combinado moderado` `PrEP 20%` `Status quo`
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Incidencia m… 0.001 0.002 0.002 0.002
## 2 R_efectivo,0… 0.039 0.063 0.081 0.082
## 3 R_efectivo,0… 0.039 0.063 0.081 0.082
## 4 R_efectivo,0… 0.039 0.063 0.081 0.082
## 5 R_efectivo,0… 0.039 0.063 0.081 0.082
## 6 R_efectivo,0… 0.039 0.063 0.081 0.082
## 7 R_efectivo,0… 0.039 0.063 0.081 0.082
## 8 R_efectivo,0… 0.039 0.063 0.081 0.082
## 9 T_elim (α = … 4.1 NA NA NA
## 10 T_elim (α = … 3 13.7 NA NA
## 11 T_elim (α = … 2.3 9.5 NA NA
## 12 T_elim (α = … 1.7 7.2 NA NA
## 13 T_elim (α = … 1.3 5.6 NA NA
## 14 T_elim (α = … 0.7 3.5 NA NA
## 15 T_elim (α = … 0.3 2.1 NA NA
## # ℹ 1 more variable: `TasP 80%` <dbl>14.3.1 Preparar datos para las gráficas
tabla_Reff_abs_alpha_plot <- tabla_Reff_abs_alpha %>%
dplyr::mutate(
Escenario = factor(
Escenario,
levels = c(
"Status quo",
"TasP 80%",
"PrEP 20%",
"Combinado moderado",
"Combinado fuerte"
)
),
alpha_lab = paste0(alpha_pct, "%") # etiqueta bonita para el eje x
)14.3.2 Heatmap de 𝑅efectivo,0
ggplot(tabla_Reff_abs_alpha_plot,
aes(x = alpha_lab, y = Escenario, fill = R_eff_0)) +
geom_tile(color = "white", linewidth = 0.4) +
geom_text(
aes(label = ifelse(is.na(R_eff_0),
"NA",
sprintf("%.2f", R_eff_0))),
color = "white",
size = 3
) +
scale_fill_viridis_c(
option = "plasma",
name = expression(R[efectivo,0]),
na.value = "grey70"
) +
labs(
x = expression(alpha~"(umbral como % de"~I[1]^{ref}*")"),
y = "Escenario",
title = expression("Número reproductivo efectivo promedio en el año inicial"~R[efectivo,0]),
subtitle = "Valores < 1 indican epidemia en descenso; valores > 1, epidemia en crecimiento"
) +
theme_minimal(base_size = 12) +
theme(
panel.grid = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
plot.title = element_text(face = "bold"),
legend.position = "right",
legend.title = element_text(face = "bold")
)
14.3.3 Heatmap de tiempo hasta eliminación 𝑇 elim T elim
ggplot(tabla_Reff_abs_alpha_plot,
aes(x = alpha_lab, y = Escenario, fill = years_to_elim)) +
geom_tile(color = "white", linewidth = 0.4) +
geom_text(
aes(label = ifelse(is.na(years_to_elim),
"NA",
sprintf("%.1f", years_to_elim))),
color = "white",
size = 3
) +
scale_fill_viridis_c(
option = "magma",
name = "Años hasta eliminación",
na.value = "grey70"
) +
labs(
x = expression(alpha~"(umbral como % de"~I[1]^{ref}*")"),
y = "Escenario",
title = expression("Tiempo hasta eliminación"~T[elim]),
subtitle = "NA indica que no se alcanza el umbral en el horizonte simulado"
) +
theme_minimal(base_size = 12) +
theme(
panel.grid = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
plot.title = element_text(face = "bold"),
legend.position = "right",
legend.title = element_text(face = "bold")
)
15 14. Número reproductivo efectivo por grupo de riesgo (LR y HR) en el escenario base
En el análisis anterior se definió un proxy del número reproductivo efectivo global \(R_{\text{efectivo}}(t)\) agregando todos los compartimentos infectados. Sin embargo, en este modelo la población de HSH está estratificada en dos grupos de riesgo: bajo riesgo (LR) y alto riesgo (HR). Es útil cuantificar, para el escenario Status quo, el número reproductivo efectivo dentro de cada grupo en el año inicial del horizonte de simulación.
15.1 14.1 Definiciones por grupo de riesgo
Denotamos por \(I_{LR}(t)\) e \(I_{HR}(t)\) el número total de infectados en cada grupo de riesgo:
\[ \begin{aligned} I_{LR}(t) &= A1_{LR}(t) + A2_{LR}(t) + A3_{LR}(t) \\ &\quad + C1_{LR}(t) + C2_{LR}(t) + C3_{LR}(t) \\ &\quad + D1_{LR}(t) + D2_{LR}(t) + D3_{LR}(t), \\ \\ I_{HR}(t) &= A1_{HR}(t) + A2_{HR}(t) + A3_{HR}(t) \\ &\quad + C1_{HR}(t) + C2_{HR}(t) + C3_{HR}(t) \\ &\quad + D1_{HR}(t) + D2_{HR}(t) + D3_{HR}(t). \end{aligned} \]
Por construcción del modelo, las nuevas infecciones que pasan a la clase infectada en cada grupo (es decir, flujo de entrada a los compartimentos A1) vienen dadas por las fuerzas de infección en susceptibles:
\[ \begin{aligned} F_{LR}(t) &\equiv \text{new\_inf}_{LR}(t) \\ &= \lambda_{1,LR}(t)\,S_{1,LR}(t) + \lambda_{2,LR}(t)\,S_{2,LR}(t) + \lambda_{3,LR}(t)\,S_{3,LR}(t), \\ \\ F_{HR}(t) &\equiv \text{new\_inf}_{HR}(t) \\ &= \lambda_{1,HR}(t)\,S_{1,HR}(t) + \lambda_{2,HR}(t)\,S_{2,HR}(t) + \lambda_{3,HR}(t)\,S_{3,HR}(t). \end{aligned} \]
Si agregamos todas las ecuaciones de los compartimentos infectados de un grupo \(g \in \{LR, HR\}\), podemos escribir de forma esquemática:
\[ \frac{d I_g(t)}{dt} = F_g(t) - G_g(t), \]
donde:
- \(F_g(t)\): flujo de nuevas infecciones que entran al conjunto de infectados del grupo \(g\).
- \(G_g(t)\): flujo de salida de infectados (progresión a muerte, salida del sistema, etc.) desde el conjunto de infectados del grupo \(g\).
De esta identidad se obtiene:
\[ G_g(t) = F_g(t) - \frac{d I_g(t)}{dt}. \]
Siguiendo la interpretación estándar de modelos compartimentales, definimos el número reproductivo efectivo instantáneo dentro del grupo \(g\) como
\[ R_{\text{efectivo}}^{(g)}(t) = \frac{F_g(t)}{G_g(t)} = \frac{\text{new\_inf}_g(t)} {\text{new\_inf}_g(t) - \dfrac{d I_g(t)}{dt}}. \]
Intuitivamente:
- Si \(d I_g(t)/dt \approx 0\), el grupo \(g\) está aproximadamente en equilibrio y \(R_{\text{efectivo}}^{(g)}(t) \approx 1\).
- Si \(d I_g(t)/dt > 0\), el número de infectados en el grupo \(g\) crece en el tiempo y \(R_{\text{efectivo}}^{(g)}(t) > 1\).
- Si \(d I_g(t)/dt < 0\), la epidemia decrece en ese grupo y \(R_{\text{efectivo}}^{(g)}(t) < 1\).
15.2 14.2 Promedio en el año inicial (\(R_{\text{efectivo},0}^{LR}\) y \(R_{\text{efectivo}}\))
Para evitar la inestabilidad numérica asociada a evaluar en un único instante \(t = 0, se define el número reproductivo efectivo **promedio en el primer año** del horizonte de simulación como:
\[ R_{\text{efectivo},0}^{(g)} \approx \frac{ \displaystyle \sum_{t_i \in [0,1)} R_{\text{efectivo}}^{(g)}(t_i)\,\Delta t_i }{ \displaystyle \sum_{t_i \in [0,1)} \Delta t_i }, \qquad g \in \{LR, HR\}, \]
donde \(t_i\) son los puntos de tiempo de la simulación y \(\Delta t_i\) el paso de integración (en este modelo, \(dt = 0.1\) años).
En el escenario Status quo, los valores \(R_{\text{efectivo},0}^{LR}\) y \(R_{\text{efectivo},0}^{HR}\) resumen el potencial de transmisión inicial en los subgrupos de bajo y alto riesgo, respectivamente:
- Si \(R_{\text{efectivo},0}^{HR} > R_{\text{efectivo},0}^{LR}\), se confirma que la transmisión se sostiene principalmente en el subgrupo de alto riesgo.
- Si ambos son cercanos a 1, la epidemia está aproximadamente en equilibrio en cada subgrupo bajo las condiciones actuales del Status quo
Para calcular estos indicadores, primero se extiende la función
add_indicators_single() para desagregar las nuevas
infecciones en LR y HR. Si ya se ha redefinido esta función de la misma
forma en una sección anterior, esta redefinición es redundante y puede
omitirse; se incluye aquí para que el código quede autocontenido.
# 13.3.1 Extender add_indicators_single para separar LR y HR
add_indicators_single <- function(sim, parms) {
times <- sim$time
n <- nrow(sim)
# Paso de tiempo entre puntos de la simulación
if (n > 1) {
dt_vec <- c(diff(times), tail(diff(times), 1))
} else {
dt_vec <- 0
}
# Vectores para:
# - nuevas infecciones totales
# - nuevas infecciones en LR y HR
# - incidencia instantánea global
new_inf_total <- numeric(n)
new_inf_LR <- numeric(n)
new_inf_HR <- numeric(n)
inc_rate_inst <- numeric(n)
for (i in seq_len(n)) {
# Estado en t_i como vector nombrado
state_vec <- as.numeric(sim[i, state_names])
names(state_vec) <- state_names
# Fuerzas de infección λ en t_i
lambdas <- foi_susceptibles(
t = times[i],
state = state_vec,
parms = parms
)
tmp <- with(as.list(c(parms, as.list(state_vec), lambdas)), {
# Nuevas infecciones en LR (destino LR)
new_inf_LR_i <- lambda1_LR * S1_LR +
lambda2_LR * S2_LR +
lambda3_LR * S3_LR
# Nuevas infecciones en HR (destino HR)
new_inf_HR_i <- lambda1_HR * S1_HR +
lambda2_HR * S2_HR +
lambda3_HR * S3_HR
# Total de nuevas infecciones
new_inf_i <- new_inf_LR_i + new_inf_HR_i
# Incidencia instantánea global
inc_rate_i <- ifelse(sim$N_tot[i] > 0,
new_inf_i / sim$N_tot[i],
NA_real_)
list(
new_inf_total = new_inf_i,
new_inf_LR = new_inf_LR_i,
new_inf_HR = new_inf_HR_i,
inc_rate = inc_rate_i
)
})
new_inf_total[i] <- tmp$new_inf_total
new_inf_LR[i] <- tmp$new_inf_LR
new_inf_HR[i] <- tmp$new_inf_HR
inc_rate_inst[i] <- tmp$inc_rate
}
sim %>%
dplyr::mutate(
dt = dt_vec,
new_inf_year = new_inf_total, # compatibilidad con secciones previas
new_inf_LR = new_inf_LR,
new_inf_HR = new_inf_HR,
inc_rate_inst = inc_rate_inst,
inc_per1000_inst = inc_rate_inst * 1000
)
}# 13.3.2 Cálculo de R_efectivo,0 por grupo LR / HR en el primer año
compute_Reff_by_group_statusquo <- function(parms, init) {
# 1) Simular el modelo con estos parámetros
sim <- sim_vi_hsh(parms = parms, init = init)
sim_i <- add_indicators_single(sim, parms)
time <- sim_i$time
dt_vec <- sim_i$dt
# 2) Infectados totales por grupo de riesgo
I_LR <- sim_i$I_LR
I_HR <- sim_i$I_HR
# 3) Nuevas infecciones por grupo de destino
new_LR <- sim_i$new_inf_LR
new_HR <- sim_i$new_inf_HR
# 4) Derivada numérica dI_g/dt para cada grupo (cociente de diferencias)
dI_LR_dt <- c(
diff(I_LR) / diff(time),
tail(diff(I_LR) / diff(time), 1)
)
dI_HR_dt <- c(
diff(I_HR) / diff(time),
tail(diff(I_HR) / diff(time), 1)
)
# 5) Flujos de salida G_g(t) = F_g(t) - dI_g/dt
removal_LR <- new_LR - dI_LR_dt
removal_HR <- new_HR - dI_HR_dt
# 6) Número reproductivo efectivo instantáneo por grupo
R_eff_LR_t <- ifelse(removal_LR > 0,
new_LR / removal_LR,
NA_real_)
R_eff_HR_t <- ifelse(removal_HR > 0,
new_HR / removal_HR,
NA_real_)
# 7) Promedio en el primer año [0,1)
year1 <- time >= 0 & time < 1
# LR
if (any(year1 & is.finite(R_eff_LR_t))) {
R_eff0_LR <- sum(
R_eff_LR_t[year1] * dt_vec[year1],
na.rm = TRUE
) /
sum(
dt_vec[year1][is.finite(R_eff_LR_t[year1])],
na.rm = TRUE
)
} else {
R_eff0_LR <- NA_real_
}
# HR
if (any(year1 & is.finite(R_eff_HR_t))) {
R_eff0_HR <- sum(
R_eff_HR_t[year1] * dt_vec[year1],
na.rm = TRUE
) /
sum(
dt_vec[year1][is.finite(R_eff_HR_t[year1])],
na.rm = TRUE
)
} else {
R_eff0_HR <- NA_real_
}
list(
R_eff0_LR = R_eff0_LR,
R_eff0_HR = R_eff0_HR
)
}# 13.3.3 Cálculo de R_efectivo,0^LR y R_efectivo,0^HR en el Status quo
res_Reff_groups <- compute_Reff_by_group_statusquo(
parms = parms_base, # o parms_calibrated si se ajustaron los parámetros
init = init
)
tibble::tibble(
grupo_riesgo = c("LR", "HR"),
R_efectivo_0 = c(res_Reff_groups$R_eff0_LR,
res_Reff_groups$R_eff0_HR)
)## # A tibble: 2 × 2
## grupo_riesgo R_efectivo_0
## <chr> <dbl>
## 1 LR 0.0666
## 2 HR 0.11915.2.1 Función para añadir new_inf por grupo y R efectivo (global, LR, HR)
library(dplyr)
add_indicators_Reff_by_group <- function(sim, parms) {
times <- sim$time
n <- nrow(sim)
# dt entre tiempos (último repetido)
if (n > 1) {
dt_vec <- c(diff(times), tail(diff(times), 1))
} else {
dt_vec <- 0
}
new_inf_LR <- numeric(n)
new_inf_HR <- numeric(n)
new_inf_tot <- numeric(n)
# 1) Calcular nuevas infecciones LR/HR en cada tiempo
for (i in seq_len(n)) {
state_vec <- as.numeric(sim[i, state_names])
names(state_vec) <- state_names
lambdas <- foi_susceptibles(
t = times[i],
state = state_vec,
parms = parms
)
tmp <- with(as.list(c(parms, as.list(state_vec), lambdas)), {
# nuevas infecciones por grupo
new_LR <- lambda1_LR * S1_LR + lambda2_LR * S2_LR + lambda3_LR * S3_LR
new_HR <- lambda1_HR * S1_HR + lambda2_HR * S2_HR + lambda3_HR * S3_HR
list(new_LR = new_LR, new_HR = new_HR)
})
new_inf_LR[i] <- tmp$new_LR
new_inf_HR[i] <- tmp$new_HR
new_inf_tot[i] <- tmp$new_LR + tmp$new_HR
}
# 2) Derivadas aproximadas dI/dt (LR, HR, total)
dI_LR <- c(diff(sim$I_LR) / dt_vec[-n], NA_real_)
dI_HR <- c(diff(sim$I_HR) / dt_vec[-n], NA_real_)
dI_tot <- c(diff(sim$I_tot) / dt_vec[-n], NA_real_)
# repetir el último valor para el final
if (n > 1) {
dI_LR[n] <- dI_LR[n - 1]
dI_HR[n] <- dI_HR[n - 1]
dI_tot[n] <- dI_tot[n - 1]
}
# 3) Flujos de salida G = F - dI/dt
G_LR <- new_inf_LR - dI_LR
G_HR <- new_inf_HR - dI_HR
G_tot <- new_inf_tot - dI_tot
# 4) Números reproductivos efectivos instantáneos
R_eff_LR <- ifelse(G_LR > 0, new_inf_LR / G_LR, NA_real_)
R_eff_HR <- ifelse(G_HR > 0, new_inf_HR / G_HR, NA_real_)
R_eff_tot <- ifelse(G_tot > 0, new_inf_tot / G_tot, NA_real_)
# 5) Incidencia instantánea
inc_rate_inst <- ifelse(sim$N_tot > 0, new_inf_tot / sim$N_tot, NA_real_)
inc_per1000_inst <- inc_rate_inst * 1000
sim %>%
mutate(
dt = dt_vec,
new_inf_LR = new_inf_LR,
new_inf_HR = new_inf_HR,
new_inf_tot = new_inf_tot,
dI_LR_dt = dI_LR,
dI_HR_dt = dI_HR,
dI_tot_dt = dI_tot,
R_eff_LR = R_eff_LR,
R_eff_HR = R_eff_HR,
R_eff_tot = R_eff_tot,
inc_rate_inst = inc_rate_inst,
inc_per1000_inst = inc_per1000_inst
)
}
# Simulación escenario base
sim_base <- sim_vi_hsh(parms = parms_base, init = init)
# Añadir indicadores por grupo
sim_base_ind <- add_indicators_Reff_by_group(sim_base, parms_base)15.2.2 Resumen en el primer año: R_efectivo_0 (global, LR, HR) e I1
# Filtrar el primer año
sim_yr1 <- sim_base_ind %>% filter(time < 1)
# Promedios ponderados por dt de R_eff
R0_LR <- with(sim_yr1, sum(R_eff_LR * dt, na.rm = TRUE) / sum(dt))
R0_HR <- with(sim_yr1, sum(R_eff_HR * dt, na.rm = TRUE) / sum(dt))
R0_tot <- with(sim_yr1, sum(R_eff_tot * dt, na.rm = TRUE) / sum(dt))
# Incidencia media año 1 (I1) por grupo y total
I1_LR <- with(sim_yr1,
sum(new_inf_LR * dt, na.rm = TRUE) / sum(N_LR * dt, na.rm = TRUE)
)
I1_HR <- with(sim_yr1,
sum(new_inf_HR * dt, na.rm = TRUE) / sum(N_HR * dt, na.rm = TRUE)
)
I1_tot <- with(sim_yr1,
sum(new_inf_tot * dt, na.rm = TRUE) / sum(N_tot * dt, na.rm = TRUE)
)15.2.3 Tabla comparativa global vs LR vs HR
tabla_Reff0_base <- tibble::tibble(
grupo = c("Global", "LR", "HR"),
R_efectivo_0 = c(R0_tot, R0_LR, R0_HR),
I1 = c(I1_tot, I1_LR, I1_HR)
) %>%
mutate(
R_efectivo_0 = round(R_efectivo_0, 4),
I1 = round(I1, 6)
)
tabla_Reff0_base## # A tibble: 3 × 3
## grupo R_efectivo_0 I1
## <chr> <dbl> <dbl>
## 1 Global 0.0822 0.0025
## 2 LR 0.0666 0.00179
## 3 HR 0.120 0.0053815.2.4 Extender a todos los escenarios:
library(dplyr)
add_Reff_by_group <- function(sim, parms) {
times <- sim$time
n <- nrow(sim)
# dt entre tiempos (repetimos el último para igualar longitud)
if (n > 1) {
dt_vec <- c(diff(times), tail(diff(times), 1))
} else {
dt_vec <- 0
}
new_inf_LR <- numeric(n)
new_inf_HR <- numeric(n)
new_inf_tot <- numeric(n)
for (i in seq_len(n)) {
# estado en t_i
state_vec <- as.numeric(sim[i, state_names])
names(state_vec) <- state_names
lambdas <- foi_susceptibles(
t = times[i],
state = state_vec,
parms = parms
)
tmp <- with(as.list(c(parms, as.list(state_vec), lambdas)), {
# nuevas infecciones por grupo
new_LR <- lambda1_LR * S1_LR + lambda2_LR * S2_LR + lambda3_LR * S3_LR
new_HR <- lambda1_HR * S1_HR + lambda2_HR * S2_HR + lambda3_HR * S3_HR
list(new_LR = new_LR, new_HR = new_HR)
})
new_inf_LR[i] <- tmp$new_LR
new_inf_HR[i] <- tmp$new_HR
new_inf_tot[i] <- tmp$new_LR + tmp$new_HR
}
# derivadas aproximadas dI/dt
dI_LR <- c(diff(sim$I_LR) / dt_vec[-n], NA_real_)
dI_HR <- c(diff(sim$I_HR) / dt_vec[-n], NA_real_)
dI_tot <- c(diff(sim$I_tot) / dt_vec[-n], NA_real_)
if (n > 1) {
dI_LR[n] <- dI_LR[n - 1]
dI_HR[n] <- dI_HR[n - 1]
dI_tot[n] <- dI_tot[n - 1]
}
# flujos de salida
G_LR <- new_inf_LR - dI_LR
G_HR <- new_inf_HR - dI_HR
G_tot <- new_inf_tot - dI_tot
# números reproductivos efectivos instantáneos
R_eff_LR <- ifelse(G_LR > 0, new_inf_LR / G_LR, NA_real_)
R_eff_HR <- ifelse(G_HR > 0, new_inf_HR / G_HR, NA_real_)
R_eff_tot <- ifelse(G_tot > 0, new_inf_tot / G_tot, NA_real_)
# incidencia instantánea (global)
inc_rate_inst <- ifelse(sim$N_tot > 0, new_inf_tot / sim$N_tot, NA_real_)
inc_per1000_inst <- inc_rate_inst * 1000
sim %>%
mutate(
dt = dt_vec,
new_inf_LR = new_inf_LR,
new_inf_HR = new_inf_HR,
new_inf_tot = new_inf_tot,
dI_LR_dt = dI_LR,
dI_HR_dt = dI_HR,
dI_tot_dt = dI_tot,
R_eff_LR = R_eff_LR,
R_eff_HR = R_eff_HR,
R_eff_tot = R_eff_tot,
inc_rate_inst = inc_rate_inst,
inc_per1000_inst = inc_per1000_inst
)
}
resumen_Reff_I1_yr1 <- function(sim_with_Reff) {
sim_yr1 <- sim_with_Reff %>% filter(time < 1)
# Promedios ponderados de R_eff
R0_LR <- with(sim_yr1, sum(R_eff_LR * dt, na.rm = TRUE) / sum(dt))
R0_HR <- with(sim_yr1, sum(R_eff_HR * dt, na.rm = TRUE) / sum(dt))
R0_tot <- with(sim_yr1, sum(R_eff_tot * dt, na.rm = TRUE) / sum(dt))
# Incidencias medias
I1_LR <- with(sim_yr1,
sum(new_inf_LR * dt, na.rm = TRUE) / sum(N_LR * dt, na.rm = TRUE)
)
I1_HR <- with(sim_yr1,
sum(new_inf_HR * dt, na.rm = TRUE) / sum(N_HR * dt, na.rm = TRUE)
)
I1_tot <- with(sim_yr1,
sum(new_inf_tot * dt, na.rm = TRUE) / sum(N_tot * dt, na.rm = TRUE)
)
tibble::tibble(
grupo = c("Global", "LR", "HR"),
R_efectivo_0 = c(R0_tot, R0_LR, R0_HR),
I1 = c(I1_tot, I1_LR, I1_HR)
)
}
escenarios <- c("Status quo",
"TasP 80%",
"PrEP 20%",
"Combinado moderado",
"Combinado fuerte")
mu_val <- parms_base$mu
rho_levels <- c(
0,
mu_val,
1.5 * mu_val,
2.0 * mu_val
)
library(purrr)
library(tidyr)
tabla_Reff_todos <- purrr::map_dfr(escenarios, function(sc) {
purrr::map_dfr(rho_levels, function(rho) {
parms_sc <- make_parms_scenario(parms_base, sc)
parms_sc$rho_LR <- rho
parms_sc$rho_HR <- rho
sim <- sim_vi_hsh(parms = parms_sc, init = init)
simR <- add_Reff_by_group(sim, parms_sc)
resumen <- resumen_Reff_I1_yr1(simR)
resumen %>%
mutate(
escenario = sc,
rho = rho
)
})
})
tabla_Reff_todos_limpia <- tabla_Reff_todos %>%
mutate(
R_efectivo_0 = round(R_efectivo_0, 4),
I1 = round(I1, 6),
rho = round(rho, 5),
escenario = factor(
escenario,
levels = c("Status quo",
"TasP 80%",
"PrEP 20%",
"Combinado moderado",
"Combinado fuerte")
)
) %>%
arrange(rho, escenario, match(grupo, c("Global", "LR", "HR")))
tabla_Reff_todos_limpia## # A tibble: 60 × 5
## grupo R_efectivo_0 I1 escenario rho
## <chr> <dbl> <dbl> <fct> <dbl>
## 1 Global 0.0822 0.0025 Status quo 0
## 2 LR 0.0666 0.00179 Status quo 0
## 3 HR 0.120 0.00538 Status quo 0
## 4 Global 0.0791 0.00236 TasP 80% 0
## 5 LR 0.0641 0.00169 TasP 80% 0
## 6 HR 0.115 0.00509 TasP 80% 0
## 7 Global 0.0807 0.00246 PrEP 20% 0
## 8 LR 0.0655 0.00176 PrEP 20% 0
## 9 HR 0.117 0.00529 PrEP 20% 0
## 10 Global 0.0633 0.00192 Combinado moderado 0
## # ℹ 50 more rowslibrary(ggplot2)
library(dplyr)
plot_heat_R0 <- tabla_Reff_todos_limpia %>%
mutate(
escenario = factor(
escenario,
levels = c("Status quo",
"TasP 80%",
"PrEP 20%",
"Combinado moderado",
"Combinado fuerte")
),
grupo = factor(grupo, levels = c("Global", "LR", "HR")),
rho_lab = paste0("rho = ", round(rho, 5))
)
ggplot(plot_heat_R0,
aes(x = escenario, y = grupo, fill = R_efectivo_0)) +
geom_tile(color = "white") +
scale_fill_viridis_c(
option = "C",
direction = -1,
name = expression(R[efectivo[0]]),
na.value = "grey90"
) +
facet_wrap(~ rho_lab) +
labs(
x = "Escenario",
y = "Grupo de riesgo",
title = expression("N\u00FAmero reproductivo efectivo medio en el a\u00F1o 0"),
subtitle = "Por escenario, grupo de riesgo y nivel de reposición demográfica (rho)"
) +
theme_minimal(base_size = 11) +
theme(
legend.position = "right",
axis.text.x = element_text(angle = 45, hjust = 1),
strip.text = element_text(face = "bold"),
panel.grid = element_blank()
)
16 15. Número reproductivo básico \(R_0\) mediante matriz de próxima generación (Es un indicador de gravedad basal de la epidemia)
En modelos compartimentales, el número reproductivo básico \(R_0\) se define como el número esperado de infecciones secundarias producidas por un individuo infeccioso típico introducido en una población completamente susceptible, en ausencia de intervenciones adicionales.
Para derivar \(R_0\) usamos el enfoque de matriz de próxima generación simplificado a dos grupos de riesgo:
- \(g \in \{\text{LR}, \text{HR}\}\): bajo y alto riesgo.
- Cada individuo infectado pasa por las fases aguda, crónica y SIDA,
pero para \(R_0\) se colapsan en una
“trayectoria de infección típica” con:
- una duración total de infección \(D\);
- una infectividad promedio por pareja \(\beta_g\) (probabilidad de transmisión por pareja a lo largo de toda la infección).
16.1 13.1. Duración total de la infección y pesos por fase
En el modelo, las duraciones medias de cada fase (en años) son:
- Fase aguda: \(\text{dur\_acute\_y}\),
- Fase crónica: \(\text{dur\_chronic\_y}\),
- Fase SIDA: \(\text{dur\_aids\_y}\).
La duración total de la infección se aproxima como
\[ D = \text{dur\_acute\_y} + \text{dur\_chronic\_y} + \text{dur\_aids\_y}. \]
Definimos los pesos de tiempo relativo en cada fase:
\[ w_A = \frac{\text{dur\_acute\_y}}{D}, \qquad w_C = \frac{\text{dur\_chronic\_y}}{D}, \qquad w_D = \frac{\text{dur\_aids\_y}}{D}, \]
de manera que \(w_A + w_C + w_D = 1\).
16.2 13.2. Infectividad promedio por pareja \(\beta_g\)
Sea:
- \(\sigma_{\text{acute}}\): probabilidad de transmisión por pareja en fase aguda;
- \(\sigma_{\text{chronic}}\): probabilidad de transmisión por pareja en fase crónica;
- \(\sigma_{\text{aids}}\): probabilidad de transmisión por pareja en SIDA;
- \(u\): proporción de actos con condón (\(\text{condom\_use}\));
- \(\varepsilon\): eficacia del condón (\(\text{condom\_eff}\)).
El factor de protección global por condón es:
\[ \text{condom\_factor} = 1 - u\,\varepsilon. \]
Si ignoramos, para el cálculo de \(R_0\), efectos de diagnóstico y TAR (es decir, consideramos una infección “sin control” durante toda su trayectoria), la infectividad promedio por pareja de un individuo típico es:
\[ \beta = \text{condom\_factor}\, \big( w_A \,\sigma_{\text{acute}} + w_C \,\sigma_{\text{chronic}} + w_D \,\sigma_{\text{aids}} \big). \] En este esquema simplificado tomamos \(\beta_{\text{LR}} = \beta_{\text{HR}} = \beta\), y dejamos que la diferencia de riesgo entre LR y HR venga dada por:
- el número de parejas/año \(c_{\text{LR}}, c_{\text{HR}}\),
- la mezcla sexual \(H\) entre grupos.
16.3 13.3. Matriz de mezcla sexual \(H\)
Como antes, definimos:
- \(N_{\text{LR}}\), \(N_{\text{HR}}\): tamaños de los grupos LR y HR (a partir de men_18_59, msm_frac y lr_frac).
- \(H\): matriz de mezcla 2×2: \[ H = \begin{pmatrix} H_{\text{LR,LR}} & H_{\text{LR,HR}} \\ H_{\text{HR,LR}} & H_{\text{HR,HR}} \end{pmatrix}, \] donde, aproximadamente, \(H_{g,g'}\) es la probabilidad de que una pareja de un individuo del grupo \(g\) provenga del grupo \(g'\).
La función mixing_matrix() ya implementa esta estructura
y garantiza reciprocidad de contactos.
16.4 13.4. Matriz de próxima generación \(K\) (2×2)
Denotemos por:
- \(c_{\text{LR}}\), \(c_{\text{HR}}\): número de parejas por año
de un individuo de bajo y alto riesgo (
partners_year["LR"]y["HR"]); - \(\beta_{\text{LR}}, \beta_{\text{HR}}\): infectividad promedio por pareja en cada grupo (aquí las tomamos iguales: \(\beta_{\text{LR}}=\beta_{\text{HR}}=\beta\));
- \(D\): duración total de la infección.
El elemento \(K_{i,j}\) de la matriz de próxima generación se interpreta como el número esperado de infecciones en el grupo \(i\) causadas, a lo largo de toda su infección, por un individuo infectado en el grupo \(j\):
\[ K_{i,j} = c_j \, H_{i,j} \, \beta_j \, D, \qquad i,j \in \{\text{LR}, \text{HR}\}. \]
Específicamente:
\[ \begin{aligned} K_{\text{LR,LR}} &= c_{\text{LR}}\, H_{\text{LR,LR}}\, \beta_{\text{LR}}\, D,\\[4pt] K_{\text{LR,HR}} &= c_{\text{HR}}\, H_{\text{LR,HR}}\, \beta_{\text{HR}}\, D,\\[4pt] K_{\text{HR,LR}} &= c_{\text{LR}}\, H_{\text{HR,LR}}\, \beta_{\text{LR}}\, D,\\[4pt] K_{\text{HR,HR}} &= c_{\text{HR}}\, H_{\text{HR,HR}}\, \beta_{\text{HR}}\, D. \end{aligned} \]
La matriz de próxima generación resulta:
\[ K = \begin{pmatrix} K_{\text{LR,LR}} & K_{\text{LR,HR}} \\ K_{\text{HR,LR}} & K_{\text{HR,HR}} \end{pmatrix}. \]
El número reproductivo básico global es el radio espectral de \(K\), es decir, el valor propio de mayor módulo:
\[ R_0 = \rho(K). \]
En una matriz 2×2 con entradas \(a=K_{\text{LR,LR}}, d=K_{\text{HR,HR}}, b=K_{\text{LR,HR}}, c=K_{\text{HR,LR}}\), los valores propios son:
\[ \lambda_{\pm} = \frac{(a + d) \pm \sqrt{(a-d)^2 + 4bc}}{2}, \]
y entonces
\[ R_0 = \max(\lambda_+, \lambda_-). \]
Además, podemos definir números reproductivos por grupo de origen como:
- \(R_0^{\text{LR}} = K_{\text{LR,LR}} + K_{\text{HR,LR}}\): infecciones (LR + HR) generadas por un infectado LR;
- \(R_0^{\text{HR}} = K_{\text{LR,HR}} + K_{\text{HR,HR}}\): infecciones (LR + HR) generadas por un infectado HR.
Estos no son “\(R_0\) distintos”, sino componentes de la matriz que ayudan a entender cuál grupo contribuye más a la transmisión.
#Función general: calcular matriz K y R0 a partir de parms
compute_R0_components <- function(parms) {
with(as.list(parms), {
# 1) Duración total de la infección y pesos por fase
D <- dur_acute_y + dur_chronic_y + dur_aids_y
wA <- dur_acute_y / D
wC <- dur_chronic_y / D
wD <- dur_aids_y / D
# 2) Factor de protección por condón
condom_factor <- 1 - condom_use * condom_eff
if (!is.finite(condom_factor) || condom_factor < 0) condom_factor <- 0
# 3) Infectividad promedio por pareja (ignorando dx/TAR)
beta <- condom_factor * (
wA * sigma_acute +
wC * sigma_chronic +
wD * sigma_aids
)
beta_LR <- beta
beta_HR <- beta
# 4) Tamaños de grupo y matriz de mezcla H
N_hsh <- men_18_59 * msm_frac
N_LR <- N_hsh * lr_frac
N_HR <- N_hsh - N_LR
H <- mixing_matrix(N_LR, N_HR, mixing_m)
# 5) Nº de parejas/año por grupo
c_LR <- partners_year["LR"]
c_HR <- partners_year["HR"]
# 6) Elementos de K (2x2)
K_LR_LR <- c_LR * H["LR", "LR"] * beta_LR * D
K_LR_HR <- c_HR * H["LR", "HR"] * beta_HR * D
K_HR_LR <- c_LR * H["HR", "LR"] * beta_LR * D
K_HR_HR <- c_HR * H["HR", "HR"] * beta_HR * D
K <- matrix(
c(K_LR_LR, K_HR_LR,
K_LR_HR, K_HR_HR),
nrow = 2, byrow = FALSE,
dimnames = list(c("LR", "HR"), c("LR", "HR"))
)
# 7) Autovalores y R0
eig <- eigen(K)
eig_vals <- Re(eig$values)
R0 <- max(eig_vals)
# 8) R0 por grupo de origen
R0_LR <- K["LR", "LR"] + K["HR", "LR"]
R0_HR <- K["LR", "HR"] + K["HR", "HR"]
contrib_LR <- R0_LR / (R0_LR + R0_HR)
contrib_HR <- R0_HR / (R0_LR + R0_HR)
list(
K = K,
R0 = R0,
eigen = eig,
R0_LR = R0_LR,
R0_HR = R0_HR,
contrib = c(LR = contrib_LR, HR = contrib_HR)
)
})
}
#Aplicar al escenario base (Status quo)
R0_base <- compute_R0_components(parms_base)
R0_base$R0 # R0 global## [1] 0.2852254R0_base$R0_LR # infecciones generadas por un infectado LR## [1] 0.2261272R0_base$R0_HR # infecciones generadas por un infectado HR## [1] 0.3155871R0_base$contrib # proporción de R0 atribuible a LR y HR## LR HR
## 0.4174289 0.5825711R0_base$K # matriz K completa## LR HR
## LR 0.13695027 0.07617618
## HR 0.08917692 0.23941087#Extender a TODOS los escenarios
make_parms_scenario <- function(base_parms, scenario) {
p <- base_parms
p$scenario <- scenario
# 1) Escenario base: Status quo
if (scenario == "Status quo") {
return(p)
}
# 2) TasP 80%
if (scenario == "TasP 80%") {
p$tested_12m_prop <- c(LR = 0.80, HR = 0.80)
p$psi1_LR <- derive_psi_from_prop(p$tested_12m_prop["LR"])
p$psi1_HR <- derive_psi_from_prop(p$tested_12m_prop["HR"])
p$psi4_LR <- p$psi4_LR * 1.5
p$psi4_HR <- p$psi4_HR * 1.5
p$psi7_LR <- p$psi7_LR * 1.5
p$psi7_HR <- p$psi7_HR * 1.5
p$psi10_LR <- p$psi10_LR * 1.5
p$psi10_HR <- p$psi10_HR * 1.5
p$alpha5 <- p$alpha5 * 1.5
p$alpha8 <- p$alpha8 * 1.5
p$alpha11 <- p$alpha11 * 1.5
p$art_cov_target <- 0.95
return(p)
}
# 3) PrEP 20%
if (scenario == "PrEP 20%") {
p$prep_cov_target <- 0.20
p$phi2 <- 0.25
p$phi3 <- 0.0
p$prep_adherence <- min(0.85, p$prep_adherence * 0.95)
return(p)
}
# 4) Combinado moderado
if (scenario == "Combinado moderado") {
p$tested_12m_prop <- c(LR = 0.70, HR = 0.75)
p$psi1_LR <- derive_psi_from_prop(p$tested_12m_prop["LR"])
p$psi1_HR <- derive_psi_from_prop(p$tested_12m_prop["HR"])
p$psi4_LR <- p$psi4_LR * 1.3
p$psi4_HR <- p$psi4_HR * 1.3
p$psi7_LR <- p$psi7_LR * 1.3
p$psi7_HR <- p$psi7_HR * 1.3
p$psi10_LR <- p$psi10_LR * 1.3
p$psi10_HR <- p$psi10_HR * 1.3
p$alpha5 <- p$alpha5 * 1.3
p$alpha8 <- p$alpha8 * 1.3
p$alpha11 <- p$alpha11 * 1.3
p$art_cov_target <- 0.93
p$prep_cov_target <- 0.20
p$phi2 <- 0.20
p$phi3 <- 0.0
p$prep_adherence <- min(0.80, p$prep_adherence * 0.90)
p$condom_use <- min(0.80, p$condom_use + 0.15)
return(p)
}
# 5) Combinado fuerte
if (scenario == "Combinado fuerte") {
p$tested_12m_prop <- c(LR = 0.90, HR = 0.90)
p$psi1_LR <- derive_psi_from_prop(p$tested_12m_prop["LR"])
p$psi1_HR <- derive_psi_from_prop(p$tested_12m_prop["HR"])
p$psi4_LR <- p$psi4_LR * 2.0
p$psi4_HR <- p$psi4_HR * 2.0
p$psi7_LR <- p$psi7_LR * 2.0
p$psi7_HR <- p$psi7_HR * 2.0
p$psi10_LR <- p$psi10_LR * 2.0
p$psi10_HR <- p$psi10_HR * 2.0
p$alpha5 <- p$alpha5 * 1.8
p$alpha8 <- p$alpha8 * 1.8
p$alpha11 <- p$alpha11 * 1.8
p$art_cov_target <- 0.98
p$prep_cov_target <- 0.30
p$phi2 <- 0.35
p$phi3 <- 0.0
p$prep_adherence <- min(0.90, p$prep_adherence * 0.95)
p$condom_use <- min(0.95, p$condom_use + 0.25)
p$partners_year["HR"] <- p$partners_year["HR"] * 0.7
p$beh_red_dx <- min(0.50, p$beh_red_dx + 0.20)
p$beh_red_art <- min(0.70, p$beh_red_art + 0.20)
p$beh_red_aids <- min(0.95, p$beh_red_aids + 0.05)
return(p)
}
stop("Escenario no reconocido: ", scenario)
}
library(dplyr)
library(purrr)
library(tidyr)
tabla_R0_esc <- map_dfr(escenarios, function(sc) {
parms_sc <- make_parms_scenario(parms_base, sc)
# Para R0, normalmente ponemos rho = 0 (cohorte cerrada)
parms_sc$rho_LR <- 0
parms_sc$rho_HR <- 0
R0_sc <- compute_R0_components(parms_sc)
tibble(
escenario = sc,
R0 = R0_sc$R0,
R0_LR = R0_sc$R0_LR,
R0_HR = R0_sc$R0_HR,
contrib_LR = R0_sc$contrib["LR"],
contrib_HR = R0_sc$contrib["HR"]
)
}) %>%
mutate(
R0 = round(R0, 3),
R0_LR = round(R0_LR, 3),
R0_HR = round(R0_HR, 3),
contrib_LR = round(contrib_LR, 3),
contrib_HR = round(contrib_HR, 3),
escenario = factor(
escenario,
levels = c("Status quo", "TasP 80%", "PrEP 20%",
"Combinado moderado", "Combinado fuerte")
)
)
tabla_R0_esc## # A tibble: 5 × 6
## escenario R0 R0_LR R0_HR contrib_LR contrib_HR
## <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Status quo 0.285 0.226 0.316 0.417 0.583
## 2 TasP 80% 0.285 0.226 0.316 0.417 0.583
## 3 PrEP 20% 0.285 0.226 0.316 0.417 0.583
## 4 Combinado moderado 0.237 0.188 0.262 0.417 0.583
## 5 Combinado fuerte 0.16 0.162 0.158 0.506 0.49416.4.1 Heatmap de la matriz K por escenario
library(dplyr)
library(purrr)
library(tibble)
tabla_K_esc <- map_dfr(escenarios, function(sc) {
parms_sc <- make_parms_scenario(parms_base, sc)
parms_sc$rho_LR <- 0
parms_sc$rho_HR <- 0
K_sc <- compute_R0_components(parms_sc)$K
# Paso estable: data.frame -> tibble
dfK <- as.data.frame(as.table(K_sc), stringsAsFactors = FALSE)
# Ahora sí, renombrar explícitamente
colnames(dfK) <- c("destino", "origen", "valor")
as_tibble(dfK) %>%
mutate(escenario = sc)
})
tabla_K_esc_limpia <- tabla_K_esc %>%
mutate(
escenario = factor(
escenario,
levels = c("Status quo", "TasP 80%", "PrEP 20%",
"Combinado moderado", "Combinado fuerte")
),
origen = factor(origen, levels = c("LR", "HR")),
destino = factor(destino, levels = c("LR", "HR"))
)
ggplot(tabla_K_esc_limpia,
aes(x = origen, y = destino, fill = valor)) +
geom_tile(color = "white") +
scale_fill_viridis_c(
option = "C",
direction = -1,
name = "Infecciones esperadas"
) +
facet_wrap(~ escenario) +
labs(
x = "Grupo de origen (infectado)",
y = "Grupo de destino (nuevas infecciones)",
title = "Matriz de próxima generación K por escenario",
subtitle = "Número esperado de infecciones a lo largo de toda la infección de un individuo"
) +
theme_minimal(base_size = 11) +
theme(
panel.grid = element_blank(),
legend.position = "right",
strip.text = element_text(face = "bold")
)
16.4.2 Diagrama de Sankey de flujos LR/HR
library(networkD3)
# Tomamos, por ejemplo, el escenario Status quo
K_sq <- compute_R0_components(
make_parms_scenario(parms_base, "Status quo")
)$K
# Preparamos nodos y enlaces
nodes <- data.frame(
name = c("Inf_LR", "Inf_HR", "Inf_LR→Sus_LR", "Inf_LR→Sus_HR",
"Inf_HR→Sus_LR", "Inf_HR→Sus_HR")
)
links <- data.frame(
source = c(0, 0, 1, 1), # índices de nodos origen
target = c(2, 3, 4, 5), # índices de nodos destino
value = c(K_sq["LR","LR"], K_sq["LR","HR"],
K_sq["HR","LR"], K_sq["HR","HR"])
)
# Sankey
sankeyNetwork(
Links = links,
Nodes = nodes,
Source = "source",
Target = "target",
Value = "value",
NodeID = "name",
fontSize = 12,
nodeWidth = 30
)library(networkD3)
make_sankey_K <- function(scenario, base_parms) {
# 1. Parámetros del escenario
parms_sc <- make_parms_scenario(base_parms, scenario)
# Para R0 / K asumimos cohorte sin reposición
parms_sc$rho_LR <- 0
parms_sc$rho_HR <- 0
# 2. Matriz K del escenario
K_sc <- compute_R0_components(parms_sc)$K
# 3. Nodos
nodes <- data.frame(
name = c("Inf_LR", "Inf_HR",
"Inf_LR → Sus_LR", "Inf_LR → Sus_HR",
"Inf_HR → Sus_LR", "Inf_HR → Sus_HR")
)
# 4. Enlaces (links)
links <- data.frame(
source = c(0, 0, 1, 1), # índices de nodos origen
target = c(2, 3, 4, 5), # índices de nodos destino
value = c(
K_sc["LR", "LR"], # Inf_LR → Sus_LR
K_sc["LR", "HR"], # Inf_LR → Sus_HR
K_sc["HR", "LR"], # Inf_HR → Sus_LR
K_sc["HR", "HR"] # Inf_HR → Sus_HR
)
)
# 5. Construir Sankey
sankeyNetwork(
Links = links,
Nodes = nodes,
Source = "source",
Target = "target",
Value = "value",
NodeID = "name",
fontSize = 14,
nodeWidth = 30,
sinksRight = TRUE
)
}
escenarios <- c("Status quo",
"TasP 80%",
"PrEP 20%",
"Combinado moderado",
"Combinado fuerte")
sankey_list <- lapply(escenarios, make_sankey_K, base_parms = parms_base)
names(sankey_list) <- escenarios
sankey_list[["Status quo"]]sankey_list[["TasP 80%"]]sankey_list[["PrEP 20%"]]sankey_list[["Combinado moderado"]]sankey_list[["Combinado fuerte"]]17 16. Intervenciones conductuales aisladas: parejas HR y uso de condón
En este apartado se aborda el Objetivo 3 del estudio:
Evaluar el impacto de la implementación de intervenciones conductuales de forma aislada (reducción del número anual de parejas sexuales en HSH de alto riesgo o aumento de la tasa general de uso de condón) para determinar la intensidad mínima requerida de cada medida para lograr la eliminación de la transmisión del VIH para 2050.
Se consideran dos parámetros conductuales clave, modificados uno por uno sobre el escenario de Status quo calibrado:
- \(c_{\text{HR}}\): número de
parejas por año en HSH de alto riesgo
(
partners_year["HR"]). - \(u\): proporción de actos sexuales
protegidos con condón (
condom_use).
Como criterio de eliminación se usa el mismo que en la sección anterior:
Umbral de incidencia instantánea: \[ \tau = \alpha \, I_1^{\text{ref}}, \] donde \(I_1^{\text{ref}}\) es la incidencia media del primer año en el Status quo calibrado y \(\alpha\) es una fracción (por defecto \(\alpha = 0{,}1\)).
Tiempo hasta eliminación: \[ T_{\text{elim}}(\theta) = \min\left\{ t \ge 0 : \text{inc}(t;\theta) < \tau \right\}. \]
En este análisis se consideran solo cambios conductuales, manteniendo:
- Estrategia de prueba y tratamiento (TasP) del Status quo calibrado.
- Cobertura de PrEP en cero (o en el valor base).
- Resto de parámetros iguales a
parms_baseya calibrados.
# Fracción alpha para el umbral de eliminación (10% por defecto)
alpha_behav <- 0.10
# Horizonte de simulación en años (ya definido en parms_base)
horiz_sim <- parms_base$t_end - parms_base$t_start # por ej. 40 años
# Secuencia de posibles números de parejas en HR
partners_seq <- seq(4, 20, by = 1) # puedes ajustar el rango
# Secuencia de posibles niveles de uso de condón
condom_seq <- seq(0.10, 0.95, by = 0.05) # de 10% a 95%
library(dplyr)
library(tidyr)
library(purrr)
library(ggplot2)
library(viridis)## Loading required package: viridisLite##
## Attaching package: 'viridis'## The following object is masked from 'package:scales':
##
## viridis_pal17.1 16.1 Análisis 1D (parejas HR y condón) ajustado
## 14.1 Análisis 1D: solo parejas HR o solo uso de condón
# (a) Variación del número de parejas en HRMSM (c_HR),
# manteniendo el resto de parámetros en Status quo calibrado
sens_partners <- map_dfr(partners_seq, function(n_hr) {
p <- parms_base # Status quo calibrado (ya con sigmas calibradas)
p$partners_year["HR"] <- n_hr
stats <- compute_R_eff_abs_elim(
parms = p,
init = init,
inc1_ref = inc1_ref,
alpha = alpha_behav
)
tibble(
tipo = "Parejas_HR",
partners_hr = n_hr,
condom_use = p$condom_use,
R_eff_0 = stats$R_eff_0,
years_to_elim = stats$years_to_elim,
elim_dentro_hz = !is.na(stats$years_to_elim) &&
stats$years_to_elim <= horiz_sim
)
})
# (b) Variación del uso de condón, manteniendo c_HR en el valor calibrado
sens_condom <- map_dfr(condom_seq, function(cu) {
p <- parms_base # Status quo calibrado
p$condom_use <- cu
stats <- compute_R_eff_abs_elim(
parms = p,
init = init,
inc1_ref = inc1_ref,
alpha = alpha_behav
)
tibble(
tipo = "Uso_condon",
partners_hr = p$partners_year["HR"],
condom_use = cu,
R_eff_0 = stats$R_eff_0,
years_to_elim = stats$years_to_elim,
elim_dentro_hz = !is.na(stats$years_to_elim) &&
stats$years_to_elim <= horiz_sim
)
})
sens_partners## # A tibble: 17 × 6
## tipo partners_hr condom_use R_eff_0 years_to_elim elim_dentro_hz
## <chr> <dbl> <dbl> <dbl> <dbl> <lgl>
## 1 Parejas_HR 4 0.17 0.0547 NA FALSE
## 2 Parejas_HR 5 0.17 0.0570 NA FALSE
## 3 Parejas_HR 6 0.17 0.0594 NA FALSE
## 4 Parejas_HR 7 0.17 0.0618 NA FALSE
## 5 Parejas_HR 8 0.17 0.0643 NA FALSE
## 6 Parejas_HR 9 0.17 0.0667 NA FALSE
## 7 Parejas_HR 10 0.17 0.0692 NA FALSE
## 8 Parejas_HR 11 0.17 0.0718 NA FALSE
## 9 Parejas_HR 12 0.17 0.0743 NA FALSE
## 10 Parejas_HR 13 0.17 0.0769 NA FALSE
## 11 Parejas_HR 14 0.17 0.0795 NA FALSE
## 12 Parejas_HR 15 0.17 0.0822 NA FALSE
## 13 Parejas_HR 16 0.17 0.0849 NA FALSE
## 14 Parejas_HR 17 0.17 0.0876 NA FALSE
## 15 Parejas_HR 18 0.17 0.0903 NA FALSE
## 16 Parejas_HR 19 0.17 0.0931 NA FALSE
## 17 Parejas_HR 20 0.17 0.0959 NA FALSEsens_condom## # A tibble: 18 × 6
## tipo partners_hr condom_use R_eff_0 years_to_elim elim_dentro_hz
## <chr> <dbl> <dbl> <dbl> <dbl> <lgl>
## 1 Uso_condon 15 0.1 0.0898 NA FALSE
## 2 Uso_condon 15 0.15 0.0844 NA FALSE
## 3 Uso_condon 15 0.2 0.0790 NA FALSE
## 4 Uso_condon 15 0.25 0.0737 NA FALSE
## 5 Uso_condon 15 0.3 0.0685 NA FALSE
## 6 Uso_condon 15 0.35 0.0634 NA FALSE
## 7 Uso_condon 15 0.4 0.0584 NA FALSE
## 8 Uso_condon 15 0.45 0.0534 NA FALSE
## 9 Uso_condon 15 0.5 0.0486 NA FALSE
## 10 Uso_condon 15 0.55 0.0438 NA FALSE
## 11 Uso_condon 15 0.6 0.0391 NA FALSE
## 12 Uso_condon 15 0.65 0.0345 NA FALSE
## 13 Uso_condon 15 0.7 0.0299 NA FALSE
## 14 Uso_condon 15 0.75 0.0255 NA FALSE
## 15 Uso_condon 15 0.8 0.0211 NA FALSE
## 16 Uso_condon 15 0.85 0.0168 3.8 TRUE
## 17 Uso_condon 15 0.9 0.0125 1.7 TRUE
## 18 Uso_condon 15 0.95 0.00835 0.5 TRUE17.2 16.2 Gráficos 1D ajustados al horizonte simulado
# (1) R_efectivo,0 vs número de parejas en HR
ggplot(sens_partners,
aes(x = partners_hr, y = R_eff_0)) +
geom_line(linewidth = 1.2, color = "#440154FF") +
geom_point(aes(color = elim_dentro_hz),
size = 2) +
scale_color_manual(
values = c(`TRUE` = "#2ECC71", `FALSE` = "#E74C3C"),
name = "Eliminación\ndentro del horizonte"
) +
geom_hline(yintercept = 1, linetype = "dashed", color = "grey50") +
labs(
x = "Número de parejas/año en HRMSM",
y = expression(R[efectivo,0]),
title = expression(R[efectivo,0]~
"según número de parejas en HRMSM"),
subtitle = "Línea punteada: umbral R[efectivo,0] == 1"
) +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold"))
# (2) Años hasta eliminación vs número de parejas en HR
ggplot(sens_partners,
aes(x = partners_hr, y = years_to_elim)) +
geom_line(linewidth = 1.2, color = "#31688EFF") +
geom_point(aes(color = elim_dentro_hz),
size = 2) +
scale_color_manual(
values = c(`TRUE` = "#2ECC71", `FALSE` = "#E74C3C"),
name = "Eliminación\ndentro del horizonte"
) +
geom_hline(yintercept = horiz_sim,
linetype = "dashed", color = "red") +
labs(
x = "Número de parejas/año en HRMSM",
y = "Años hasta eliminación (T_elim)",
title = "Tiempo hasta eliminación según reducción de parejas en HRMSM",
subtitle = paste0("Línea roja: horizonte de simulación = ",
horiz_sim, " años")
) +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold"))## Warning: Removed 17 rows containing missing values or values outside the scale range
## (`geom_line()`).## Warning: Removed 17 rows containing missing values or values outside the scale range
## (`geom_point()`).
# (3) R_efectivo,0 vs uso de condón
ggplot(sens_condom,
aes(x = condom_use * 100, y = R_eff_0)) +
geom_line(linewidth = 1.2, color = "#440154FF") +
geom_point(aes(color = elim_dentro_hz),
size = 2) +
scale_color_manual(
values = c(`TRUE` = "#2ECC71", `FALSE` = "#E74C3C"),
name = "Eliminación\ndentro del horizonte"
) +
geom_hline(yintercept = 1, linetype = "dashed", color = "grey50") +
labs(
x = "Uso de condón en actos sexuales (%)",
y = expression(R[efectivo,0]),
title = expression(R[efectivo,0]~
"según uso de condón"),
subtitle = "Línea punteada: umbral R[efectivo,0] == 1"
) +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold"))
# (4) Años hasta eliminación vs uso de condón
ggplot(sens_condom,
aes(x = condom_use * 100, y = years_to_elim)) +
geom_line(linewidth = 1.2, color = "#31688EFF") +
geom_point(aes(color = elim_dentro_hz),
size = 2) +
scale_color_manual(
values = c(`TRUE` = "#2ECC71", `FALSE` = "#E74C3C"),
name = "Eliminación\ndentro del horizonte"
) +
geom_hline(yintercept = horiz_sim,
linetype = "dashed", color = "red") +
labs(
x = "Uso de condón en actos sexuales (%)",
y = "Años hasta eliminación (T_elim)",
title = "Tiempo hasta eliminación según incremento del uso de condón",
subtitle = paste0("Línea roja: horizonte de simulación = ",
horiz_sim, " años")
) +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold"))## Warning: Removed 15 rows containing missing values or values outside the scale range
## (`geom_line()`).## Warning: Removed 15 rows containing missing values or values outside the scale range
## (`geom_point()`).
#14.3 Análisis 2D (grid parejas HR × condón) ajustado
## 14.3 Análisis 2D: combinación de parejas HR y uso de condón
grid_behav <- tidyr::expand_grid(
partners_hr = partners_seq,
condom_use = condom_seq
)
sens_grid <- grid_behav %>%
dplyr::rowwise() %>%
dplyr::mutate(
stats = list({
p <- parms_base
p$partners_year["HR"] <- partners_hr
p$condom_use <- condom_use
compute_R_eff_abs_elim(
parms = p,
init = init,
inc1_ref = inc1_ref,
alpha = alpha_behav
)
}),
R_eff_0 = stats$R_eff_0,
years_to_elim = stats$years_to_elim,
elim_dentro_hz = !is.na(years_to_elim) &&
years_to_elim <= horiz_sim
) %>%
dplyr::ungroup() %>%
dplyr::mutate(
condom_pct = condom_use * 100
)
sens_grid## # A tibble: 306 × 7
## partners_hr condom_use stats R_eff_0 years_to_elim elim_dentro_hz
## <dbl> <dbl> <list> <dbl> <dbl> <lgl>
## 1 4 0.1 <named list [3]> 0.0596 NA FALSE
## 2 4 0.15 <named list [3]> 0.0561 NA FALSE
## 3 4 0.2 <named list [3]> 0.0526 NA FALSE
## 4 4 0.25 <named list [3]> 0.0492 NA FALSE
## 5 4 0.3 <named list [3]> 0.0458 NA FALSE
## 6 4 0.35 <named list [3]> 0.0425 NA FALSE
## 7 4 0.4 <named list [3]> 0.0392 NA FALSE
## 8 4 0.45 <named list [3]> 0.0360 NA FALSE
## 9 4 0.5 <named list [3]> 0.0328 NA FALSE
## 10 4 0.55 <named list [3]> 0.0296 NA FALSE
## # ℹ 296 more rows
## # ℹ 1 more variable: condom_pct <dbl>17.3 16.4 Heatmaps 2D: R_efectivo,0 y T_elim en función de c_HR y uso de condón
# (1) Heatmap de R_efectivo,0
ggplot(sens_grid,
aes(x = partners_hr,
y = condom_pct,
fill = R_eff_0)) +
geom_raster() +
scale_fill_viridis_c(
option = "plasma",
name = expression(R[efectivo,0]),
na.value = "grey80"
) +
geom_contour(aes(z = R_eff_0),
breaks = 1,
color = "white",
linewidth = 0.6) +
labs(
x = "Número de parejas/año en HRMSM",
y = "Uso de condón (%)",
title = expression("Mapa de calor de"~R[efectivo,0]),
subtitle = "Línea blanca: frontera aproximada R[efectivo,0] == 1"
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(face = "bold")
)## Warning: `stat_contour()`: Zero contours were generated## Warning in min(x): no non-missing arguments to min; returning Inf## Warning in max(x): no non-missing arguments to max; returning -Inf
# (2) Heatmap de años hasta eliminación (T_elim)
ggplot(sens_grid,
aes(x = partners_hr,
y = condom_pct,
fill = years_to_elim)) +
geom_raster() +
scale_fill_viridis_c(
option = "magma",
name = "Años hasta eliminación",
na.value = "grey80"
) +
geom_contour(
aes(z = years_to_elim),
breaks = horiz_sim,
color = "cyan",
linewidth = 0.7
) +
labs(
x = "Número de parejas/año en HRMSM",
y = "Uso de condón (%)",
title = "Mapa de calor del tiempo hasta eliminación (T_elim)",
subtitle = paste0("Línea cian: combinaciones que logran eliminación\n",
"dentro del horizonte de simulación (", horiz_sim, " años)")
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(face = "bold")
)## Warning: Removed 243 rows containing non-finite outside the scale range
## (`stat_contour()`).## Warning: `stat_contour()`: Zero contours were generated## Warning in min(x): no non-missing arguments to min; returning Inf## Warning in max(x): no non-missing arguments to max; returning -Inf
18 17. Intervenciones biomédico-estructurales aisladas
Se consideran dos familias de intervenciones, aplicadas individualmente sobre el escenario de status quo calibrado:
TasP/Test-and-Treat aislado: aumento en la proporción anual de HSH tamizados y en las tasas de diagnóstico e inicio de TAR.
PrEP aislada: aumento en la tasa de entrada a PrEP entre susceptibles, con alta adherencia. ## 15.1 Parámetros generales del análisis
# Horizonte de simulación en años (ya fijado en parms_base)
horiz_sim <- parms_base$t_end - parms_base$t_start # p.ej. 40 años
# Nivel de severidad del umbral de eliminación (fracción de I1_ref)
alpha_biomed <- 0.10 # 10% de la incidencia media año 1 en Status quo
library(dplyr)
library(tidyr)
library(purrr)
library(ggplot2)
library(viridis)18.1 17.1 TasP aislado: prueba y tratamiento
Aquí se considera una familia de escenarios en los que se mejora solo:
La proporción de HSH que se realiza prueba anual (tested_12m_prop).
Las tasas de diagnóstico en infectados (psi4, psi7, psi10).
Las tasas de inicio de TAR (alpha5, alpha8, alpha11).
Se parametriza la intensidad de TasP mediante la proporción de tamizaje anual objetivo 𝑝test∈[0.2,0.95]ptest ∈[0.2,0.95] (igual para LR y HR). Desde el valor base (20%) hasta niveles casi universales (95%). La idea es que, para cada 𝑝testp test: Se fijan tested_12m_prop = c(LR, HR) = p_test. Se deriva la tasa de prueba anual psi1_* mediante: 𝜓1= −log (1−𝑝test) Se escala proporcionalmente las tasas de diagnóstico e inicio de TAR en función de scale_tasp=𝑝test𝑝base,scale_tasp=pbase ptest ,donde 𝑝base=0.20 (status quo).
# Secuencia de coberturas de tamizaje anual (LR = HR)
test_cov_seq <- seq(0.20, 0.95, by = 0.05) # 20% a 95%
# Proporción de tamizaje base (Status quo calibrado)
p_test_base <- parms_base$tested_12m_prop["LR"]
# Función que construye un escenario TasP dado p_test (proporción tamizada en 12 meses)
make_parms_tasp_only <- function(base_parms, p_test) {
p <- base_parms
# 1) Proporción tamizada en 12 meses
p$tested_12m_prop <- c(LR = p_test, HR = p_test)
p$psi1_LR <- derive_psi_from_prop(p_test)
p$psi1_HR <- derive_psi_from_prop(p_test)
# 2) Factor multiplicativo respecto al status quo
scale_tasp <- p_test / p_test_base
# Escalar tasas de diagnóstico
p$psi4_LR <- base_parms$psi4_LR * scale_tasp
p$psi4_HR <- base_parms$psi4_HR * scale_tasp
p$psi7_LR <- base_parms$psi7_LR * scale_tasp
p$psi7_HR <- base_parms$psi7_HR * scale_tasp
p$psi10_LR <- base_parms$psi10_LR * scale_tasp
p$psi10_HR <- base_parms$psi10_HR * scale_tasp
# Escalar tasas de inicio de TAR
p$alpha5 <- base_parms$alpha5 * scale_tasp
p$alpha8 <- base_parms$alpha8 * scale_tasp
p$alpha11 <- base_parms$alpha11 * scale_tasp
# Meta descriptiva de cobertura de TAR (no entra directamente en las EDO)
p$art_cov_target <- min(0.98, 0.80 + 0.18 * (p_test - p_test_base) / (0.95 - p_test_base))
p
}
# Barrido de escenarios TasP
sens_tasp <- map_dfr(test_cov_seq, function(p_test) {
p_tasp <- make_parms_tasp_only(parms_base, p_test)
stats <- compute_R_eff_abs_elim(
parms = p_tasp,
init = init,
inc1_ref = inc1_ref,
alpha = alpha_biomed
)
tibble(
tipo = "TasP",
test_cov = p_test,
test_cov_pct = p_test * 100,
R_eff_0 = stats$R_eff_0,
inc1 = stats$inc1,
years_to_elim = stats$years_to_elim,
elim_dentro_hz = !is.na(stats$years_to_elim) &&
stats$years_to_elim <= horiz_sim
)
})
sens_tasp## # A tibble: 16 × 7
## tipo test_cov test_cov_pct R_eff_0 inc1 years_to_elim elim_dentro_hz
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <lgl>
## 1 TasP 0.2 20 0.0822 0.00250 NA FALSE
## 2 TasP 0.25 25 0.0806 0.00243 NA FALSE
## 3 TasP 0.3 30 0.0791 0.00237 NA FALSE
## 4 TasP 0.35 35 0.0777 0.00231 NA FALSE
## 5 TasP 0.4 40 0.0765 0.00225 NA FALSE
## 6 TasP 0.45 45 0.0753 0.00220 NA FALSE
## 7 TasP 0.5 50 0.0741 0.00215 NA FALSE
## 8 TasP 0.55 55 0.0730 0.00211 NA FALSE
## 9 TasP 0.6 60 0.0720 0.00207 NA FALSE
## 10 TasP 0.65 65 0.0710 0.00203 NA FALSE
## 11 TasP 0.7 70 0.0701 0.00200 NA FALSE
## 12 TasP 0.75 75 0.0692 0.00196 NA FALSE
## 13 TasP 0.8 80 0.0683 0.00193 NA FALSE
## 14 TasP 0.85 85 0.0674 0.00190 NA FALSE
## 15 TasP 0.9 90 0.0666 0.00187 NA FALSE
## 16 TasP 0.95 95 0.0658 0.00184 NA FALSE# Cobertura mínima que logra eliminación dentro del horizonte, si existe
tasp_min_cov <- sens_tasp %>%
filter(elim_dentro_hz) %>%
summarise(min_cov = min(test_cov_pct)) %>%
pull(min_cov)## Warning: There was 1 warning in `summarise()`.
## ℹ In argument: `min_cov = min(test_cov_pct)`.
## Caused by warning in `min()`:
## ! no non-missing arguments to min; returning Inftasp_min_cov # puede ser NA si ningún nivel logra eliminar## [1] Inf18.2 17.2 Gráficos para TasP aislado
## R_efectivo,0 vs cobertura de testeo anual (TasP)
ggplot(sens_tasp,
aes(x = test_cov_pct, y = R_eff_0)) +
geom_line(linewidth = 1.2, color = "#440154FF") +
geom_point(aes(color = elim_dentro_hz), size = 2) +
scale_color_manual(
values = c(`TRUE` = "#2ECC71", `FALSE` = "#E74C3C"),
name = "Eliminación\ndentro del horizonte"
) +
geom_hline(yintercept = 1,
linetype = "dashed",
color = "grey50") +
labs(
x = "Cobertura de test VIH en 12 meses (%)",
y = expression(R[efectivo,0]),
title = expression(R[efectivo,0]~
"según cobertura de test (TasP aislado)"),
subtitle = "Línea punteada: umbral R[efectivo,0] == 1"
) +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold"))
## Tiempo hasta eliminación vs cobertura de testeo (TasP)
ggplot(sens_tasp,
aes(x = test_cov_pct, y = years_to_elim)) +
geom_line(linewidth = 1.2, color = "#31688EFF") +
geom_point(aes(color = elim_dentro_hz), size = 2) +
scale_color_manual(
values = c(`TRUE` = "#2ECC71", `FALSE` = "#E74C3C"),
name = "Eliminación\ndentro del horizonte"
) +
geom_hline(yintercept = horiz_sim,
linetype = "dashed",
color = "red") +
labs(
x = "Cobertura de test VIH en 12 meses (%)",
y = "Años hasta eliminación (T_elim)",
title = "Tiempo hasta eliminación según cobertura TasP aislada",
subtitle = paste0("Línea roja: horizonte de simulación = ",
horiz_sim, " años")
) +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold"))## Warning: Removed 16 rows containing missing values or values outside the scale range
## (`geom_line()`).## Warning: Removed 16 rows containing missing values or values outside the scale range
## (`geom_point()`).
18.3 17.3 PrEP aislada En esta familia de escenarios se aumenta solo la intensidad de PrEP entre susceptibles, manteniendo:
Tasas de testeo y TAR en valores de status quo calibrado.
Conducta sexual (número de parejas y uso de condón) sin cambios adicionales.
Se usa como variable de control una cobertura objetivo de PrEP 𝑐PrEP (fracción de susceptibles en PrEP), en el rango 0–40 %. En la práctica, la cobertura dinámica está determinada por la tasa de entrada a PrEP 𝜙2 y la tasa de abandono 𝜙3 . Aquí se adopta una relación simple: 𝜙3=0 (abandono despreciable en el escenario ideal). 𝜙2≈1,2×𝑐PrEP(hazard anual aproximado que tiende a producir esa cobertura).
Además, se asume una adherencia alta pero realista, por ejemplo 85 %.
## 15.3 PrEP aislada
# Secuencia de coberturas objetivo de PrEP (0 a 40%)
prep_cov_seq <- seq(0.00, 0.40, by = 0.02)
# Función que construye un escenario de PrEP aislada para una cobertura objetivo dada
make_parms_prep_only <- function(base_parms, prep_cov_target) {
p <- base_parms
# No cambiar testeo ni TAR respecto al status quo
# (p$tested_12m_prop, psi4/7/10, alpha* permanecen como en parms_base)
# Cobertura "objetivo" (descriptiva)
p$prep_cov_target <- prep_cov_target
# Tasa de entrada a PrEP: aproximación lineal basada en escenarios previos
# (para 20% se usaba phi2 ≈ 0.25; para 30% φ2 ≈ 0.35)
p$phi2 <- 1.2 * prep_cov_target # hazard anual aproximado
p$phi3 <- 0.0 # abandono despreciable
# Adherencia alta pero < 100%
p$prep_adherence <- 0.85
p
}
# Barrido de escenarios PrEP aislada
sens_prep <- map_dfr(prep_cov_seq, function(cov_prep) {
p_prep <- make_parms_prep_only(parms_base, cov_prep)
stats <- compute_R_eff_abs_elim(
parms = p_prep,
init = init,
inc1_ref = inc1_ref,
alpha = alpha_biomed
)
tibble(
tipo = "PrEP",
prep_cov = cov_prep,
prep_cov_pct = cov_prep * 100,
R_eff_0 = stats$R_eff_0,
inc1 = stats$inc1,
years_to_elim = stats$years_to_elim,
elim_dentro_hz = !is.na(stats$years_to_elim) &&
stats$years_to_elim <= horiz_sim
)
})
sens_prep## # A tibble: 21 × 7
## tipo prep_cov prep_cov_pct R_eff_0 inc1 years_to_elim elim_dentro_hz
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <lgl>
## 1 PrEP 0 0 0.0823 0.00250 NA FALSE
## 2 PrEP 0.02 2 0.0822 0.00250 NA FALSE
## 3 PrEP 0.04 4 0.0820 0.00250 NA FALSE
## 4 PrEP 0.06 6 0.0818 0.00249 NA FALSE
## 5 PrEP 0.08 8 0.0817 0.00249 NA FALSE
## 6 PrEP 0.1 10 0.0815 0.00248 NA FALSE
## 7 PrEP 0.12 12 0.0814 0.00248 NA FALSE
## 8 PrEP 0.14 14 0.0812 0.00247 NA FALSE
## 9 PrEP 0.16 16 0.0811 0.00247 NA FALSE
## 10 PrEP 0.18 18 0.0809 0.00247 NA FALSE
## # ℹ 11 more rows# Cobertura mínima de PrEP que logra eliminación dentro del horizonte (si existe)
prep_min_cov <- sens_prep %>%
filter(elim_dentro_hz) %>%
summarise(min_cov = min(prep_cov_pct)) %>%
pull(min_cov)## Warning: There was 1 warning in `summarise()`.
## ℹ In argument: `min_cov = min(prep_cov_pct)`.
## Caused by warning in `min()`:
## ! no non-missing arguments to min; returning Infprep_min_cov # puede ser NA si ningún nivel lo logra## [1] Inf18.4 17.4. Gráficos para PrEP aislada
## R_efectivo,0 vs cobertura de PrEP
ggplot(sens_prep,
aes(x = prep_cov_pct, y = R_eff_0)) +
geom_line(linewidth = 1.2, color = "#440154FF") +
geom_point(aes(color = elim_dentro_hz), size = 2) +
scale_color_manual(
values = c(`TRUE` = "#2ECC71", `FALSE` = "#E74C3C"),
name = "Eliminación\ndentro del horizonte"
) +
geom_hline(yintercept = 1,
linetype = "dashed",
color = "grey50") +
labs(
x = "Cobertura objetivo de PrEP en susceptibles (%)",
y = expression(R[efectivo,0]),
title = expression(R[efectivo,0]~
"según cobertura de PrEP (aislada)"),
subtitle = "Línea punteada: umbral R[efectivo,0] == 1"
) +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold"))
## Tiempo hasta eliminación vs cobertura de PrEP
ggplot(sens_prep,
aes(x = prep_cov_pct, y = years_to_elim)) +
geom_line(linewidth = 1.2, color = "#31688EFF") +
geom_point(aes(color = elim_dentro_hz), size = 2) +
scale_color_manual(
values = c(`TRUE` = "#2ECC71", `FALSE` = "#E74C3C"),
name = "Eliminación\ndentro del horizonte"
) +
geom_hline(yintercept = horiz_sim,
linetype = "dashed",
color = "red") +
labs(
x = "Cobertura objetivo de PrEP en susceptibles (%)",
y = "Años hasta eliminación (T_elim)",
title = "Tiempo hasta eliminación según cobertura de PrEP aislada",
subtitle = paste0("Línea roja: horizonte de simulación = ",
horiz_sim, " años")
) +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold"))## Warning: Removed 21 rows containing missing values or values outside the scale range
## (`geom_line()`).## Warning: Removed 21 rows containing missing values or values outside the scale range
## (`geom_point()`).
19 18. Combinaciones de intervención y rutas factibles hacia la eliminación (Objetivo 6)
En este modelo:
- El horizonte de simulación se definió a partir de un año base
(
base_year), y un tiempo finalt_end(en años) enparms_base: \[ \text{horizonte simulado} = [\text{año inicial}, \text{año final}] = [\min(\text{year}), \max(\text{year})], \] dondeyearproviene deyearly_all2. - Los escenarios evaluados son:
- Status quo
- TasP 80%
- PrEP 20%
- Combinado moderado
- Combinado fuerte
Para responder al objetivo 6, combinamos tres tipos de información:
- Carga epidémica acumulada por escenario (infecciones totales, infecciones evitadas y % de reducción frente a Status quo).
- Indicadores de eliminación, obtenidos a partir de:
- el número reproductivo efectivo promedio en el año inicial \(R_{\text{efectivo},0}\),
- el tiempo hasta eliminación \(T_{\text{elim}}\) para un umbral relativo: \[ \tau = \alpha\,I_1^{\text{ref}}, \] donde \(\alpha\) es un nivel de política (por ejemplo, 10 % de la incidencia de referencia).
- Clasificación de escenarios según si logran alcanzar el umbral de eliminación dentro del horizonte simulado y cuántas infecciones logran evitar.
En lo que sigue se usa la tabla de resultados anuales
yearly_all2, la tabla de escenarios
tabla_escenarios (sección 10.4) y la tabla de \(R_{\text{efectivo},0}\) y \(T_{\text{elim}}\) con distintos niveles de
\(\alpha\)
(tabla_Reff_abs_alpha, sección 12).
19.1 18.1. Horizonte temporal efectivo de la simulación
#Rango de años realmente simulados
year_start <- min(yearly_all2$year, na.rm = TRUE)
year_end <- max(yearly_all2$year, na.rm = TRUE)
cat("Horizonte temporal de la simulación: de", year_start, "a", year_end, "\n")## Horizonte temporal de la simulación: de 2025 a 206519.2 ##18.2. Resumen conjunto: infecciones evitadas + R_efectivo,0 + T_elim Elegimos un nivel de 𝛼 α (por ejemplo, 10 %) como criterio principal de eliminación:
Umbral de eliminación:𝜏=0,10×𝐼1ref donde:𝐼1ref es la incidencia media del primer año en Status quo.
𝑇elim es el primer tiempo en que la incidencia instantánea cae por debajo de 𝜏
#Nivel de política para el umbral de eliminación (en % de la incidencia de referencia)
alpha_policy_pct <- 10
alpha_policy <- alpha_policy_pct / 100
#Tabla de R_efectivo,0 y T_elim para ese alpha puntual
tabla_elim_policy <- tabla_Reff_abs_alpha %>%
dplyr::filter(alpha_pct == alpha_policy_pct) %>%
dplyr::select(
Escenario,
R_eff_0 = R_eff_0,
T_elim = years_to_elim
)
tabla_elim_policy## # A tibble: 5 × 3
## Escenario R_eff_0 T_elim
## <chr> <dbl> <dbl>
## 1 Status quo 0.0822 NA
## 2 TasP 80% 0.0791 NA
## 3 PrEP 20% 0.0807 NA
## 4 Combinado moderado 0.0633 NA
## 5 Combinado fuerte 0.0390 4.1Ahora unimos estos indicadores con la tabla de escenarios tabla_escenarios, que ya resume:
Prevalencia final,
Incidencia final,
Infecciones totales,
Infecciones evitadas (absolutas y en %).
#Unir tabla de escenarios (carga acumulada) con R_eff_0 y T_elim
tabla_ruta_eliminacion <- tabla_escenarios %>%
dplyr::rename(Escenario = Escenario) %>%
dplyr::left_join(
tabla_elim_policy,
by = "Escenario"
) %>%
dplyr::mutate(
logro_eliminacion = !is.na(T_elim) & T_elim <= (parms_base$t_end + 1e-6)
)
tabla_ruta_eliminacion## # A tibble: 5 × 10
## Escenario Año `Prevalencia (%)` `Incidencia (/1000 p-año)`
## <chr> <dbl> <dbl> <dbl>
## 1 Status quo 2065 10.4 9.51
## 2 TasP 80% 2065 9.66 8.88
## 3 PrEP 20% 2065 3.28 2.34
## 4 Combinado moderado 2065 1.01 0.77
## 5 Combinado fuerte 2065 0.22 0.704
## # ℹ 6 more variables: `Inf. totales (×10.000)` <dbl>,
## # `Inf. evitadas (×10.000)` <dbl>, `Inf. evitadas (%)` <dbl>, R_eff_0 <dbl>,
## # T_elim <dbl>, logro_eliminacion <lgl>19.3 18.3. Escenarios que logran eliminación y ranking por infecciones evitadas
#Escenarios que sí alcanzan el umbral de eliminación dentro del horizonte
escenarios_eliminacion <- tabla_ruta_eliminacion %>%
dplyr::filter(logro_eliminacion) %>%
dplyr::arrange(T_elim)
escenarios_eliminacion## # A tibble: 1 × 10
## Escenario Año `Prevalencia (%)` `Incidencia (/1000 p-año)`
## <chr> <dbl> <dbl> <dbl>
## 1 Combinado fuerte 2065 0.22 0.704
## # ℹ 6 more variables: `Inf. totales (×10.000)` <dbl>,
## # `Inf. evitadas (×10.000)` <dbl>, `Inf. evitadas (%)` <dbl>, R_eff_0 <dbl>,
## # T_elim <dbl>, logro_eliminacion <lgl>#Escenarios ordenados por infecciones evitadas (%), excluyendo Status quo
escenarios_mas_efectivos <- tabla_ruta_eliminacion %>%
dplyr::filter(
Escenario != "Status quo",
!is.na(`Inf. evitadas (%)`)
) %>%
dplyr::arrange(dplyr::desc(`Inf. evitadas (%)`))
escenarios_mas_efectivos## # A tibble: 4 × 10
## Escenario Año `Prevalencia (%)` `Incidencia (/1000 p-año)`
## <chr> <dbl> <dbl> <dbl>
## 1 Combinado fuerte 2065 0.22 0.704
## 2 Combinado moderado 2065 1.01 0.77
## 3 PrEP 20% 2065 3.28 2.34
## 4 TasP 80% 2065 9.66 8.88
## # ℹ 6 more variables: `Inf. totales (×10.000)` <dbl>,
## # `Inf. evitadas (×10.000)` <dbl>, `Inf. evitadas (%)` <dbl>, R_eff_0 <dbl>,
## # T_elim <dbl>, logro_eliminacion <lgl>20 Apéndice. Explicación didáctica de las ecuaciones del modelo
Esta sección resume, con un lenguaje más pedagógico, las ecuaciones que se implementan en el código del modelo. La idea es que puedas usarla casi tal cual en una presentación de maestría o ante un comité de ética.
20.1 A.1 Compartimentos y población total
Para cada grupo de riesgo \(g \in \{\text{LR}, \text{HR}\}\) se consideran 12 compartimentos:
- \(S_1^g\): susceptibles sin prueba reciente.
- \(S_2^g\): susceptibles con prueba VIH negativa reciente.
- \(S_3^g\): susceptibles en profilaxis preexposición (PrEP).
- \(A_1^g, A_2^g, A_3^g\): infección aguda (no diagnosticada, diagnosticada sin TAR, en TAR).
- \(C_1^g, C_2^g, C_3^g\): infección crónica (no diagnosticada, diagnosticada sin TAR, en TAR).
- \(D_1^g, D_2^g, D_3^g\): SIDA (no diagnosticada, diagnosticada sin TAR, en TAR).
El tamaño poblacional del grupo \(g\) en el tiempo \(t\) se define como:
\[ N_g(t) = \sum_{k=1}^3 S_k^g(t) + \sum_{k=1}^3 A_k^g(t) + \sum_{k=1}^3 C_k^g(t) + \sum_{k=1}^3 D_k^g(t). \]
La población total del modelo es simplemente:
\[ N(t) = N_{\text{LR}}(t) + N_{\text{HR}}(t). \]
20.2 A.2 Mezcla sexual y fuerza de infección
El modelo supone que cada grupo de riesgo tiene un número promedio de parejas sexuales por año \(c_g\) y que las parejas se forman según una matriz de mezcla \(H\), que describe la fracción de parejas que se dan entre y dentro de grupos:
\[ H = \begin{pmatrix} H_{\text{LR,LR}} & H_{\text{LR,HR}} \\ H_{\text{HR,LR}} & H_{\text{HR,HR}} \end{pmatrix}, \]
donde, por ejemplo, \(H_{\text{LR,HR}}\) es la proporción de parejas de personas LR con parejas HR. La suma por fila es 1:
\[ H_{\text{LR,LR}} + H_{\text{LR,HR}} = 1,\quad H_{\text{HR,LR}} + H_{\text{HR,HR}} = 1. \]
Para cada grupo \(g\) se define una infectividad promedio por pareja \(\beta_g(t)\), que promedia las probabilidades de transmisión por pareja en cada estadio (agudo, crónico, SIDA), ponderadas por el número de personas en cada compartimento y por la reducción de infectividad asociada a la TAR y al uso de condón.
En forma simplificada, si \(\sigma_s\) es la probabilidad de transmisión por pareja en el estadio \(s\) y \(\omega_s\) es un factor de reducción (por ejemplo, por TAR efectiva), la infectividad promedio en el grupo \(g\) puede escribirse como:
\[ \beta_g(t) = \frac{1}{N_g(t)} \sum_{s} \sigma_s\,\omega_s\,I_s^g(t), \]
donde \(I_s^g(t)\) denota el número de personas infecciosas en el estadio \(s\) dentro del grupo \(g\).
La fuerza de infección (hazard de infección) para los susceptibles de cada grupo se define como:
\[ \lambda_g(t) = c_g\Big[ H_{\text{g,LR}}\,\beta_{\text{LR}}(t) + H_{\text{g,HR}}\,\beta_{\text{HR}}(t) \Big], \quad g \in \{\text{LR}, \text{HR}\}. \]
En el código, esta cantidad se calcula en la función
foi_susceptibles(), que devuelve \(\lambda_{\text{LR}}(t)\) y \(\lambda_{\text{HR}}(t)\), y luego se
reparte entre \(S_1, S_2\) y \(S_3\) según los supuestos del escenario
(por ejemplo, menor riesgo entre quienes están en PrEP).
20.3 A.3 Parámetros de progresión y mortalidad
El modelo asume tiempos medios de progresión entre estadios (agudo, crónico, SIDA) expresados en años:
- \(\text{dur\_acute\_y}\): duración media de la fase aguda.
- \(\text{dur\_chronic\_y}\): duración media de la fase crónica.
- \(\text{dur\_aids\_y}\): duración media desde SIDA hasta la muerte.
A partir de ellos se obtienen las tasas de progresión: \[ \theta_{\text{acute}} = \frac{1}{\text{dur\_acute\_y}},\quad \theta_{\text{chronic}} = \frac{1}{\text{dur\_chronic\_y}},\quad \theta_{\text{aids}} = \frac{1}{\text{dur\_aids\_y}}. \]
Estas tasas se aplican de forma similar a los compartimentos \(A\), \(C\) y \(D\), independientemente de si la persona está diagnosticada o en TAR. La mortalidad general se representa por \(\mu\) y se añade un exceso de mortalidad específico por estadio \(b_s\). Por ejemplo, en SIDA sin TAR se usa un valor alto (por ejemplo, \(b_{10} = b_{11} \approx 1{,}68\) por año), lo que implica una supervivencia media menor a un año, coherente con la historia natural del VIH sin tratamiento.
Con estos parámetros, la tasa de salida total de un compartimento infeccioso típico es la suma de:
\[ \text{tasa de salida} = \text{progresión} + \text{mortalidad general} + \text{exceso de mortalidad}. \]
No se identifican valores abiertamente inverosímiles en las tasas usadas en el código; sin embargo, pueden refinarse con nueva evidencia epidemiológica o revisión de la literatura especializada.
20.4 A.4 Ecuaciones diferenciales para los susceptibles (ejemplo)
La ecuación genérica para los susceptibles de tipo 1 en el grupo de bajo riesgo LR puede escribirse como:
\[ \frac{d S_1^{\text{LR}}(t)}{dt} = \rho_{\text{LR}} - \lambda_1^{\text{LR}}(t)\, S_1^{\text{LR}}(t) - \psi_1^{\text{LR}}\, S_1^{\text{LR}}(t) - \mu\, S_1^{\text{LR}}(t), \]
donde:
- \(\rho_{\text{LR}}\) es la tasa de entrada demográfica al grupo LR.
- \(\lambda_1^{\text{LR}}(t)\) es la fuerza de infección efectiva sobre \(S_1^{\text{LR}}\), derivada de \(\lambda_{\text{LR}}(t)\).
- \(\psi_1^{\text{LR}}\) es la tasa a la que los susceptibles del tipo 1 se realizan prueba (y, en caso de ser negativos, pasan a \(S_2\)).
- \(\mu\) es la mortalidad general.
Ecuaciones análogas gobiernan \(S_2^g\) y \(S_3^g\), incorporando además las tasas de inicio y abandono de PrEP (\(\phi_2, \phi_3\)) y el retorno eventual desde PrEP a susceptibilidad estándar si el programa se interrumpe.
20.5 A.5 Ecuaciones para infección aguda, crónica y SIDA (esquema general)
De forma esquemática, para los compartimentos agudos en LR se tiene:
\[ \frac{d A_1^{\text{LR}}}{dt} = \lambda_{\text{LR}}(t)\big(S_1^{\text{LR}} + S_2^{\text{LR}} + S_3^{\text{LR}}\big) - (\psi_4^{\text{LR}} + \theta_4 + \mu + b_4)\, A_1^{\text{LR}}, \]
\[ \frac{d A_2^{\text{LR}}}{dt} = \psi_4^{\text{LR}} A_1^{\text{LR}} - (\alpha_5 + \theta_5 + \mu + b_5)\, A_2^{\text{LR}}, \]
\[ \frac{d A_3^{\text{LR}}}{dt} = \alpha_5 A_2^{\text{LR}} - (\theta_6 + \mu + b_6)\, A_3^{\text{LR}}. \]
- \(\psi_4^{\text{LR}}\) representa la tasa de diagnóstico en fase aguda.
- \(\alpha_5\) es la tasa de inicio de TAR en agudo diagnosticado.
- \(\theta_4, \theta_5, \theta_6\) son las tasas de progresión en los tres compartimentos agudos.
- \(b_4, b_5, b_6\) son los excesos de mortalidad específicos de cada compartimento.
El mismo patrón se replica para la fase crónica (\(C_1, C_2, C_3\)) y la fase de SIDA (\(D_1, D_2, D_3\)), con sus respectivas tasas de diagnóstico y de inicio de TAR (\(\psi_7, \psi_{10}, \alpha_8, \alpha_{11}\)). El código organiza estas ecuaciones en bloques dentro de la función que define el sistema de EDO.
20.6 A.6 Indicadores epidemiológicos
A partir de las trayectorias simuladas se calculan indicadores estándar:
Prevalencia en el grupo \(g\): \[ \text{prev}_g(t) = \frac{I_g(t)}{N_g(t)}, \] donde \(I_g(t)\) es la suma de todos los compartimentos infecciosos en \(g\).
Incidencia instantánea (tasa de casos nuevos por persona-año): \[ \text{inc}(t) = \frac{\text{nuevas infecciones en } [t,t+\Delta t)}{N(t)\,\Delta t}. \]
En la práctica se integra numéricamente la incidencia a lo largo del tiempo para obtener:
- Incidencia media en el primer año.
- Número acumulado de infecciones nuevas en el horizonte de simulación.
- Infecciones evitadas al comparar un escenario de intervención con el escenario de referencia (status quo).
20.7 A.7 Criterio operativo de eliminación y número reproductivo efectivo aproximado
Se define un umbral de eliminación basado en una fracción \(\alpha\) de la incidencia inicial de referencia \(I_{\text{ref}}(0)\):
\[ \text{Incidencia objetivo} = \alpha \times I_{\text{ref}}(0). \]
El tiempo hasta la eliminación es el primer instante \(t\) en que la incidencia del escenario considerado cae por debajo de ese umbral.
Además, se usa un proxy del número reproductivo efectivo en el año 0, \(R_{\text{efectivo},0}\), calculado como la razón entre la incidencia media del primer año en el escenario de intervención y la incidencia media del primer año en el escenario de referencia. Valores menores que 1 sugieren que la epidemia tendería a decrecer bajo ese escenario:
\[ R_{\text{efectivo},0} \approx \frac{\text{Incidencia media año 1 (escenario)}}{\text{Incidencia media año 1 (referencia)}}. \]
Este enfoque es coherente con la teoría de modelos compartimentales, pero mantiene una interpretación intuitiva y fácilmente comunicable en una sustentación de maestría o ante un comité de ética.