Detección de cascadas terapéuticas asociadas a gabapentinoides y benzodiacepinas en adultos mayores de 65 años en España usando BIFAP mediante análisis de simetría de secuencia de prescripción

Versión de revisión v0.1 — no citar

Instantánea destinada a revisión interna. Para comentar, inicia sesión en Hypothesis con una cuenta que pertenezca al grupo privado GABALOOP. Esta web limita la barra lateral a ese grupo; el canal Public no está habilitado para esta revisión.

Introducción

La polifarmacia, comúnmente definida como el uso concomitante de 5 o más medicamentos, se ha asociado con múltiples problemas de la asistencia sanitaria, especialmente en los ancianos [1,2]. Estos incluyen un mayor riesgo de reacciones adversas a los medicamentos (RAM), pérdida de eficacia debido a una menor adherencia a la medicación y mayores costes sanitarios. El riesgo de polifarmacia aumenta con el número de enfermedades simultáneas que padece el paciente. En España, aproximadamente el 37.5% de la población tiene multimorbilidad (2 o más enfermedades crónicas), el 10.6% padece ≥5 enfermedades crónicas y el número promedio de problemas crónicos en adultos mayores (65-74 años) es de 2.8 [3]. Según la Base de Datos para la Investigación Farmacoepidemiológica en AP (BIFAP), la polifarmacia en la población española mayor de 14 años se ha triplicado en el período 2005-2015 pasando del 2.5 al 8.9%. En la franja de edad del 65 a 79 años esta prevalencia ha pasado del 9.1 al 23.4%; y en mayores de 80 años del 11.7 al 36.7%. El porcentaje de pacientes que consumía diez medicamentos o más aumentó diez veces en este mismo período [4].

Las cascadas terapéuticas (CT), una causa conocida de polifarmacia [5], ocurren cuando un nuevo medicamento se prescribe para tratar los signos o síntomas que surgen de un efecto adverso de otro fármaco, en lugar de suspender o reducir la dosis del medicamento que produce la RAM [6]. La facilidad en la identificación de las reacciones adversas puede dar lugar a CT intencionadas o no intencionadas [7], y que a su vez puedan ser o no apropiadas. Se han descrito muchos tipos de CT relacionadas con medicamentos utilizados en enfermedades crónicas, incluyendo cardiopatías, hipertensión, enfermedad pulmonar obstructiva, diabetes, demencia, enfermedad de Parkinson, epilepsia, esquizofrenia, depresión mayor y dolor crónico [8,9], como confirma una revisión paraguas y actualización sistemática reciente que sintetizó 190 estudios primarios sobre CT [10]. El riesgo de CT aumenta en los ancianos porque tienen más probabilidades de recibir múltiples terapias farmacológicas para múltiples enfermedades crónicas y porque las condiciones asociadas con el envejecimiento pueden enmascarar las posibles relaciones entre los medicamentos responsables y las RAM [8]. Una de las formas en las que se pueden identificar las CT es a través de la detección de medicación potencialmente inadecuada (MPI) y los potenciales efectos adversos producidos. Las MPI más frecuentemente detectadas en población anciana se relacionan con benzodiacepinas, gabapentinoides y otros depresores del sistema nervioso central, así como sus combinaciones [11–14].

Existen CT recientemente descritas: la prescripción de diuréticos de asa por el efecto paradójico de edemas en extremidades inferiores generado por los gabapentinoides; otra que también ha sido descrita es la prescripción de medicación antivertiginosa por la inestabilidad asociada a las benzodiacepinas [15–17]. La asociación gabapentinoides-diurético de asa (GABA-DA) se ha confirmado en al menos tres países: Canadá [16], Estados Unidos y Dinamarca [15], y, mediante un diseño ASSP que replica de forma conjunta ambas cascadas —el mismo abordaje empleado en este estudio—, en Japón, donde también se describió la cascada benzodiacepina-antivertiginoso [18]. Aunque esto se ha explorado cuando los gabapentinoides se han prescrito en el dolor lumbar, no se ha estudiado en aquellos que lo usen de forma generalizada u otras formas como dolor neuropático de localizaciones diversas, ni con otros medicamentos de la misma clase aprobados en España.

El observatorio del uso de medicamentos en España reconoce la importancia de la creciente prescripción de ansiolíticos e hipnóticos y sedantes, así como de determinados antiepilépticos como pregabalina [19,20]. El aumento de la conciencia sobre las CT entre los profesionales a través de la evidencia derivada de la investigación es crucial para su prevención [8,21], por lo que conocer el alcance y potencial impacto de estas en la población española sería de especial importancia. La base de datos de BIFAP es potencialmente un entorno donde analizar este tipo de estudios de farmacovigilancia dada la potencial baja frecuencia de algunos eventos [22].

El objetivo de este estudio es estimar las tasas crudas y ajustadas de prescripción secuencial de las combinaciones de diurético tras gabapentinoides (GABA-DA) y de antivertiginosos tras benzodiacepinas (BZD-AV) en población española mayor de 65 años con polifarmacia, comparada con el resto de la población incluida en BIFAP; como objetivo secundario, se describieron las características sociodemográficas de los sujetos incluidos y se estudiaron los factores individuales asociados al orden temporal de las prescripciones.

Material y Métodos

Diseño

Estudio longitudinal multicéntrico retrospectivo de análisis de secuencias simétricas de prescripción (ASSP). El abordaje ASSP es un diseño que incluye solo sujetos que presentan el evento de interés (presencia de fármaco considerado evento) y se ha usado en más de 50 estudios para detectar prescripción en cascada [23–25]. Entre sus ventajas se encuentra que es un diseño más sencillo, de más fácil computación y robusto para su ajuste de factores confusores que dependen del tiempo. El análisis ASSP valora la temporalización relativa del inicio de la exposición al fármaco (sospechoso de causar la RAM) frente al fármaco evento en pacientes que reciben ambos fármacos incidentes, en lugar de usuarios prevalentes (Supplementary Picture 1).

Población y muestra

Se incluyeron los pacientes con más de 360 días de registro continuado en BIFAP antes de la primera prescripción del fármaco de exposición, que definió la fecha de entrada en la cohorte, y con al menos 180 días de seguimiento posterior. Se exigió un período de 180 días sin exposición previa. El inicio del fármaco resultado/evento se aceptó en los 90 días anteriores o posteriores a la entrada en la cohorte. Las prescripciones de exposición y evento en el mismo día se excluyeron del ratio de secuencia y se informaron por separado (Supplementary Picture 1).

Variables

1. Principal: se determinó la cohorte de exposición a partir del listado de fármacos de exposición (gabapentinoides y benzodiacepinas) y los eventos a partir del listado de fármacos clasificados como eventos (diuréticos de asa y antivertiginosos) en los períodos especificados en el diseño (Supplementary Picture 1). Se utilizó la fecha de prescripción como referencia para los cálculos.
2. Secundarias y de ajuste: sociodemográficas (id anonimizado, sexo, edad); clínicas (número de fármacos crónicos —tomados durante más de 3 meses continuados previos a la fecha de inclusión—, presencia de determinados fármacos confusores en una ventana de 90 días pre-post al momento de entrada, número de fármacos crónicos al finalizar el seguimiento a 180 días, y comorbilidades del listado O’Halloran con su fecha de episodio); codificación CIAP asociada al fármaco de exposición y su fecha de apertura; y, del evento, codificación CIAP asociada, presencia del código K07 posterior a la entrada en la cohorte, y fecha de nuevos episodios codificados que incluyan “EDEMA” en el descriptivo en los tres meses posteriores a la entrada (sujeta a verificación en historia clínica y validación por al menos un investigador del proyecto).
3. De ajuste multinivel: comunidad autónoma (CCAA) del centro de prescripción, obtenida a partir de los datos de Dataset3. El identificador de cupo médico (CIAS) no está disponible en los extractos facilitados por BIFAP.

Fuente de datos

