---
# PASO 4: TÍTULO DEL CUADRO DE MANDOS
title: "Comparación de indicadores por país (fuente: OWID)"
output:
  flexdashboard::flex_dashboard:
    orientation: columns
runtime: shiny
---

<!--
================================================================================
  PLANTILLA DE ANÁLISIS DE LA RELACIÓN ENTRE DOS INDICADORES

  Estudia la relación entre dos atributos numéricos de una tabla de
  "observaciones x atributos": el diagrama de dispersión con su recta de
  regresión (primera pestaña), y los residuos de esa recta junto con la tabla
  del análisis del ajuste (segunda pestaña: recta, R2, sigma y los tres tests
  de aedv_lm_analisis, presentados con aedv_tabla_lm).

  Permite transformar cada indicador (raíz, logaritmo, Yeo-Johnson) para
  buscar una relación lineal donde no la hay entre las variables originales.

  Es una PLANTILLA: los indicadores del desplegable NO se enumeran a mano, se
  deducen de la propia tabla (todas sus columnas numéricas). Para aplicarla a
  otros datos solo hay que tocar los cuatro PASOS señalados; el resto —las
  opciones de los controles, los rótulos de los ejes y de la tabla, y el
  contenido de los tres paneles— se adapta solo.

  Requiere las funciones aedv_yeo_johnson(), aedv_lambda_optimo(),
  aedv_lm_analisis() y aedv_tabla_lm() del fichero utilidades.R del curso.
================================================================================
-->

```{r setup, include=FALSE}
library(flexdashboard) # Estructura del cuadro de mando en columnas
library(shiny)         # Reactividad
library(openxlsx)      # Lectura de Excel
library(plotly)        # Gráficas dinámicas: ggplotly convierte los ggplot
library(gt)            # Tabla de los tests coloreada (aedv_tabla_lm)
library(tidyverse)    # Manipulación de datos y ggplot2 (se carga el último
                       # para que dplyr::filter() no quede tapado por plotly)

source("https://ctim.es/AEDV/utilidades.R") # Funciones auxiliares del curso
```

```{r datos, include=FALSE}
# ==============================================================================
# PASO 1: LOS DATOS
#
# Se leen aquí, UNA SOLA VEZ al arrancar. Ha de ser una tabla en formato ANCHO
# de "observaciones x atributos": una fila por observación (país, municipio,
# provincia...) y una columna por indicador numérico.
# ==============================================================================
DATOS <- read.xlsx("https://ctim.es/AEDV/data/owid_country.xlsx", sheet = 1) |>
  as_tibble()

# ==============================================================================
# PASO 2: LA COLUMNA QUE IDENTIFICA CADA OBSERVACIÓN
#
# El nombre que se enseña al pasar el ratón por un punto y en la tabla de la
# derecha. Es la ÚNICA columna que hay que nombrar: los indicadores del
# desplegable son todas las columnas numéricas de DATOS, se deducen solas.
# ==============================================================================
COLUMNA_ETIQUETA <- "location"

# ==============================================================================
# PASO 3: OPCIONES DEL CUADRO DE MANDOS
# ==============================================================================

# Indicadores con los que arranca el cuadro de mandos. Si alguno no existe en
# DATOS se usa el primero disponible, así que la plantilla no se rompe al
# cambiar de datos y olvidarse de tocar estas dos líneas.
X_INICIAL <- "life_expectancy"
Y_INICIAL <- "gdp_per_capit"

# Nivel de significación con el que se decide si cada test se supera
ALFA <- 0.05
```

```{r preparacion, include=FALSE}
# ==============================================================================
# A PARTIR DE AQUÍ NO HAY QUE TOCAR NADA
# ==============================================================================

# LOS INDICADORES SON LAS COLUMNAS NUMÉRICAS de la tabla. Deducirlos así, y no
# por su posición, es lo que hace que la plantilla se adapte sola: da igual
# cuántas columnas de texto haya ni en qué orden estén.
INDICADORES <- DATOS |> select(where(is.numeric)) |> names()

# Valor inicial de un desplegable, con el primer indicador como red de seguridad
inicial <- function(nombre) if (nombre %in% INDICADORES) nombre else INDICADORES[1]

# Transformaciones ofrecidas para cada indicador. Se aplican conservando el
# signo del dato, para que valgan también con valores negativos.
# Yeo-Johnson es distinta y solo se ofrece en y: necesita un lambda, que se
# calcula buscando el mejor ajuste lineal con la x ya transformada.
TRANSFORMACIONES_X <- c("none", "sqrt", "log")
TRANSFORMACIONES_Y <- c("none", "sqrt", "log", "YeoJohnson")

transforma <- function(v, tipo) {
  switch(tipo,
    none = v,
    sqrt = sign(v) * sqrt(abs(v)),
    log  = sign(v) * log(abs(v)),
    v
  )
}

# Rótulo de un eje: el nombre del indicador y, entre paréntesis, lo que se le ha
# hecho: la normalización previa (si está marcada) y su transformación
rotulo <- function(indicador, tipo, normalizado = FALSE) {
  partes <- c(if (normalizado) "normalizado", if (tipo != "none") tipo)
  if (length(partes) == 0) indicador
  else paste0(indicador, " (", paste(partes, collapse = ", "), ")")
}
```

