library(shiny) library(ggplot2) library(dplyr) library(scales) # ============================================================ # UI # ============================================================ ui <- fluidPage( tags$head(tags$title("Modelos de Mortalidad — JML")), titlePanel( tagList( h2("Modelos de Mortalidad/Supervivencia ", style = "margin-bottom: 4px; color: #1a3a6b; font-size: 22px;"), h3("JuanMLejarza"), h5("Dormoy 1 · Dormoy 2 · Sang · Gompertz · Makeham · predefinidos · personalizada", style = "color: gray; margin-top: 0; font-size: 14px;") ) ), sidebarLayout( sidebarPanel( radioButtons("modo", "Función de cohorte l(x):", choices = c( "Dormoy 1 (S)" = "dormoy1", "Dormoy 2 (S1, S2)" = "dormoy2", "Sang (B, ω)" = "sang", "Gompertz (g, c)" = "gompertz", "Makeham (s, g, c)" = "makeham", "Modelo predefinido" = "predef", "Expresión personalizada" = "custom" ), selected = "dormoy2" ), hr(), # Dormoy 1 conditionalPanel( condition = "input.modo == 'dormoy1'", h4("Parámetros Dormoy 1"), sliderInput("s_d1", "S (supervivencia base):", min = 0.98, max = 0.999, value = 0.997, step = 0.001), helpText("l(x) = l₀ · S^x. Tanto constante μ = -ln(S)."), numericInput("radix_d1", "Población inicial (l₀):", value = 100000, min = 1000, max = 1000000) ), # Dormoy 2 conditionalPanel( condition = "input.modo == 'dormoy2'", h4("Parámetros Dormoy 2"), sliderInput("s1_d2", "S1 (supervivencia base):", min = 0.98, max = 0.999, value = 0.997, step = 0.001), sliderInput("s2_d2", "S2 (aceleración edad):", min = 0.999, max = 0.99999, value = 0.9999, step = 0.00001), helpText("l(x) = l₀ · S1^x · S2^(x²). Tanto lineal μ = -2x·ln(S2) - ln(S1)."), numericInput("radix_d2", "Población inicial (l₀):", value = 100000, min = 1000, max = 1000000) ), # Sang conditionalPanel( condition = "input.modo == 'sang'", h4("Parámetros Ley de Sang"), sliderInput("b_sang", "B (base de la exponencial):", min = 0.80, max = 0.99, value = 0.97, step = 0.01), numericInput("omega_sang", "ω (edad límite):", value = 110, min = 2, max = 400), helpText("l(x) = l₀ · (B^x - B^ω)/(1 - B^ω)."), numericInput("radix_sang", "Población inicial (l₀):", value = 100000, min = 1000, max = 1000000) ), # Gompertz — parametrizado con g y c conditionalPanel( condition = "input.modo == 'gompertz'", h4("Parámetros Gompertz"), numericInput("g_gompertz", "g (factor de supervivencia):", value = 0.97, min = 0.01, max = 0.99, step = 0.01), numericInput("c_gompertz", "c (base exponencial, >1):", value = 1.03, min = 1.001, max = 2, step = 0.01), helpText("l(x) = l₀ · g^(c^x - 1). Tanto μ(x) = -ln(g)·ln(c)·c^x."), numericInput("radix_gompertz", "Población inicial (l₀):", value = 1000, min = 1000, max = 1000000) ), # Makeham — parametrizado con s, g y c conditionalPanel( condition = "input.modo == 'makeham'", h4("Parámetros Makeham"), numericInput("s_makeham", "s (factor de supervivencia base):", value = 0.99, min = 0.01, max = 0.99, step = 0.01), numericInput("g_makeham", "g (factor de envejecimiento):", value = 0.97, min = 0.01, max = 0.99, step = 0.01), numericInput("c_makeham", "c (base exponencial, >1):", value = 1.03, min = 1.001, max = 2, step = 0.01), helpText("l(x) = l₀ · s^x · g^(c^x - 1). Tanto μ(x) = -ln(s) - ln(g)·ln(c)·c^x."), numericInput("radix_makeham", "Población inicial (l₀):", value = 1000, min = 1000, max = 1000000) ), # Predefinido conditionalPanel( condition = "input.modo == 'predef'", h4("Modelos predefinidos"), selectInput("modelo_predef", "Elige el modelo:", choices = c( "Cuadrática: 10000 - x²" = "cuadratica", "De Moivre: (100 - x)/190" = "moivre190", "De Moivre: (100 - x)/150" = "moivre150", "Lineal normalizada: (100-x)/100" = "lineal100", "Exponencial: 1000·e^(-0.01x)" = "exponencial" ) ), numericInput("radix_predef", "Población inicial (l₀):", value = 100000, min = 1000, max = 1000000) ), # Personalizado conditionalPanel( condition = "input.modo == 'custom'", h4("Expresión personalizada"), textInput("lx_expr", "Escribe l(x):", value = "10000 - x^2", placeholder = "Ej: 1000 * exp(-0.01 * x)"), helpText("Usa x como variable, punto decimal y exp() para exponenciales."), numericInput("radix_custom", "Población inicial (l₀):", value = 100000, min = 1000, max = 1000000) ), hr(), sliderInput("edad", "Edad actual (x):", min = 0, max = 100, value = 40), sliderInput("omega", "Edad máxima (ω):", min = 50, max = 400, value = 100), actionButton("calcular", "Calcular", class = "btn-primary") ), mainPanel( tabsetPanel( tabPanel("Gráficos", plotOutput("grafico_supervivencia"), plotOutput("grafico_esperanza"), plotOutput("grafico_mortalidad") ), tabPanel("Resultados", h4("Indicadores Clave"), tableOutput("indicadores"), h4("Tabla de Mortalidad (resumen)"), tableOutput("tabla_resumen") ), tabPanel("eₓ vía Tₓ", withMathJax(), h4("Definición de Tₓ"), helpText("$$T_x = \\int_x^{\\omega} l(s)\\,ds = \\int_0^{\\omega-x} l_{x+t}\\,dt$$"), h4("Esperanza de vida"), helpText("$$e_x = \\frac{T_x}{l_x}$$"), hr(), h4("Valores a la edad seleccionada"), tableOutput("tabla_Tx"), hr(), h4("Gráfico de Tₓ frente a la edad"), plotOutput("grafico_Tx") ), tabPanel("Fórmulas", withMathJax(), h4("Función de cohorte actual"), uiOutput("formula_lx"), hr(), h4("Esperanza de vida:"), helpText("$$e_x = \\frac{T_x}{l_x} = \\frac{\\int_x^{\\omega} l(s)\\,ds}{l_x}$$"), h4("Tasa instantánea de mortalidad:"), helpText("$$\\mu_x = -\\frac{d}{dx}\\ln l_x$$") ), tabPanel("Fórmulas eₓ vía μₓ", withMathJax(), HTML('