Se incluyeron los registros disponibles de atención primaria del territorio nacional en BIFAP, una base de datos que contiene información anónima sobre diagnósticos clínicos, notas médicas en texto libre, pruebas de laboratorio, exploraciones complementarias y prescripciones (incluida dispensación), registrada por médicos de atención primaria del Sistema Nacional de Salud y validada a través de múltiples estudios de investigación farmacoepidemiológica [22]. BIFAP incluye 20 millones de historias clínicas con 2.400 millones de prescripciones [26].

Análisis estadístico

Las características sociodemográficas y clínicas se describieron mediante mediana y rango intercuartílico para las variables continuas, y frecuencia y porcentaje para las variables categóricas.

Para el análisis de secuencia simétrica de prescripción (ASSP) se calculó el ratio de secuencia (SR) como medida de asimetría entre nuevos usuarios de ambos fármacos [27]. Para cada paciente biexpuesto (que recibió tanto el fármaco exposición como el fármaco evento en la ventana [–90, +90] días), se identificó la primera prescripción cronológica del fármaco evento dentro de esa ventana —es decir, la de fecha calendario más temprana (valor algebraico mínimo de los días transcurridos desde el índice), y no necesariamente la de menor distancia absoluta al índice—, tras eliminar duplicados exactos. Esta distinción es relevante: un paciente con una prescripción en el día –20 y otra en el día +5 se clasifica como B→A (pre), no como A→B (post), aunque la segunda prescripción esté más próxima a la fecha índice. Los pacientes cuyo inicio observable coincidió con la fecha índice (mismo día, día 0) se excluyeron del SR y se informaron por separado. El SR se obtuvo dividiendo el número de pacientes cuyo evento se inició después de la exposición (secuencia A→B) entre aquellos cuyo evento se inició antes (secuencia B→A). Se excluyeron como usuarios prevalentes quienes tenían una prescripción del evento en los 91–270 días anteriores al índice.

El SR crudo se calculó globalmente y por subgrupo (Polimedicado ≥65 años frente a Resto), con intervalos de confianza al 95% mediante aproximación de Wald en escala logarítmica. También se calculó la razón entre SR de ambos subgrupos (RSR). La corrección por tendencia secular mediante el Null-Effect Sequence Ratio (NESR) del método de Hallas [23] requiere las tasas de prescripción del fármaco evento en la población total de BIFAP por período calendario, que no forman parte de los extractos facilitados para este proyecto y cuya obtención no está prevista. Por ello se fijó NESR = 1 (sin corrección de tendencia) con carácter definitivo, no como un valor provisional a la espera de datos adicionales, y todos los resultados se denominan SR crudos.

Entre los pacientes biexpuestos clasificables como pre o post se ajustó una regresión logística multinivel, con intercepto aleatorio por comunidad autónoma (CCAA; CIAS no disponible). El desenlace fue secuencia post frente a pre. La comparación principal fue Polimedicado ≥65 frente a Resto, ajustada por sexo, categoría de edad, número de fármacos crónicos, número de enfermedades crónicas, fármacos confusores y enfermedades específicas de cada cascada. Los OR describen el orden temporal entre biexpuestos y no la incidencia ni un efecto causal. La heterogeneidad geográfica se resumió con la mediana de oportunidad relativa (MOR) y el coeficiente de correlación intraclase (ICC). Como validación post-hoc de la especificación, se comparó mediante AIC/BIC el modelo completo frente a versiones sin las enfermedades específicas de indicación y sin el efecto aleatorio por CCAA, y se contrastó este último mediante una prueba de razón de verosimilitudes (LRT); la discriminación se resumió con el área bajo la curva ROC (AUC) y el pseudo-R² de Nagelkerke.

