94  Incorporare dati storici di controllo in una RCT

Prerequisiti

Concetti e competenze chiave

Preparazione del Notebook

import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import scipy.stats as stats
import arviz as az
import warnings
from cmdstanpy import cmdstan_path, CmdStanModel

seed: int = sum(map(ord, "stan_rct"))
rng: np.random.Generator = np.random.default_rng(seed=seed)
az.style.use("arviz-darkgrid")
%config InlineBackend.figure_format = "retina"

# Get the home directory
home_directory = os.path.expanduser("~")
# Construct the path to the Quarto project directory
project_directory = os.path.join(home_directory, "_repositories", "psicometria")

Introduzione

Questo capitolo fornisce una trattazione semplificata di un importante problema affrontato da Frank Harrell in un suo intervento intitolato Incorporating Historical Control Data Into an RCT.

94.1 Studi Controllati Randomizzati

Nella ricerca psicologica, ci troviamo di fronte a notevoli ostacoli nel reclutare un numero sufficiente di partecipanti per condurre studi controllati randomizzati (RCT), un problema che si acuisce quando si prevede l’assegnazione dei partecipanti a gruppi di controllo che ricevono trattamenti standard. Un’altra difficoltà sorge quando i potenziali partecipanti sono riluttanti a iscriversi agli studi a causa della possibilità di non ricevere il trattamento sperimentale. In questo contesto, l’utilizzo di Dati Storici (HD) per informare su possibili esiti nei gruppi di controllo assume un’importanza vitale. Tuttavia, l’integrazione di tali dati richiede strategie sofisticate per adeguare i bias e le disparità tra i diversi disegni di studio.

Stuart Pocock ha proposto un metodo nel quadro frequentista che valorizza la dimensione campionaria degli HD, pur ammettendo che questi possano riflettere realtà diverse rispetto agli esiti attesi nei gruppi di controllo degli RCT prospettici. La discrepanza include sia la vera performance sconosciuta del gruppo di controllo sia il bias inerente agli HD.

Björn Holzhauer ha ampliato questa visione attraverso lo sviluppo di approcci Bayesiani per l’appropriazione di dati, in particolare riguardo ai tassi di pericolo esponenziali, mediante l’utilizzo di simulazioni MCMC Bayesiane. Queste metodologie consentono l’elaborazione parallela di più modelli e l’inclusione diretta dei dati grezzi degli HD nell’analisi di nuovi dati sperimentali, affrontando direttamente le possibili discrepanze negli oggetti di stima tra HD e RCT.

Una caratteristica fondamentale di queste tecniche è l’aggregazione di dati grezzi da diverse fonti, facilitando un’analisi più precisa che considera l’intero spettro delle incertezze. Questo si contrappone agli approcci tradizionali, che spesso si basano su statistiche riassuntive e tendono a trascurare importanti variabilità, conferendo una fiducia ingiustificata nei dati storici. È, inoltre, cruciale l’ajustamento per covariate al fine di gestire l’eterogeneità degli esiti tra i trattamenti, aumentando così l’affidabilità e l’accuratezza delle stime dell’effetto del trattamento.

Procedendo senza dati diretti per valutare il bias, ci affidiamo a una distribuzione a priori gaussiana con media zero e deviazione standard sigma per descriverlo. Un valore di sigma nullo implica l’assenza di bias, permettendo di integrare gli HD nell’analisi allo stesso livello dei dati di controllo dello studio. Al contrario, un sigma infinito indica una completa ignoranza riguardo al bias, rendendo gli HD non informativi e pertanto trascurabili.

Il focus principale, l’effetto del trattamento delta, viene esaminato più efficacemente attraverso l’analisi di sigma. Un sigma inferiore a 2 rende lo studio informativo su delta, mentre un sigma superiore a 2 può portare a conclusioni errate, specialmente quando gli HD hanno una media artificiosamente gonfiata. Questo sottolinea l’importanza di una scelta accurata dei parametri analitici, in particolare nell’integrazione dei dati storici con quelli dei nuovi studi clinici, per fornire una rappresentazione equilibrata e critica dell’integrazione di tali dati nell’analisi statistica contemporanea.

94.2 Creazione di un Braccio di Controllo con HD

Nel contesto di uno studio RCT condotto con un unico braccio sperimentale, dove i dati di controllo derivano unicamente da controlli storici non contemporanei, la sfida si amplifica. In queste circostanze, il bias intrinseco negli HD non può essere quantificato direttamente. Invece, ci si affida completamente alla distribuzione di incertezza prescelta per il bias. Questo approccio offre un’analisi che, pur essendo paragonabile a una verifica di sensibilità, si avvale del rigoroso quadro analitico Bayesiano. Quest’ultimo, per sua natura flessibile, permette un’estensione naturale a include la meta-analisi di dati individuali dei pazienti e l’adeguamento per covariate, arricchendo così la robustezza e l’affidabilità delle inferenze tratte dallo studio.

Importiamo e compiliamo il modello.

stan_file = os.path.join(project_directory, "stan", "rct.stan")
model = CmdStanModel(stan_file=stan_file)
print(model.code())