Column {.sidebar data-width=320}
--------------------------------------------------

```{r}
# NORMALIZACIÓN PREVIA. Aplica Yeo-Johnson a CADA variable POR SEPARADO, con el
# lambda que más aproxima esa variable a una normal (aedv_yeojohnson_lambda).
# Es un criterio DISTINTO del de la opción "YeoJohnson" del desplegable de y,
# que busca el lambda del mejor ajuste lineal ENTRE LAS DOS variables.
checkboxInput("normalizar",
              label = "Normalizar antes cada variable (Yeo-Johnson)",
              value = FALSE)

# Los desplegables se construyen con INDICADORES, deducido de los datos
selectInput("x", label = "x: indicador 1",
            choices = INDICADORES, selected = inicial(X_INICIAL))

selectInput("x_scale", label = "Transformación del indicador x",
            choices = TRANSFORMACIONES_X, selected = "none")

selectInput("y", label = "y: indicador 2",
            choices = INDICADORES, selected = inicial(Y_INICIAL))

selectInput("y_scale", label = "Transformación del indicador y",
            choices = TRANSFORMACIONES_Y, selected = "none")

# Lanza el cálculo: hasta pulsarlo no se calcula nada, ni al arrancar ni al
# cambiar los controles
actionButton("calcular", "Calcular")
```

```{r}
# CÁLCULO CENTRAL, del que beben los tres paneles: selecciona los dos
# indicadores, les aplica la transformación elegida, ajusta la recta de
# regresión y evalúa el ajuste. Devuelve una lista con nombres (no una lista
# posicional) para que cada panel tome solo lo que necesita.
#
# Es un eventReactive(): solo se recalcula al pulsar el botón «Calcular», no
# cada vez que se toca un control. Así se pueden cambiar varios controles
# seguidos y hacer el cálculo una sola vez. Al arrancar NO se calcula nada,
# para que el cuadro de mandos cargue antes: los paneles quedan vacíos hasta
# el primer clic.
calculo <- eventReactive(input$calcular, {

  tb <- DATOS |>
    select(etiqueta = all_of(COLUMNA_ETIQUETA),
           x = all_of(input$x),
           y = all_of(input$y)) |>
    na.omit()

  # PASO PREVIO OPCIONAL: normalizar cada variable por separado. Aquí solo
  # interviene la propia variable —no la otra—, y el lambda se busca para que
  # los datos transformados se parezcan lo más posible a una normal.
  lambda_norm_x <- NA_real_
  lambda_norm_y <- NA_real_
  if (isTRUE(input$normalizar)) {
    lambda_norm_x <- aedv_yeojohnson_lambda(tb$x)
    lambda_norm_y <- aedv_yeojohnson_lambda(tb$y)
    tb$x <- aedv_yeo_johnson(tb$x, lambda_norm_x)
    tb$y <- aedv_yeo_johnson(tb$y, lambda_norm_y)
  }

  tb$x <- transforma(tb$x, input$x_scale)

  # Yeo-Johnson calcula su lambda; el resto de transformaciones son directas
  lambda <- NA_real_
  if (input$y_scale == "YeoJohnson") {
    lambda <- aedv_lambda_optimo(tb$x, tb$y)
    tb$y   <- aedv_yeo_johnson(tb$y, lambda)
  } else {
    tb$y <- transforma(tb$y, input$y_scale)
  }

  # El logaritmo de un cero deja valores no finitos, que romperían la regresión
  tb <- tb |> filter(is.finite(x), is.finite(y))
  req(nrow(tb) >= 3) # con menos de 3 observaciones no hay nada que ajustar

  fit <- lm(y ~ x, data = tb)
  tb$residuo <- residuals(fit)

  # aedv_lm_analisis() ajusta esta misma recta y devuelve, en una fila, los
  # indicadores de calidad del ajuste y los p-valores de los tres tests.
  # Su test de Shapiro-Wilk admite como mucho 5.000 observaciones, así que con
  # tablas mayores se avisa en vez de dejar caer el cuadro de mandos.
  valores <- if (nrow(tb) <= 5000) aedv_lm_analisis(tb$x, tb$y) else NULL

  list(tb = tb, fit = fit, lambda = lambda, valores = valores,
       lambda_norm_x = lambda_norm_x, lambda_norm_y = lambda_norm_y)
})
```

###