Paso 1: Supervivencia

$${}_t p_x = \\frac{l_{x+t}}{l_x}$$

Paso 2: Fuerza de mortalidad

$$\\mu_{x+t} = -\\frac{d}{ds}\\ln l_s\\Big|_{s=x+t}$$

Paso 3: Densidad

$$f_T(t) = {}_t p_x \\cdot \\mu_{x+t}$$

Paso 4: Esperanza

$$e_x = \\int_0^{\\omega-x} t \\cdot {}_t p_x \\cdot \\mu_{x+t}\\, dt$$

Paso 5: Fórmula cerrada (si existe)

'), uiOutput("paso5_formula") ) ) ) ) ) # ============================================================ # SERVER # ============================================================ server <- function(input, output, session) { # ---------- Función l(x) reactiva ---------- lx_func <- reactive({ if (input$modo == "dormoy1") { S <- input$s_d1; l0 <- input$radix_d1 function(x) l0 * (S^x) } else if (input$modo == "dormoy2") { S1 <- input$s1_d2; S2 <- input$s2_d2; l0 <- input$radix_d2 function(x) l0 * (S1^x) * (S2^(x^2)) } else if (input$modo == "sang") { B <- input$b_sang; w <- input$omega_sang; l0 <- input$radix_sang function(x) l0 * ((B^x - B^w) / (1 - B^w)) } else if (input$modo == "gompertz") { g <- input$g_gompertz; c <- input$c_gompertz; l0 <- input$radix_gompertz function(x) l0 * (g^(c^x - 1)) } else if (input$modo == "makeham") { s <- input$s_makeham; g <- input$g_makeham; c <- input$c_makeham l0 <- input$radix_makeham function(x) l0 * (s^x) * (g^(c^x - 1)) } else if (input$modo == "predef") { l0 <- input$radix_predef switch(input$modelo_predef, "cuadratica" = function(x) (10000 - x^2) * (l0 / 10000), "moivre190" = function(x) ((100 - x) / 190) * l0, "moivre150" = function(x) ((100 - x) / 150) * l0, "lineal100" = function(x) ((100 - x) / 100) * l0, "exponencial" = function(x) (1000 * exp(-0.01 * x)) * (l0 / 1000) ) } else { expr <- tryCatch(parse(text = input$lx_expr), error = function(e) NULL) validate(need(!is.null(expr), "❌ Expresión no válida. Revisa la sintaxis.")) l0 <- input$radix_custom base_fun <- function(x) eval(expr, envir = list(x = x)) function(x) base_fun(x) * (l0 / base_fun(0)) } }) # ---------- Validación ---------- observe({ f <- lx_func(); omega <- input$omega test <- tryCatch(f(seq(0, omega, length.out = 50)), error = function(e) rep(NA_real_, 50)) if (any(is.na(test))) { showNotification("⚠️ l(x) produce NA. Revisa la expresión.", type = "error", duration = 5) } else if (any(test < 0)) { showNotification("⚠️ l(x) toma valores negativos.", type = "warning", duration = 5) } else if (abs(f(omega)) > 0.05 * f(0)) { showNotification(sprintf("⚠️ l(ω) = %.4g (no es ≈ 0). Ajusta ω o la expresión.", f(omega)), type = "warning", duration = 5) } }) # ---------- Tabla de mortalidad ---------- calcular_tabla <- reactive({ req(input$calcular) edad <- 0:input$omega n <- length(edad) f <- lx_func() l_x <- tryCatch(f(edad), error = function(e) rep(NA_real_, n)) validate(need(all(!is.na(l_x)), "❌ l(x) no se puede evaluar en todo el rango.")) validate(need(all(l_x >= 0), "❌ l(x) es negativa en algún punto.")) d_x <- c(-diff(l_x), l_x[n]) L_x <- numeric(n) for (i in 1:(n - 1)) L_x[i] <- (l_x[i] + l_x[i + 1]) / 2 L_x[n] <- l_x[n] * 0.5 T_x <- rev(cumsum(rev(L_x))) e_x <- T_x / l_x e_x[is.nan(e_x) | is.infinite(e_x)] <- NA mu_x <- -c(diff(log(l_x)) / diff(edad), NA) data.frame(edad = edad, l_x = l_x, d_x = d_x, L_x = L_x, T_x = T_x, e_x = e_x, mu_x = mu_x) }) # ---------- Gráficos ---------- output$grafico_supervivencia <- renderPlot({ tabla <- calcular_tabla(); l0 <- tabla$l_x[1] ggplot(tabla, aes(x = edad, y = l_x / l0)) + geom_line(color = "#1a3a6b", linewidth = 1) + geom_vline(xintercept = input$edad, color = "red", linetype = "dashed") + labs(title = "Función de Supervivencia lₓ/l₀", x = "Edad (x)", y = "Proporción de supervivientes") + theme_minimal() + scale_y_continuous(labels = percent) }) output$grafico_esperanza <- renderPlot({ tabla <- calcular_tabla() ggplot(tabla, aes(x = edad, y = e_x)) + geom_line(color = "#22497e", linewidth = 1) + geom_vline(xintercept = input$edad, color = "red", linetype = "dashed") + geom_point(data = tabla[tabla$edad == input$edad, ], color = "red", size = 3) + labs(title = paste("Esperanza de Vida eₓ (e₀ =", round(tabla$e_x[1], 2), "años)"), x = "Edad (x)", y = "Esperanza de vida restante") + theme_minimal() }) output$grafico_mortalidad <- renderPlot({ tabla <- calcular_tabla() ggplot(tabla, aes(x = edad, y = mu_x)) + geom_line(color = "#1a3a6b", linewidth = 1) + geom_vline(xintercept = input$edad, color = "red", linetype = "dashed") + labs(title = "Tasa Instantánea de Mortalidad μₓ", x = "Edad (x)", y = "Fuerza de mortalidad") + theme_minimal() }) # ---------- Indicadores ---------- output$indicadores <- renderTable({ tabla <- calcular_tabla(); idx <- input$edad + 1 data.frame( Indicador = c( "Esperanza de vida al nacer (e₀)", paste("Esperanza de vida a los", input$edad, "años (eₓ)"), "Probabilidad de supervivencia a edad actual", "Fuerza de mortalidad a edad actual" ), Valor = c( round(tabla$e_x[1], 4), round(tabla$e_x[idx], 4), paste0(round(tabla$l_x[idx] / tabla$l_x[1] * 100, 2), "%"), round(tabla$mu_x[idx], 6) ) ) }, bordered = TRUE) # ---------- Tabla resumen ---------- output$tabla_resumen <- renderTable({ tabla <- calcular_tabla() edades <- unique(c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100)) edades <- edades[edades <= input$omega] res <- tabla[tabla$edad %in% edades, c("edad", "l_x", "e_x", "mu_x")] res$l_x <- round(res$l_x); res$e_x <- round(res$e_x, 2); res$mu_x <- round(res$mu_x, 6) names(res) <- c("Edad", "Sobrevivientes", "Esperanza vida", "μₓ") res }, bordered = TRUE) # ---------- Tabla Tₓ a la edad actual ---------- output$tabla_Tx <- renderTable({ tabla <- calcular_tabla() idx <- input$edad + 1 data.frame( Indicador = c( "Edad actual x", "Sobrevivientes lₓ", "Años-persona Tₓ = ∫ₓ^ω l(s) ds", "Esperanza de vida eₓ = Tₓ / lₓ" ), Valor = c( input$edad, round(tabla$l_x[idx], 2), round(tabla$T_x[idx], 2), round(tabla$e_x[idx], 4) ) ) }, bordered = TRUE) # ---------- Gráfico Tₓ frente a edad ---------- output$grafico_Tx <- renderPlot({ tabla <- calcular_tabla() ggplot(tabla, aes(x = edad, y = T_x)) + geom_line(color = "#1a3a6b", linewidth = 1) + geom_vline(xintercept = input$edad, color = "red", linetype = "dashed") + geom_point(data = tabla[tabla$edad == input$edad, ], color = "red", size = 3) + labs(title = "Años-persona por vivir Tₓ", x = "Edad (x)", y = "Tₓ") + theme_minimal() }) # ---------- Fórmula l(x) dinámica ---------- output$formula_lx <- renderUI({ if (input$modo == "dormoy1") { withMathJax(helpText( "$$l_x = l_0 \\cdot S^x$$", sprintf("con \\(l_0 = %s,\\; S = %s\\)", input$radix_d1, input$s_d1) )) } else if (input$modo == "dormoy2") { withMathJax(helpText( "$$l_x = l_0 \\cdot S_1^x \\cdot S_2^{x^2}$$", sprintf("con \\(l_0 = %s,\\; S_1 = %s,\\; S_2 = %s\\)", input$radix_d2, input$s1_d2, input$s2_d2) )) } else if (input$modo == "sang") { withMathJax(helpText( "$$l_x = l_0 \\cdot \\frac{B^x - B^{\\omega}}{1 - B^{\\omega}}$$", sprintf("con \\(l_0 = %s,\\; B = %s,\\; \\omega = %s\\)", input$radix_sang, input$b_sang, input$omega_sang) )) } else if (input$modo == "gompertz") { withMathJax(helpText( "$$l_x = l_0 \\cdot g^{c^x - 1}$$", "$$\\mu_x = -\\ln(g)\\cdot\\ln(c)\\cdot c^x$$", sprintf("con \\(l_0 = %s,\\; g = %s,\\; c = %s\\)", input$radix_gompertz, input$g_gompertz, input$c_gompertz) )) } else if (input$modo == "makeham") { withMathJax(helpText( "$$l_x = l_0 \\cdot s^x \\cdot g^{c^x - 1}$$", "$$\\mu_x = -\\ln(s) - \\ln(g)\\cdot\\ln(c)\\cdot c^x$$", sprintf("con \\(l_0 = %s,\\; s = %s,\\; g = %s,\\; c = %s\\)", input$radix_makeham, input$s_makeham, input$g_makeham, input$c_makeham) )) } else if (input$modo == "predef") { txt <- switch(input$modelo_predef, "cuadratica" = "l_x = l_0 \\cdot \\frac{10000 - x^2}{10000}", "moivre190" = "l_x = l_0 \\cdot \\frac{100 - x}{190}", "moivre150" = "l_x = l_0 \\cdot \\frac{100 - x}{150}", "lineal100" = "l_x = l_0 \\cdot \\frac{100 - x}{100}", "exponencial" = "l_x = l_0 \\cdot e^{-0{,}01x}" ) withMathJax(helpText(paste0("$$", txt, "$$"), sprintf("con \\(l_0 = %s\\)", input$radix_predef))) } else { HTML(paste0( "