Se realizaron análisis de sensibilidad con ventanas de ±30 y ±60 días, exclusión de pacientes con eventos pre y post, clasificación del mismo día como pre o post, estratificación por ATC y año calendario, y una especificación alternativa edad×polifarmacia. Todos los análisis se realizaron con R 4.6.0. Codex (OpenAI) se utilizó para la revisión asistida y depuración de los scripts; los investigadores verificaron el código, las decisiones metodológicas y los resultados. Git se utilizó para control de versiones y trazabilidad del pipeline analítico. El código y los materiales reproducibles se mantienen en GitHub (https://github.com/juanchingao/2025-07-GABALOOP.git). Los datos individuales de BIFAP no forman parte del repositorio por sus condiciones de acceso restringido y podrán solicitarse por los cauces establecidos por BIFAP, previa autorización correspondiente.

Plan de datos

Los datos fueron anonimizados en origen y se recibieron en ficheros protegidos a través del circuito seguro de explotación de datos de BIFAP. Se almacenaron y analizaron en el servidor seguro del Laboratorio de Supercomputación de la Universidad Rey Juan Carlos. Se previó su destrucción conforme al protocolo de BIFAP una vez completados los objetivos.

Consideraciones éticas

El estudio se realizó de acuerdo con los principios básicos de protección de los derechos y la dignidad del ser humano según la Declaración de Helsinki (Fortaleza 2013) y las directrices nacionales vigentes, el Reglamento (UE) 2016/679 de Protección de Datos (RGPD) y la Ley Orgánica 3/2018 de Protección de Datos de Carácter Personal y garantía de derechos digitales. El protocolo contempló la solicitud de exención del consentimiento informado por el uso retrospectivo de datos de vida real, de acuerdo con los requisitos del comité ético de investigación clínica del Hospital Clínico San Carlos.

Resultados

In [1]:

find_project_root <- function(start = getwd()) {
  path <- normalizePath(start, winslash = "/", mustWork = FALSE)
  for (i in 1:5) {
    if (dir.exists(file.path(path, "scripts")) && dir.exists(file.path(path, "config"))) {
      return(path)
    }
    parent <- dirname(path)
    if (identical(parent, path)) break
    path <- parent
  }
  stop("No se pudo localizar la raiz del proyecto GABALOOP desde ", start)
}

project_root <- find_project_root()
source(file.path(project_root, "config", "paths.R"))
outputs_dir <- get_outputs_dir(project_root)

# Las figuras externas deben copiarse a manuscript/figures/ para que Quarto
# pueda calcular rutas relativas correctas en el HTML/DOCX de salida.
# (knitr::include_graphics con ruta absoluta produce ../outputs/figures/...
#  que desde _manuscript/ resuelve mal; la copia local lo evita.)
local({
  fig_dir <- file.path(project_root, "manuscript", "figures")
  dir.create(fig_dir, showWarnings = FALSE)
  src <- file.path(outputs_dir, "figures", "figura_flujo_poblacion.png")
  if (file.exists(src))
    file.copy(src, file.path(fig_dir, "figura_flujo_poblacion.png"), overwrite = TRUE)
  src1 <- file.path(outputs_dir, "figures", "figura1_cohorte.png")
  if (file.exists(src1))
    file.copy(src1, file.path(fig_dir, "figura1_cohorte.png"), overwrite = TRUE)
  src2 <- file.path(outputs_dir, "figures", "sr_anual_ambas_cascadas.png")
  if (file.exists(src2))
    file.copy(src2, file.path(fig_dir, "sr_anual_ambas_cascadas.png"), overwrite = TRUE)
})

[1] TRUE

Selección de la población de estudio

In [2]:

knitr::include_graphics("figures/figura_flujo_poblacion.png")

Figure 1: Diagrama de flujo de la población de estudio (Dataset1)

Del extracto original de Dataset1 (BIFAP, N = 6.173.351), se excluyeron 1.217 pacientes con fecha de nacimiento no plausible (anterior a 1910, edad >110 años al entrar en cohorte) y 138.358 pacientes expuestos a otros antiepilépticos no incluidos en la definición analítica de exposición (mayoritariamente topiramato y levetiracetam; Figure 1), quedando una población de estudio final de 6.033.776 pacientes.

Características de la población incluida

In [3]:

library(gt)
library(readr)

tabla1 <- read_csv(
  file.path(outputs_dir, "tables", "tabla1_descriptiva_poblacion.csv"),
  show_col_types = FALSE
)

tbl_gt <- tabla1 |>
  gt() |>
  cols_label(
    `Característica`         = md("**Característica**"),
    `Global (N = 6.033.776)` = md("**Global**<br>(N = 6.033.776)"),
    `Benzodiacepina`         = md("**Benzodiacepina**<br>(n = 5.431.812)"),
    `Gabapentinoide`         = md("**Gabapentinoide**<br>(n = 601.964)")
  ) |>
  tab_style(
    style = list(cell_fill(color = "#eef2f8"), cell_text(weight = "bold")),
    locations = cells_body(rows = Característica == "N")
  ) |>
  tab_footnote(
    footnote = "Los valores continuos se expresan como mediana [P25–P75]; los categóricos como n (%).",
    locations = cells_column_labels(columns = `Característica`)
  ) |>
  tab_footnote(
    footnote = "N05 y N03 incluyen los fármacos de exposición del estudio (benzodiacepinas y gabapentinoides, respectivamente); sus prevalencias reflejan la definición de la cohorte, no comorbilidad.",
    locations = cells_body(
      columns = `Característica`,
      rows = grepl("^(N05|N03) —", Característica)
    )
  ) |>
  opt_footnote_marks(marks = "standard") |>
  tab_options(
    table.border.top.color            = "#2b3a55",
    column_labels.border.bottom.color = "#2b3a55",
    table_body.hlines.color           = "#dee2e6",
    table.width                       = pct(100)
  )

# knitr::asis_output() marca el resultado como HTML crudo para que knitr
# no lo envuelva en <pre><code>; pandoc lo procesa directamente.
# En DOCX pandoc convierte la tabla HTML a tabla Word nativa.
knitr::asis_output(as_raw_html(tbl_gt))

In [4]:

Table 1: Características descriptivas de la población de estudio, por grupo de exposición

Característica*	Global (N = 6.033.776)	Benzodiacepina (n = 5.431.812)	Gabapentinoide (n = 601.964)
N	6.033.776	5.431.812	601.964
Sexo, Mujer	3,528,289 (58.5%)	3,226,541 (59.4%)	301,748 (50.1%)
Sexo, Hombre	2,505,487 (41.5%)	2,205,271 (40.6%)	300,216 (49.9%)
Edad al entrar en cohorte (años)	50 [37, 64]	49 [37, 63]	59 [46, 72]
Tiempo en BIFAP (años)	14.4 [10.0, 18.0]	14.5 [10.0, 18.0]	13.7 [9.0, 17.8]
Tiempo desde inicio de registro hasta entrada en cohorte (años)	4.1 [1.7, 8.0]	4.0 [1.6, 7.8]	4.7 [2.1, 9.2]
Nº de fármacos crónicos	5 [3, 9]	5 [3, 8]	7 [5, 11]
Nº de enfermedades crónicas	2.00 [1.00, 4.00]	2.00 [1.00, 4.00]	3.00 [1.00, 5.00]
Enfermedades crónicas seleccionadas, n (%)	NA	NA	NA
Hipertensión arterial	1.550.280 (25,7%)	1.336.466 (24,6%)	213.814 (35,5%)
Cefalea (otras)	976.741 (16,2%)	891.725 (16,4%)	85.016 (14,1%)
Trastorno de ansiedad	935.141 (15,5%)	901.101 (16,6%)	34.040 (5,7%)
Obesidad	720.848 (11,9%)	625.736 (11,5%)	95.112 (15,8%)
Depresión	667.022 (11,1%)	621.788 (11,4%)	45.234 (7,5%)
Diabetes mellitus (cualquier tipo)	525.273 (8,7%)	429.708 (7,9%)	95.565 (15,9%)
Cáncer de piel	502.327 (8,3%)	443.640 (8,2%)	58.687 (9,7%)
Anemia (otras)	494.713 (8,2%)	445.434 (8,2%)	49.279 (8,2%)
Varices	496.541 (8,2%)	433.847 (8,0%)	62.694 (10,4%)
Síndrome lumbar con ciática	495.911 (8,2%)	392.898 (7,2%)	103.013 (17,1%)
Arritmia (otras)	347.384 (5,8%)	305.844 (5,6%)	41.540 (6,9%)
Grupos ATC-2 seleccionados en fármacos crónicos, n (%)	NA	NA	NA
N05 — Psicolépticos (ansiolíticos, hipnóticos, antipsicóticos)†	5.543.977 (91,9%)	5.431.810 (100,0%)	112.167 (18,6%)
M01 — Antiinflamatorios no esteroideos	3.248.624 (53,8%)	2.884.881 (53,1%)	363.743 (60,4%)
N02 — Analgésicos	3.128.277 (51,8%)	2.704.567 (49,8%)	423.710 (70,4%)
A02 — Antiácidos y antiulcerosos	2.308.965 (38,3%)	1.988.069 (36,6%)	320.896 (53,3%)
N06 — Psicoanalépticos (antidepresivos, etc.)	1.270.676 (21,1%)	1.171.342 (21,6%)	99.334 (16,5%)
C09 — Agentes del sistema renina-angiotensina	1.223.713 (20,3%)	1.047.241 (19,3%)	176.472 (29,3%)
C10 — Hipolipemiantes	1.130.713 (18,7%)	956.427 (17,6%)	174.286 (29,0%)
B01 — Antitrombóticos	851.172 (14,1%)	721.425 (13,3%)	129.747 (21,6%)
N03 — Antiepilépticos (gabapentinoides incluidos)†	787.218 (13,0%)	185.254 (3,4%)	601.964 (100,0%)
R03 — Antiobstructivos respiratorios	696.202 (11,5%)	612.735 (11,3%)	83.467 (13,9%)
C03 — Diuréticos	523.387 (8,7%)	449.526 (8,3%)	73.861 (12,3%)
* Los valores continuos se expresan como mediana [P25–P75]; los categóricos como n (%).
† N05 y N03 incluyen los fármacos de exposición del estudio (benzodiacepinas y gabapentinoides, respectivamente); sus prevalencias reflejan la definición de la cohorte, no comorbilidad.

De los 6.033.776 pacientes de la población de estudio (Table 1), 5.431.812 (90,0%) estaban expuestos a benzodiacepinas y 601.964 (10,0%) a gabapentinoides. El grupo de gabapentinoides presentó una mayor edad al entrar en cohorte (mediana 59 [P25-P75: 46-72] años frente a 49 [37-63] en benzodiacepinas) y una distribución de sexo más equilibrada (49,9% hombres frente a 40,6% en benzodiacepinas). Asimismo, presentó una mayor carga de polifarmacia y comorbilidad: mediana de 7 [5-11] fármacos crónicos frente a 5 [3-8] en benzodiacepinas, y 3 [1-5] enfermedades crónicas frente a 2 [1-4], respectivamente.

Cohortes ASSP

In [5]:

local({
  fig_dir <- file.path(project_root, "manuscript", "figures")
  dir.create(fig_dir, showWarnings = FALSE)
  for (fn in c("assp_forest_sr.png", "sr_anual_ambas_cascadas.png", "modelo_or_forest.png",
               "assp_semanal_gaba_da.png", "assp_semanal_bzd_av.png")) {
    src <- file.path(outputs_dir, "figures", fn)
    if (file.exists(src))
      file.copy(src, file.path(fig_dir, fn), overwrite = TRUE)
  }
})

tabla_assp <- readr::read_csv(
  file.path(outputs_dir, "tables", "tabla_assp.csv"),
  show_col_types = FALSE
)

Para la cascada GABA-DA, la cohorte de análisis incluyó 572.111 pacientes expuestos a gabapentinoides tras excluir 195.820 usuarios prevalentes de diurético de asa (con prescripción en los 91–270 días previos al índice). De ellos, 14.908 (2,6%) presentaron al menos una prescripción de diurético de asa en la ventana [–90, +90] días. Para la cascada BZD-AV, la cohorte incluyó 5.326.401 pacientes expuestos a benzodiacepinas (excluidos 127.247 usuarios prevalentes de antivertiginoso), de los que 144.168 (2,7%) presentaron al menos una prescripción de antivertiginoso en la ventana.

Análisis de secuencias simétricas de prescripción

In [6]:

library(gt)

tbl_assp_gt <- tabla_assp |>
  dplyr::select(cascada, subgrupo, n_pre, n_post, n_mismo_dia, sr_crudo, ic_low, ic_high) |>
  gt(groupname_col = "cascada") |>
  cols_label(
    subgrupo    = md("**Subgrupo**"),
    n_pre       = md("**Pre**"),
    n_post      = md("**Post**"),
    n_mismo_dia = md("**Mismo día**"),
    sr_crudo    = md("**SR crudo**"),
    ic_low      = md("**IC 95% inf.**"),
    ic_high     = md("**IC 95% sup.**")
  ) |>
  fmt_number(columns = c(n_pre, n_post, n_mismo_dia), sep_mark = ".", dec_mark = ",", decimals = 0) |>
  fmt_number(columns = c(sr_crudo, ic_low, ic_high), dec_mark = ",", decimals = 2) |>
  tab_footnote(
    "NESR = 1 (sin corrección por tendencia secular; limitación permanente, ver Métodos). IC 95% por aproximación de Wald en escala logarítmica.",
    locations = cells_column_labels(columns = sr_crudo)
  ) |>
  tab_options(table.width = pct(100))

knitr::asis_output(as_raw_html(tbl_assp_gt))

In [7]:

Table 2: Ratio de secuencia (SR crudo) por cascada y subgrupo. Pre: pacientes con primera prescripción del fármaco evento antes del índice; post: después del índice; mismo día: excluidos del SR.

Subgrupo	Pre	Post	Mismo día	SR crudo1	IC 95% inf.	IC 95% sup.
GABA-DA
Global	5.200	6.251	3.457	1,20	1,16	1,25
Polimedicado ≥65	4.043	4.865	2.723	1,20	1,15	1,25
Resto	1.157	1.386	734	1,20	1,11	1,29
BZD-AV
Global	49.397	52.898	41.873	1,07	1,06	1,08
Polimedicado ≥65	17.569	19.247	10.922	1,10	1,07	1,12
Resto	31.828	33.651	30.951	1,06	1,04	1,07
1 NESR = 1 (sin corrección por tendencia secular; limitación permanente, ver Métodos). IC 95% por aproximación de Wald en escala logarítmica.

En GABA-DA se clasificaron 5.200 pacientes como pre y 6.251 como post, con SR crudo global de 1,20 (IC 95%: 1,16–1,25). En BZD-AV hubo 49.397 secuencias pre y 52.898 post, con SR de 1,07 (IC 95%: 1,06–1,08) (Table 2).

La distribución semanal de las primeras prescripciones del fármaco evento en torno a la fecha índice mostró un incremento progresivo tras la exposición en ambas cascadas, sin un salto abrupto aislado que sugiriese un artefacto de registro puntual (Figure 2, Figure 3).

In [8]:

knitr::include_graphics("figures/assp_semanal_gaba_da.png")

Figure 2: ASSP GABA-DA: número de prescripciones de diurético de asa por semana relativa a la fecha índice. La semana +0 (mismo día) se excluye del SR y se informa por separado.

In [9]:

knitr::include_graphics("figures/assp_semanal_bzd_av.png")

Figure 3: ASSP BZD-AV: número de prescripciones de antivertiginoso por semana relativa a la fecha índice. La semana +0 (mismo día) se excluye del SR y se informa por separado.

En GABA-DA los SR fueron prácticamente idénticos en Polimedicado ≥65 (SR 1,20; IC 95%: 1,15–1,26) y Resto (SR 1,20; IC 95%: 1,11–1,30); la RSR fue 1,00 (IC 95%: 0,92–1,10). En BZD-AV se observó una diferencia pequeña: SR 1,10 (IC 95%: 1,07–1,12) frente a 1,06 (IC 95%: 1,04–1,07), con RSR 1,04 (IC 95%: 1,01–1,06) (Figure 4).

In [10]:

knitr::include_graphics("figures/assp_forest_sr.png")

Figure 4: Ratio de secuencia crudo (SR) e intervalos de confianza al 95% por cascada y subgrupo. La línea discontinua indica SR = 1 (sin asimetría de secuencia). NESR = 1 (sin corrección por tendencia secular).

Se clasificaron como mismo día 3.457 pacientes (23,2% de los biexpuestos) en GABA-DA y 41.873 (29,0%) en BZD-AV. Estos casos se excluyeron del análisis principal y se exploraron mediante escenarios alternativos, dada la imposibilidad de establecer su orden temporal.

Factores asociados a la cascada terapéutica

Se ajustó un modelo de regresión logística multinivel con efecto aleatorio por CCAA (9 comunidades autónomas) para cada cascada (Figure 5). La variable respuesta fue el inicio del fármaco evento después del fármaco de exposición (y = 1: secuencia A→B) frente a antes (y = 0: secuencia B→A); se excluyeron las prescripciones en mismo día. Se incluyeron como covariables fijas el subgrupo de exposición principal, sexo, categoría de edad, número de fármacos crónicos (por unidad adicional), número de enfermedades crónicas (por unidad adicional), presencia de fármacos confusores preespecificados, y confusores por indicación específicos de cada cascada.

In [11]:

library(gt)
library(readr)

tabla_mod_csv <- read_csv(
  file.path(outputs_dir, "tables", "tabla_modelo_logistico.csv"),
  show_col_types = FALSE
)

nombres_es <- c(
  "(Intercept)"                              = "Intercepto",
  "subgrupoPolimedicado \u226565"            = "Polimedicado \u226565 (vs. Resto)",
  "sexoM"                                    = "Sexo masculino (vs. femenino)",
  "edad_cat65\u201374"                       = "Edad 65\u201374 (vs. <65)",
  "edad_cat75\u201384"                       = "Edad 75\u201384 (vs. <65)",
  "edad_cat\u226585"                         = "Edad \u226585 (vs. <65)",
  "n_farmacos_cronicos (por f\u00e1rmaco)"  = "F\u00e1rmacos cr\u00f3nicos (por unidad)",
  "n_enfermedades_cronicas (por enfermedad)" = "Enfermedades cr\u00f3nicas (por unidad)",
  "fcosconf_var"                             = "F\u00e1rmaco confusor (s\u00ed vs. no)",
  "IN_RENCRON_VAR"  = "Insuficiencia renal cr\u00f3nica",
  "K77_INSCAR_VAR"  = "Insuficiencia card\u00edaca",
  "DIABETES_VAR"    = "Diabetes mellitus",
  "K87_HTA_VAR"     = "Hipertensi\u00f3n arterial",
  "N94_NEUROP_VAR"  = "Neuropat\u00eda perif\u00e9rica",
  "H82_SDVERT_VAR"  = "S\u00edndrome vertiginoso",
  "P74_ANSIED_VAR"  = "Trastorno de ansiedad",
  "P70_DEMENC_VAR"  = "Demencia",
  "P76_DEPRES_VAR"  = "Depresi\u00f3n",
  "N87_PARKIN_VAR"  = "Enfermedad de Parkinson",
  "P72_ESQUIZ_VAR"  = "Esquizofrenia/psicosis"
)

tbl_mod_gt <- tabla_mod_csv |>
  dplyr::filter(effect == "fixed", term != "(Intercept)") |>
  dplyr::mutate(
    Term = dplyr::recode(term, !!!nombres_es),
    OR   = sprintf("%.2f", estimate),
    IC95 = sprintf("%.2f\u2013%.2f", conf.low, conf.high),
    p    = dplyr::case_when(
      p.value < 0.001 ~ "<0,001",
      TRUE            ~ sprintf("%.3f", p.value)
    )
  ) |>
  dplyr::select(cascada, Term, OR, IC95, p) |>
  gt(groupname_col = "cascada") |>
  cols_label(
    Term = md("**Variable**"),
    OR   = md("**OR**"),
    IC95 = md("**IC 95%**"),
    p    = md("***p***")
  ) |>
  tab_style(
    style = list(cell_text(weight = "bold")),
    locations = cells_body(
      columns = OR,
      rows    = p %in% c("<0,001") | as.numeric(gsub(",",".",p)) < 0.05
    )
  ) |>
  tab_options(table.width = pct(100),
              table.border.top.color = "#2b3a55",
              column_labels.border.bottom.color = "#2b3a55",
              table_body.hlines.color = "#dee2e6")

Warning in rlang::eval_tidy(expr = quo, data = data): NAs introduced by
coercion

knitr::asis_output(as_raw_html(tbl_mod_gt))

In [12]:

Table 3: Odds ratios ajustados (IC 95%) del modelo de regresión logística multinivel. Referencia: subgrupo Resto; sexo femenino; edad <65 años.

Variable	OR	IC 95%	p
GABA-DA
Polimedicado ≥65 (vs. Resto)	0.85	0.53–1.37	0.503
Sexo masculino (vs. femenino)	0.96	0.89–1.04	0.333
Edad 65–74 (vs. <65)	1.23	0.76–2.00	0.396
Edad 75–84 (vs. <65)	1.36	0.84–2.21	0.206
Edad ≥85 (vs. <65)	1.36	0.84–2.21	0.216
Fármacos crónicos (por unidad)	1.01	1.01–1.02	<0,001
Enfermedades crónicas (por unidad)	0.98	0.97–1.00	0.045
Fármaco confusor (sí vs. no)	1.07	0.99–1.17	0.092
Insuficiencia renal crónica	0.76	0.67–0.87	<0,001
Insuficiencia cardíaca	0.40	0.35–0.47	<0,001
Diabetes mellitus	0.92	0.84–1.01	0.065
Hipertensión arterial	0.93	0.85–1.01	0.082
Neuropatía periférica	0.84	0.62–1.13	0.247
BZD-AV
Polimedicado ≥65 (vs. Resto)	1.02	0.93–1.13	0.633
Sexo masculino (vs. femenino)	1.09	1.07–1.12	<0,001
Edad 65–74 (vs. <65)	0.93	0.84–1.02	0.122
Edad 75–84 (vs. <65)	0.95	0.86–1.05	0.298
Edad ≥85 (vs. <65)	1.10	0.99–1.22	0.065
Fármacos crónicos (por unidad)	1.02	1.01–1.02	<0,001
Enfermedades crónicas (por unidad)	0.98	0.98–0.99	<0,001
Fármaco confusor (sí vs. no)	0.98	0.94–1.01	0.189
Síndrome vertiginoso	0.49	0.47–0.52	<0,001
Trastorno de ansiedad	0.99	0.95–1.02	0.513
Demencia	1.17	1.09–1.26	<0,001
Depresión	1.12	1.07–1.16	<0,001
Enfermedad de Parkinson	0.89	0.65–1.22	0.460
Esquizofrenia/psicosis	0.89	0.80–0.98	0.017

In [13]:

knitr::include_graphics("figures/modelo_or_forest.png")

Figure 5: Odds ratios ajustados (IC 95%) del modelo de regresión logística multinivel ASSP, faceteado por cascada (izquierda: GABA-DA; derecha: BZD-AV; cada una con su propio modelo y confusores por indicación). Escala logarítmica; la línea discontinua indica OR = 1. Punto relleno: IC 95% excluye 1; punto hueco: IC 95% incluye 1. Efecto aleatorio por CCAA (9 comunidades autónomas).

Cascada GABA-DA. El subgrupo Polimedicado ≥65 no se asoció con el orden post frente a pre (OR = 0,85; IC 95%: 0,53–1,37; p = 0,503) (Table 3, Figure 5). Cada fármaco crónico adicional se asoció con un OR de 1,01 (IC 95%: 1,01–1,02), mientras que cada enfermedad crónica adicional presentó OR 0,98 (IC 95%: 0,97–1,00). La insuficiencia cardíaca (OR = 0,40; IC 95%: 0,35–0,47) y la insuficiencia renal crónica (OR = 0,76; IC 95%: 0,67–0,87) se asociaron con mayor frecuencia de secuencia pre. La heterogeneidad residual por CCAA fue pequeña (MOR = 1,07; ICC = 0,0014).

Cascada BZD-AV. Tampoco hubo asociación ajustada del subgrupo principal (OR = 1,02; IC 95%: 0,93–1,13; p = 0,633). El síndrome vertiginoso crónico se asoció con secuencia pre (OR = 0,49; IC 95%: 0,47–0,52). La demencia (OR = 1,17; IC 95%: 1,09–1,26), la depresión (OR = 1,12; IC 95%: 1,07–1,16), el sexo masculino (OR = 1,09; IC 95%: 1,07–1,12) y el número de fármacos crónicos (OR = 1,02 por unidad; IC 95%: 1,01–1,02) se asociaron con mayor frecuencia de secuencia post. La heterogeneidad por CCAA también fue reducida (MOR = 1,07; ICC = 0,0016).

En ambas cascadas, la pequeña diferencia cruda entre subgrupos no persistió tras el ajuste. Estos OR describen el orden temporal entre pacientes biexpuestos y no deben interpretarse como riesgo individual ni efecto causal.

Análisis de sensibilidad

Los resultados variaron con la definición temporal. Con ventanas de ±30 días, el SR fue 0,86 en ambas cascadas; con ±60 días fue 1,03 en GABA-DA y 0,97 en BZD-AV; y con ±90 días alcanzó 1,20 y 1,07, respectivamente. Clasificar las coprescripciones del día índice como post elevó los SR a 1,87 y 1,92, mientras que clasificarlas como pre los redujo a 0,72 y 0,58.

En BZD-AV se observó heterogeneidad según el fármaco evento: el grupo N07CA presentó un SR de 0,91 (IC 95%: 0,90–0,92), mientras que sulpirida (N05AA/N05AX) presentó un SR de 1,95 (IC 95%: 1,90–2,00). El análisis anual mostró una reducción general del SR BZD-AV hasta 0,90 en 2023 (Figure 6).

In [14]:

knitr::include_graphics("figures/sr_anual_ambas_cascadas.png")

Figure 6: Ratio de secuencia crudo por año calendario del índice para ambas cascadas. El tamaño del punto es proporcional al número de pacientes analizables. La línea discontinua indica el SR global. NESR = 1 (sin corrección por tendencia secular).

Discusión

Este estudio encontró una asimetría cruda modesta en la secuencia de prescripción de diurético de asa tras gabapentinoides (SR 1,20) y, en menor medida, de antivertiginoso tras benzodiacepinas (SR 1,07), compatible con la existencia de cascadas terapéuticas descritas previamente para ambas combinaciones. Sin embargo, la comparación principal del estudio —si esa asimetría es mayor en personas de 65 o más años con polifarmacia que en el resto de la población BIFAP— no se sostuvo tras el ajuste por sexo, edad, carga de comorbilidad y confusores por indicación: los OR ajustados fueron próximos a 1 y no estadísticamente significativos en ninguna de las dos cascadas (GABA-DA: OR 0,85; BZD-AV: OR 1,02). El RSR crudo, calculado sin ajuste, tampoco mostró una diferencia relevante en GABA-DA (1,00) y solo una diferencia pequeña, aunque estadísticamente significativa por el gran tamaño muestral, en BZD-AV (1,04).

Comparación con la literatura

La cascada gabapentinoide–diurético de asa se ha descrito previamente en Canadá [16], Estados Unidos y Dinamarca [15], y mediante el mismo diseño ASSP aplicado conjuntamente a ambas cascadas en una base de reclamaciones japonesa [18], estudio cuyo planteamiento metodológico sirvió de referencia directa para este trabajo; la de benzodiacepina–antivertiginoso comparte el mismo mecanismo propuesto de tratar un efecto adverso (edema, inestabilidad) con un fármaco adicional en lugar de reconsiderar el fármaco causal [5,6]. Nuestros resultados son consistentes con esos estudios en cuanto a la dirección del SR crudo (>1 en ambas cascadas), pero no reproducen una asociación diferencial clara en el subgrupo de mayor edad y polifarmacia, que era la hipótesis específica de este trabajo para la población española de BIFAP. Un estudio de predictores en veteranos estadounidenses encontró que, a diferencia de la edad avanzada (asociada a menor riesgo), la polifarmacia extensa (≥10 fármacos) sí se asoció con mayor riesgo de la cascada GABA-DA [28], un patrón que no logramos replicar con nuestra definición combinada de edad y polifarmacia (≥5 fármacos crónicos).

Una explicación posible es que el fenómeno de cascada, si existe, no sea específico de los pacientes mayores polimedicados sino que afecte de forma más homogénea a todos los biexpuestos, lo que diluiría la señal en una comparación de subgrupos en lugar de eliminarla. Otra explicación, no excluyente, es que buena parte de la asimetría cruda observada refleje confusión por indicación o canalización clínica (channeling) —pacientes con más contactos con el sistema sanitario en general reciben más prescripciones de todo tipo en ventanas próximas al inicio de un tratamiento crónico— más que una cascada terapéutica propiamente dicha; esta interpretación es compatible con la sensibilidad del SR a la definición de ventana temporal observada en los análisis de sensibilidad, y con el hecho de que confusores por indicación inversa (insuficiencia cardíaca e insuficiencia renal en GABA-DA; síndrome vertiginoso crónico en BZD-AV) mostraran los coeficientes más marcados del modelo. Esta cautela es coherente con hallazgos recientes en atención primaria neerlandesa, donde solo el 28,6% de las secuencias de prescripción identificadas mediante una lista de cascadas de consenso internacional (ThinkCascades) se confirmaron como cascadas verdaderas tras revisión manual de la historia clínica [29], lo que subraya que un SR crudo elevado no equivale por sí solo a evidencia clínica de causalidad. Un estudio posterior en residencias de mayores estadounidenses estimó una incidencia de la cascada GABA-DA del 4,2% e identificó predictores clínicos individuales de su aparición [30], una aproximación centrada en la predicción a nivel de paciente que complementa la comparación de subgrupos por edad y polifarmacia realizada en este trabajo.

Implicaciones clínicas

Los hallazgos no respaldan que la vigilancia de estas cascadas deba concentrarse selectivamente en pacientes mayores con polifarmacia frente al resto de la población tratada con gabapentinoides o benzodiacepinas: la asimetría temporal cruda, cuando existe, parece un fenómeno más general. Esto no resta relevancia clínica a la revisión sistemática de nuevas prescripciones de diuréticos de asa tras gabapentinoides o de antivertiginosos tras benzodiacepinas como oportunidad de deprescripción [7,8,21], especialmente porque los gabapentinoides y las benzodiacepinas son de uso extendido en atención primaria en España [4] y ambas combinaciones son plausibles y ya descritas en otros sistemas sanitarios. Esta recomendación se refuerza por una revisión de historias clínicas en veteranos estadounidenses que mostró que los clínicos rara vez consideraron el gabapentinoide como causa del edema al prescribir el diurético, y que la cascada se asoció con daños posteriores como deterioro de la función renal, alteraciones electrolíticas y caídas [31].

Limitaciones

La principal limitación es la ausencia de tasas poblacionales de prescripción del fármaco evento, necesarias para estimar el NESR. Esta ausencia se asume como una limitación permanente del estudio: dichas tasas no forman parte de los extractos de BIFAP disponibles para este proyecto y no se espera obtenerlas, por lo que ninguna versión futura de este análisis, con los datos actuales, podrá corregir el SR por tendencia secular. En consecuencia, los SR presentados no están corregidos por tendencia secular y no debe asumirse que NESR = 1 sea necesariamente conservador; los resultados deben interpretarse como una medida cruda de asimetría temporal, no como una estimación ajustada del riesgo de cascada terapéutica.

La amplitud de la ventana, el tratamiento de las coprescripciones del día índice y la heterogeneidad entre grupos ATC modificaron sustancialmente las estimaciones. Esta dependencia de las decisiones analíticas limita la interpretación causal y la generalización del SR agregado, especialmente en BZD-AV. La estratificación por año calendario no sustituye al NESR poblacional.

Otra limitación podría encontrarse debido al amplio abanico de indicaciones que estos fármacos podrían tener en este grupo poblacional. Si de toda la proporción de eventos asociables al fármaco esta prescripción en cascada supusiese un bajo número de ellos, es posible que la asociación final esté atenuada [24].

Se ha demostrado previamente que el ASSP tiene una sensibilidad del 61%, una especificidad del 93%, un valor predictivo positivo del 77% y un valor predictivo negativo del 87%, utilizando eventos adversos identificados en ensayos controlados aleatorios como el gold estándar para eventos positivos y controles negativos de eventos adversos identificados en la información del producto de fármacos no relacionados a priori [25,32].

Conclusión

En la población de BIFAP se observó una asimetría temporal cruda modesta, compatible con las cascadas GABA-DA y BZD-AV descritas en otros países, pero no una asociación independiente entre pertenecer al subgrupo de 65 o más años con polifarmacia y el orden temporal de la prescripción, una vez ajustado por comorbilidad y confusores por indicación. La imposibilidad permanente de corregir por tendencia secular (NESR), la sensibilidad del SR a la definición de ventana y del día índice, y la heterogeneidad farmacológica dentro de cada cascada limitan la interpretación causal de estos resultados y aconsejan tratarlos como una estimación exploratoria de la magnitud del fenómeno en atención primaria española, más que como una cuantificación definitiva del riesgo de cascada terapéutica.

Referencias

1. Jyrkkä J, Enlund H, Korhonen MJ, Sulkava R, Hartikainen S. Polypharmacy status as an indicator of mortality in an elderly population. Drugs & Aging. 2009;26:1039–48.

2. Calderón-Larrañaga A et al. Multimorbidity, polypharmacy, referrals, and adverse drug events: Are we doing things well? British Journal of General Practice. 2012;62:e821–6.

3. Prados-Torres A et al. Cohort profile: The epidemiology of chronic diseases and multimorbidity. The EpiChron cohort study. International Journal of Epidemiology. 2018;47:382–384f.

4. Hernández-Rodríguez MÁ et al. Evolution of polypharmacy in a spanish population (2005-2015): A database study. Pharmacoepidemiology and Drug Safety. 2020;29:433–43.

5. Rochon PA, Gurwitz JH. Optimising drug treatment for elderly people: The prescribing cascade. BMJ. 1997;315:1096–9.

6. Rochon PA, Gurwitz JH. The prescribing cascade revisited. The Lancet. 2017;389:1778–80.

7. McCarthy LM, Visentin JD, Rochon PA. Assessing the scope and appropriateness of prescribing cascades. Journal of the American Geriatrics Society. 2019;67:1023–6.

8. Brath H et al. What is known about preventing, detecting, and reversing prescribing cascades: A scoping review. Journal of the American Geriatrics Society. 2018;66:2079–85.

9. Kalisch LM, Caughey GE, Roughead EE, Gilbert AL. The prescribing cascade. Australian Prescriber. 2011;34.

10. Carollo M, Crisafulli S, Zerio M, Maccarrone F, Paladin A, Selleri M, et al. Prescribing cascades: An umbrella review and updated systematic review. Drugs & Aging. 2026;43:397–424. https://doi.org/10.1007/s40266-026-01295-9

11. Rogero-Blanco E et al. Drug interactions detected by a computer-assisted prescription system in primary care patients in spain: MULTIPAP study. European Journal of General Practice. 2021;27:90–6.

12. Rogero-Blanco E et al. Use of an electronic clinical decision support system in primary care to assess inappropriate polypharmacy in young seniors with multimorbidity: Observational, descriptive, cross-sectional study. JMIR Medical Informatics. 2020;8:e14130.

13. Lopez-Rodriguez JA et al. Potentially inappropriate prescriptions according to explicit and implicit criteria in patients with multimorbidity and polypharmacy. MULTIPAP: A cross-sectional study. PLoS ONE. 2020;15:e0237186.

14. Fahrni ML, Azmy MT, Usir E, Aziz NA, Hassan Y. Inappropriate prescribing defined by STOPP and START criteria and its association with adverse drug events among hospitalized older patients: A multicentre, prospective study. PLoS ONE. 2019;14:e0219898.

15. Vouri SM et al. Association between gabapentinoids and oedema treated with loop-diuretics: A pooled sequence symmetry analysis from the united states and denmark. British Journal of Clinical Pharmacology. 2022; https://doi.org/10.1111/bcp.15447

16. Read SH et al. Evidence of a gabapentinoid and diuretic prescribing cascade among older adults with lower back pain. Journal of the American Geriatrics Society. 2021;69:2842–50.

17. Chen Z, Liu Z, Zeng L, Huang L, Zhang L. Research on prescribing cascades: A scoping review. Frontiers in Pharmacology. 2023;14:1147921.

18. Omata R, Asami A, Hara A, Urushihara H. Prescribing cascades of loop diuretics and anti-vertigo drugs following treatment with gabapentinoids and benzodiazepines: Prescription sequence symmetry analysis of a large-scale claims database including japanese older adults. Drugs - Real World Outcomes. 2024;11:529–40. https://doi.org/10.1007/s40801-024-00446-x

19. Agencia Española de Medicamentos y Productos Sanitarios. Observatorio de uso de medicamentos [Internet]. 2019. https://www.aemps.gob.es/medicamentos-de-uso-humano/observatorio-de-uso-de-medicamentos/

20. Ventura P. Utilización de medicamentos antiepilépticos en españa durante el periodo 2008-2016. 2016.

21. Lavan AH, Gallagher PF, O’Mahony D. Methods to reduce prescribing errors in elderly patients with multimorbidity. Clinical Interventions in Aging. 2016;11:857–66.

22. Maciá-Martínez M-Á et al. Base de datos para la investigación farmacoepidemiológica en atención primaria (BIFAP): A data resource for pharmacoepidemiology in spain. Pharmacoepidemiology and Drug Safety. 2020;29:1236–45.

23. Hallas J. Evidence of depression provoked by cardiovascular medication: A prescription sequence symmetry analysis. Epidemiology. 1996;7:478.

24. Lai EC-C et al. Sequence symmetry analysis in pharmacovigilance and pharmacoepidemiologic studies. European Journal of Epidemiology. 2017;32:567–82.

25. Morris EJ et al. Evaluating the use of prescription sequence symmetry analysis as a pharmacovigilance tool: A scoping review. Research in Social and Administrative Pharmacy. 2022;18:3079–93.

26. BIFAP: Base de datos para la investigación farmacoepidemiológica en atención primaria [Internet]. http://www.bifap.org/index.html\#content4-3u

27. Petri H et al. Prescription sequence analysis: A new and fast method for assessing certain adverse reactions of prescription drugs in large populations. Statistics in Medicine. 1988;7:1171–5.

28. Growdon ME, Jing B, Morris EJ, Deardorff WJ, Boscardin WJ, Byers AL, et al. Which older adults are at highest risk of prescribing cascades? A national study of the gabapentinoid-loop diuretic cascade. Journal of the American Geriatrics Society. 2024;72:1728–40. https://doi.org/10.1111/jgs.18892

29. Doherty AS, Perry M, Moriarty F, et al. Characterising prescribing cascades in older community-dwelling adults attending general practice in the netherlands: A retrospective cohort study. BMC Primary Care. 2026;27:156. https://doi.org/10.1186/s12875-026-03273-x

30. Hayes KN, Belanger E, Oganisian A, et al. Predictors of a gabapentinoid-loop-diuretic prescribing cascade in u.s. Nursing home residents. Journal of the American Geriatrics Society. 2026;74:336–44. https://doi.org/10.1111/jgs.70219

31. Growdon ME, Tjota N, Campbell R, Gayda P, Jing B, Deardorff WJ, et al. Decision-making and downstream outcomes of the gabapentinoid-diuretic prescribing cascade. JAMA Network Open. 2025;8:e2545274. https://doi.org/10.1001/jamanetworkopen.2025.45274

32. Wahab IA, Pratt NL, Wiese MD, Kalisch LM, Roughead EE. The validity of sequence symmetry analysis (SSA) for adverse drug reaction signal detection. Pharmacoepidemiology and Drug Safety. 2013;22:496–502.

Supplementary files

Supplementary Table 1. Fármacos de exposición incluidos en el análisis.

In [15]:

lookup_exp <- readr::read_csv(
  file.path(project_root, "data", "lookup", "anexo_I_exposicion_gaba_bzd.csv"),
  show_col_types = FALSE
) |>
  dplyr::bind_rows(tibble::tibble(
    subgrupo_atc = "N05BA", descripcion = "Derivado de la benzodiacepina",
    principio_activo = "Diazepam"
  )) |>
  dplyr::distinct(subgrupo_atc, principio_activo, .keep_all = TRUE) |>
  dplyr::mutate(
    Cascada = dplyr::if_else(
      principio_activo %in% c("Gabapentina", "Pregabalina"), "GABA-DA", "BZD-AV"
    )
  ) |>
  dplyr::arrange(Cascada, subgrupo_atc, principio_activo) |>
  dplyr::select(Cascada, `Código ATC` = subgrupo_atc,
                `Principio activo` = principio_activo, Descripción = descripcion)

supp_exp_gt <- lookup_exp |>
  gt::gt(groupname_col = "Cascada") |>
  gt::cols_label(
    `Código ATC` = gt::md("**Código ATC**"),
    `Principio activo` = gt::md("**Principio activo**"),
    Descripción = gt::md("**Descripción**")
  ) |>
  gt::tab_options(table.width = gt::pct(100))

knitr::asis_output(gt::as_raw_html(supp_exp_gt))

Código ATC	Principio activo	Descripción
BZD-AV
N03AE	Clonazepam	Derivado de la benzodiazepina
N05BA	Alprazolam	Ansioliticos
N05BA	Diazepam	Derivado de la benzodiacepina
N05BX	Bentazepam	Ansioliticos
N05BX	Bromazepam	Ansioliticos
N05BX	Buspirona	Ansioliticos
N05BX	Clobazam	Ansioliticos
N05BX	Clorazepato dipotasico	Ansioliticos
N05BX	Clordiazepoxido	Ansioliticos
N05BX	Clotiazepam	Ansioliticos
N05BX	Halazepam	Ansioliticos
N05BX	Hidroxizina	Ansioliticos
N05BX	Ketazolam	Ansioliticos
N05BX	Lorazepam	Ansioliticos
N05BX	Oxazepam	Ansioliticos
N05BX	Pinazepam	Ansioliticos
N05CA	Brotizolam	Hipnoticos y sedantes
N05CX	Clometiazol	Hipnoticos y sedantes
N05CX	Flunitrazepam	Hipnoticos y sedantes
N05CX	Flurazepam	Hipnoticos y sedantes
N05CX	Loprazolam	Hipnoticos y sedantes
N05CX	Lormetazepam	Hipnoticos y sedantes
N05CX	Midazolam	Hipnoticos y sedantes
N05CX	Quazepam	Hipnoticos y sedantes
N05CX	Triazolam	Hipnoticos y sedantes
N05CX	Zolpidem	Hipnoticos y sedantes
N05CX	Zopiclona	Hipnoticos y sedantes
GABA-DA
N03AX	Gabapentina	Otros antiepilepticos
N03AX	Pregabalina	Otros antiepilepticos

Supplementary Table 2. Fármacos evento incluidos en cada cascada.

In [16]:

lookup_evt <- readr::read_csv(
  file.path(project_root, "data", "lookup",
            "anexo_II_eventos_diureticos_antivertiginosos.csv"),
  show_col_types = FALSE
) |>
  dplyr::filter(
    subgrupo_atc == "C03CA" |
      subgrupo_atc %in% c("N05AA", "N05AX", "N07CA")
  ) |>
  dplyr::mutate(
    Cascada = dplyr::if_else(subgrupo_atc == "C03CA", "GABA-DA", "BZD-AV")
  ) |>
  dplyr::arrange(Cascada, subgrupo_atc, principio_activo) |>
  dplyr::select(Cascada, `Código ATC` = subgrupo_atc,
                `Principio activo` = principio_activo, Descripción = descripcion)

supp_evt_gt <- lookup_evt |>
  gt::gt(groupname_col = "Cascada") |>
  gt::cols_label(
    `Código ATC` = gt::md("**Código ATC**"),
    `Principio activo` = gt::md("**Principio activo**"),
    Descripción = gt::md("**Descripción**")
  ) |>
  gt::tab_options(table.width = gt::pct(100))

knitr::asis_output(gt::as_raw_html(supp_evt_gt))

Código ATC	Principio activo	Descripción
BZD-AV
N05AA	Sulpirida	Preparados contra el vertigo
N05AX	Sulpirida	Preparados contra el vertigo
N07CA	Betahistina	Preparados contra el vertigo
N07CA	Cinarizina	Preparados contra el vertigo
N07CA	Flunarizina	Preparados contra el vertigo
GABA-DA
C03CA	Bumetanida	Sulfonamida monofarmaco
C03CA	Furosemida	Sulfonamida monofarmaco
C03CA	Piretanida	Sulfonamida monofarmaco
C03CA	Torasemida	Sulfonamida monofarmaco

Supplementary Table 3. Fármacos confusores preespecificados incluidos en el ajuste.

In [17]:

lookup_conf <- readr::read_csv(
  file.path(project_root, "data", "lookup", "anexo_III_confusores.csv"),
  show_col_types = FALSE
) |>
  dplyr::arrange(subgrupo_atc, principio_activo) |>
  dplyr::select(`Código ATC` = subgrupo_atc,
                `Principio activo` = principio_activo)

supp_conf_gt <- lookup_conf |>
  gt::gt() |>
  gt::cols_label(
    `Código ATC` = gt::md("**Código ATC**"),
    `Principio activo` = gt::md("**Principio activo**")
  ) |>
  gt::tab_options(table.width = gt::pct(75))

knitr::asis_output(gt::as_raw_html(supp_conf_gt))

Código ATC	Principio activo
C02DB	Hidralazina
C08CA	Amlodipino
C08GA	Combinaciones
D11AX01	Minoxidil
M01AH	Celecoxib Etoricoxib
N04BC	Pramiprexol
N06AB10	Escitalopram
N06AX11	Mirtazapina

Supplementary Table 4. Validación post-hoc de la especificación del modelo de regresión logística multinivel.

In [18]:

library(gt)
library(readr)

tabla_val <- read_csv(
  file.path(outputs_dir, "tables", "tabla_validacion_modelos.csv"),
  show_col_types = FALSE
)

tbl_val_gt <- tabla_val |>
  dplyr::mutate(
    lrt = dplyr::if_else(
      !is.na(lrt_ccaa_chisq),
      sprintf(
        "\u03c7\u00b2=%.2f, p%s", lrt_ccaa_chisq,
        dplyr::if_else(lrt_ccaa_p < 0.001, "<0,001", sprintf("=%.3f", lrt_ccaa_p))
      ),
      "\u2014"
    ),
    auc_r2 = dplyr::if_else(
      !is.na(auc),
      sprintf("%.3f / %.3f", auc, r2_nagelkerke),
      "\u2014"
    )
  ) |>
  dplyr::select(cascada, modelo, n, aic, bic, lrt, auc_r2) |>
  gt(groupname_col = "cascada") |>
  cols_label(
    modelo = md("**Modelo**"),
    n      = md("**n**"),
    aic    = md("**AIC**"),
    bic    = md("**BIC**"),
    lrt    = md("**LRT CCAA**"),
    auc_r2 = md("**AUC / R\u00b2 Nagelkerke**")
  ) |>
  fmt_number(columns = c(aic, bic), decimals = 1, sep_mark = ".", dec_mark = ",") |>
  fmt_number(columns = n, decimals = 0, sep_mark = ".") |>
  tab_footnote(
    "M3: modelo completo (enfermedades espec\u00edficas + efecto aleatorio CCAA); M2: sin enfermedades espec\u00edficas; M1: sin efecto aleatorio CCAA (glm). AUC y R\u00b2 de Nagelkerke calculados sobre M3.",
    locations = cells_column_labels(columns = modelo)
  ) |>
  tab_options(table.width = pct(100),
              table.border.top.color = "#2b3a55",
              column_labels.border.bottom.color = "#2b3a55",
              table_body.hlines.color = "#dee2e6")

knitr::asis_output(as_raw_html(tbl_val_gt))

In [19]:

Table 4: Validación post-hoc de la especificación del modelo: comparación de AIC/BIC entre el modelo completo (M3), sin enfermedades específicas de indicación (M2) y sin efecto aleatorio por CCAA (M1); prueba de razón de verosimilitudes (LRT) para el efecto aleatorio CCAA; discriminación (AUC) y pseudo-R² de Nagelkerke del modelo M3.

Modelo1	n	AIC	BIC	LRT CCAA	AUC / R² Nagelkerke
GABA-DA
M3: completo (enf. + CCAA)	11.451	15.553,6	15.663,8	χ²=1.95, p=0.081	0.579 / 0.029
M2: sin enfermedades específicas (+ CCAA)	11.451	15.723,4	15.796,9	—	—
M1: sin CCAA (glm)	11.451	15.553,6	15.656,4	—	—
BZD-AV
M3: completo (enf. + CCAA)	102.295	140.041,5	140.194,1	χ²=66.35, p<0,001	0.563 / 0.021
M2: sin enfermedades específicas (+ CCAA)	102.295	141.227,5	141.322,8	—	—
M1: sin CCAA (glm)	102.295	140.105,9	140.248,9	—	—
1 M3: modelo completo (enfermedades específicas + efecto aleatorio CCAA); M2: sin enfermedades específicas; M1: sin efecto aleatorio CCAA (glm). AUC y R² de Nagelkerke calculados sobre M3.

La especificación de las covariables fue a priori (razonamiento clínico y protocolo), no seleccionada por criterio estadístico; la comparación de modelos en Table 4 es una validación posterior, no una selección de modelo. Incluir las enfermedades específicas de indicación mejoró el ajuste en ambas cascadas. El efecto aleatorio por CCAA estuvo claramente justificado en BZD-AV (LRT p < 0,001) pero fue marginal en GABA-DA (LRT p = 0,081; el BIC favoreció el modelo sin CCAA); se mantuvo en ambas cascadas por coherencia metodológica y por ser el efecto aleatorio pre-especificado en el protocolo. La discriminación fue baja en ambos modelos (AUC 0,58 y 0,56; pseudo-R² de Nagelkerke 0,03 y 0,02, respectivamente), lo cual es esperable dado que el modelo describe el orden temporal entre biexpuestos (un diseño ASSP) y no pretende predecir la prescripción del fármaco evento.

Supplementary Picture 1. Representación esquemática de la entrada en la cohorte y de las ventanas temporales del análisis de secuencias simétricas de prescripción.