12:00:16 - cmdstanpy - INFO - compiling stan file /Users/corradocaudek/_repositories/psicometria/stan/rct.stan to exe file /Users/corradocaudek/_repositories/psicometria/stan/rct
12:00:26 - cmdstanpy - INFO - compiled model executable: /Users/corradocaudek/_repositories/psicometria/stan/rct

data {
  int<lower=0> Nb;  // # obs in RCT treatment B
  int<lower=0> Nh;  // # obs in historical control data
  vector[Nb] yb;    // vector of tx=B data
  vector[Nh] yh;    // vector of historical data
  real<lower=0> sigma;  // standard deviation of prior for bias
}
parameters {
  real mua;  // unknown mean for tx=A
  real mub;  // unknown mean for tx=B
  real bias; // unknown bias
}
transformed parameters {
  real delta;
  delta = mua - mub;
}
model {
  yb   ~ normal(mub, 1.0);
  yh   ~ normal(mua + bias, 1.0);
  bias ~ normal(0., sigma);
}

Creiamo il dizionario che contiene i dati richiesti dal modello.

Na = 20
Nb = 40
Nh = 500
ya = rng.normal(loc=10, scale=1, size=Na)
yb = rng.normal(loc=5, scale=1, size=Nb)
yh = rng.normal(loc=20, scale=1, size=Nh)

stan_data = {"Nb": Nb, "Nh": Nh, "yb": yb, "yh": yh, "sigma": 1.0}

Eseguiamo il campionamento MCMC.

fit = model.sample(data=stan_data)

Esaminiamo i risultati ottenuti.

print(fit.summary())

             Mean      MCSE    StdDev         5%         50%        95%  \
lp__  -252.847000  0.035606  1.225860 -255.33900 -252.524000 -251.48100   
mua     20.034200  0.028866  1.027280   18.35000   20.040900   21.69820   
mub      5.216960  0.004029  0.160862    4.94483    5.216310    5.47481   
bias    -0.044748  0.028869  1.026390   -1.71126   -0.042136    1.62836   
delta   14.817300  0.029354  1.038220   13.10640   14.832900   16.49740   

         N_Eff  N_Eff/s    R_hat  
lp__   1185.29  1727.83  1.00003  
mua    1266.49  1846.20  1.00055  
mub    1594.19  2323.89  1.00245  
bias   1264.05  1842.64  1.00058  
delta  1250.95  1823.54  1.00088  

az.summary(fit, var_names=(["bias", "delta"]), hdi_prob=0.95)

	mean	sd	hdi_2.5%	hdi_97.5%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
bias	-0.045	1.026	-1.990	1.917	0.029	0.021	1265.0	1491.0	1.0
delta	14.817	1.038	12.911	16.837	0.029	0.021	1251.0	1285.0	1.0

_ = az.plot_trace(fit, var_names=(["bias", "delta"]))

Si noti l’utilità dei dati storici nella riduzione dell’incertezza associata a delta. Quando diminuisce sigma, diminuisce anche l’incertezza associata all’effetto del trattamento.

stan_data = {
    "Nb" : Nb, 
    "Nh" : Nh, 
    "yb" : yb, 
    "yh" : yh,
    "sigma" : 0.25
}

fit1 = model.sample(data=stan_data)

az.summary(fit1, var_names=(["bias", "delta"]), hdi_prob=0.95)

	mean	sd	hdi_2.5%	hdi_97.5%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
bias	-0.006	0.250	-0.505	0.468	0.007	0.005	1337.0	1404.0	1.0
delta	14.782	0.304	14.161	15.335	0.008	0.006	1483.0	1538.0	1.0

94.3 Conclusioni

Per stabilire la distribuzione a priori del bias, è fondamentale considerare le informazioni disponibili riguardanti l’evoluzione delle pratiche terapeutiche, il bias di selezione dei pazienti e l’andamento della patologia di interesse. L’approccio Bayesiano ci libera dalla necessità di conoscere con precisione l’entità del bias, indirizzandoci invece a definire la sua distribuzione di incertezza. In caso questa distribuzione sia modellata come gaussiana, ci si avvale dell’ipotesi di simmetria per focalizzarsi sulla deviazione standard, sigma. Un metodo efficace per determinare sigma consiste nel fissarlo in modo che, per esempio, la probabilità che il bias si collochi entro un intervallo di \([-c, c]\) sia del 95%, per poi calcolare retrospettivamente il valore di sigma necessario a soddisfare questa condizione. Questo approccio garantisce una gestione più mirata e scientificamente fondata dell’incertezza legata al bias, fondamentale per l’integrazione ottimale dei dati storici nelle analisi statistiche avanzate.

Informazioni sull’Ambiente di Sviluppo

%load_ext watermark
%watermark -n -u -v -iv -w -m -p cmdstanpy

Last updated: Tue Jul 30 2024

Python implementation: CPython
Python version       : 3.12.4
IPython version      : 8.26.0

cmdstanpy: 1.2.4

Compiler    : Clang 16.0.6 
OS          : Darwin
Release     : 23.6.0
Machine     : arm64
Processor   : arm
CPU cores   : 8
Architecture: 64bit

arviz     : 0.18.0
pandas    : 2.2.2
scipy     : 1.14.0
numpy     : 1.26.4
seaborn   : 0.13.2
matplotlib: 3.9.1

Watermark: 2.4.3