```{r}
# Resumen del análisis en HTML. La recta, R2, sigma y los tests van en la tabla
# de la pestaña del error (aedv_tabla_lm); aquí queda lo que no está en ella.
renderUI({

  resultado <- calculo()
  tb        <- resultado$tb

  texto <- paste(
    "<strong>ANÁLISIS DE REGRESIÓN LINEAL</strong><br>",
    "<strong>x:</strong>", rotulo(input$x, input$x_scale, input$normalizar), "<br>",
    "<strong>y:</strong>", rotulo(input$y, input$y_scale, input$normalizar), "<br>",
    "<strong>observaciones:</strong>", nrow(tb),
    "<br><strong>correlación:</strong>", round(cor(tb$x, tb$y), digits = 4)
  )

  # Los lambda de la normalización previa: uno por variable, cada uno calculado
  # con la suya y sin mirar a la otra
  if (isTRUE(input$normalizar)) {
    texto <- paste(texto,
      "<br><strong>lambda normalidad (x):</strong>",
      formatC(resultado$lambda_norm_x, format = "e", digits = 3),
      "<br><strong>lambda normalidad (y):</strong>",
      formatC(resultado$lambda_norm_y, format = "e", digits = 3))
  }

  # El lambda del mejor ajuste lineal, que es otra cosa: se calcula con las dos
  # variables a la vez y solo existe si se ha elegido YeoJohnson en el eje y
  if (input$y_scale == "YeoJohnson") {
    texto <- paste(texto, "<br><strong>lambda ajuste lineal (y):</strong>",
                   formatC(resultado$lambda, format = "e", digits = 3))
  }

  HTML(texto)
})
```


Column {.tabset}
--------------------------------------------------

<!--
  LA COLUMNA CENTRAL VA EN DOS PESTAÑAS: el diagrama de dispersión, y el error
  del ajuste (los residuos) con la tabla del análisis de la recta. Las dos beben
  del mismo cálculo, así que cambiar de pestaña no vuelve a ajustar la recta.
-->

### Diagrama de dispersión

```{r}
# Diagrama de dispersión con la recta de regresión
renderPlotly({

  resultado <- calculo()

  p <- resultado$tb |>
    ggplot(aes(x, y, label = etiqueta)) +
      geom_point() +
      geom_smooth(method = lm, se = FALSE, formula = y ~ x) + # recta de regresión
      labs(x = rotulo(input$x, input$x_scale, input$normalizar),
           y = rotulo(input$y, input$y_scale, input$normalizar))

  ggplotly(p)
})
```

### Error del ajuste y tests

```{r}
# Residuos (error del ajuste) frente a x. Es la lectura gráfica de los tests de
# Shapiro-Wilk y Breusch-Pagan de la tabla de debajo: una nube sin forma y de
# anchura constante alrededor del cero es lo que se espera de un buen ajuste;
# un embudo o una curva delatan lo que esos tests detectan.
output$grafico_error <- renderPlotly({

  resultado <- calculo()

  p <- resultado$tb |>
    ggplot(aes(x, residuo, label = etiqueta)) +
      geom_hline(yintercept = 0, linetype = "dashed", colour = "steelblue") +
      geom_point(alpha = 0.6) +
      labs(x = rotulo(input$x, input$x_scale, input$normalizar), y = "residuo")

  ggplotly(p)
})

# La recta, R2, sigma y los tres tests de aedv_lm_analisis(), con el p-valor en
# verde si el test se supera y en rojo si no. aedv_tabla_lm() devuelve una tabla
# gt, que se muestra con renderUI.
output$tabla_lm <- renderUI({

  resultado <- calculo()
  validate(need(!is.null(resultado$valores),
                "Tests no disponibles: el test de Shapiro-Wilk admite como mucho 5.000 observaciones."))

  aedv_tabla_lm(resultado$valores, titulo = "Análisis del ajuste de la recta", alpha = ALFA)
})

# Dos salidas en la misma pestaña: la gráfica arriba, ocupando el alto que quede
# libre, y la tabla debajo con su altura natural (flex = NA). fillCol() es lo
# que reparte el alto de la pestaña entre las dos.
fillCol(flex = c(1, NA),
        plotlyOutput("grafico_error", height = "100%"),
        uiOutput("tabla_lm"))
```

Column {data-width=300}
--------------------------------------------------

### Valor de los indicadores transformados

```{r}
# Tabla de valores, ordenada de mayor a menor por el primer indicador
renderTable({

  tabla <- calculo()$tb |>
    arrange(desc(x)) |>
    transmute(
      etiqueta,
      x = prettyNum(round(x, digits = 2), big.mark = ","),
      y = prettyNum(round(y, digits = 2), big.mark = ",")
    )

  # Las cabeceras toman el nombre real de las columnas elegidas
  names(tabla) <- c(COLUMNA_ETIQUETA, input$x, input$y)
  tabla
})
```