Expresión introducida: l(x) = ", input$lx_expr, "

", "

", "Normalizada para que l(0) = l₀ = ", input$radix_custom, ".

" )) } }) # ---------- PASO 5 DINÁMICO: fórmula cerrada según modelo ---------- output$paso5_formula <- renderUI({ if (input$modo == "dormoy1") { S <- input$s_d1 ex_teorico <- -1 / log(S) withMathJax(HTML(paste0( "
$$e_x = \\frac{1}{\\mu} = -\\frac{1}{\\ln(S)}$$
", "

Con \\(S = ", S, "\\), resulta \\(e_x = ", round(ex_teorico, 4), "\\) años, ", "independiente de la edad (falta de memoria de la exponencial).

" ))) } else if (input$modo == "predef" && input$modelo_predef == "cuadratica") { x <- input$edad ex_teorico <- (100 - x) * (x + 200) / (3 * (100 + x)) withMathJax(HTML(paste0( "
$$e_x = \\frac{(100-x)(x+200)}{3(100+x)}$$
", "

Fórmula cerrada del modelo cuadrático \\(l_x = l_0(1 - x^2/10000)\\). ", "Con \\(x = ", x, "\\), resulta \\(e_x = ", round(ex_teorico, 4), "\\) años.

" ))) } else { withMathJax(HTML(paste0( "
$$e_x = \\int_0^{\\omega-x} {}_t p_x \\, dt$$
", "

Para este modelo no hay una fórmula cerrada simple. La esperanza se calcula ", "numéricamente como el área bajo la curva de supervivencia residual, o equivalentemente ", "con \\(\\int_0^{\\omega-x} t \\cdot {}_t p_x \\cdot \\mu_{x+t}\\,dt\\).

" ))) } }) } shinyApp(ui = ui, server = server)