Modelli di rischio e analisi di dati epidemiologici

Premessa

Questa nota è il primo approfondimento che fa seguito alla lettera sottoscritta da me e da altri comparsa su Epidemiologia & Prevenzione¹ a proposito del Position Paper della Società Italiana di Medicina del Lavoro, pubblicato ora anche su La Medicina del Lavoro.² Facevamo allora notare alcune mancanze e alcune formulazioni contenute nel documento che, a nostro avviso, meritavano qualche integrazione, chiarimento e approfondimento, al fine di non prestarsi a travisamenti e strumentalizzazioni, soprattutto in sede giudiziaria, a danno di lavoratori esposti ad amianto e danneggiati da tale esposizione. In particolare facevamo riferimento ai temi della relazione dose-risposta, durata dell’induzione e periodo di latenza; a quello della diagnosi clinica del mesotelioma maligno; a quello dell’uso e dell’interpretazione dei modell statistici. Segnalavamo allora come il documento avesse ignorato del tutto i risultati degli studi sperimentali, nel trattare la relazione dose-risposta e la funzione di rischio; non avesse considerato l’approccio clinico pratico nella diagnosi di mesotelioma maligno (vista anche l’assenza dal gruppo di  lavoro di competenze del clinico medico, del radiologo, del chirurgo, dell’oncologo); avesse trattato in modo poco approfondito l’uso dei modelli statistici nell’analisi dei dati epidemiologici, contribuendo a rafforzare l’idea che da essi non siano ricavabili risultati validi per il caso singolo, quando questa traslazione di risultati avviene regolarmente nel caso della Medicina basata sulle prove di efficacia (Evidence-Based Medicine, EBM), con la decisione per il singolo caso presa proprio in base a risultati ottenuti, come valori attesi, in sperimentazioni cliniche controllate di gruppo.
Osservazioni al documento sono state espresse anche, da altra prospettiva, da medici del lavoro (vd. questo fascicolo p. 73).³ Il nostro impegno deriva dal considerare responsabilità degli “esperti” contribuire a chiarire le posizioni condivise dalla comunità scientifica per quanto riguarda la relazione causale fra esposizione ad amianto e insorgenza dei tumori. Di seguito considererò il significato dell’uso di modelli di probabilità in epidemiologia, in particolare di quelli applicati all’analisi dei dati di studi longitudinali prospettici (di coorte), o analisi della sopravvivenza, svolta sia per stimare Hazard Ratio (HR) o Rischio Relativo (RR) dell’esposizione sia per  stimare la anticipazione temporale dell’esito o Rate Advancement Period (RAP) attribuibile all’esposizione. Mi soffermerò a dimostrare come i due tipi di risultato non  siano che formulazioni distinte, con diverse metriche, della stessa situazione di rischio per gli esposti. La prima più appropriata a esprimere la forza di causazione  dell’esposizione nel produrre l’insorgenza della patologia, la seconda più appropriata a esprimere l’entità del danno inferto dall’esposizione al soggetto che ha sviluppatola patologia fra gli esposti.

Introduzione

Ciò che conosciamo è sempre molto meno di quello che ignoriamo. Più aumenta il raggio della conoscenza e più si accresce la consapevolezza della profondità della nostra ignoranza. Ma allora in che cosa consiste il “progresso della conoscenza”? Mossi dalla curiosità di scoprire “la verità” e dalla necessità di soddisfare “i  bisogni”, constatiamo che spesso conoscenze fallaci e “parziali” si rivelano utili nella pratica.

1. Una teoria che considera la terra piatta, cui si applica la geometria euclidea, è di fatto bastata a progettare e costruire le piramidi, ai tempi dei faraoni, l’Empire State Building, in tempi più recenti. Da tempo sappiamo che la terra è sferica, all’incirca, ma ai costruttori non fu necessario tenere conto di tale “verità”.
2. La teoria della gravitazione universale è bastata a impostare i calcoli su cui basare l’esplorazione dello spazio, fino al successo della discesa sul suolo lunare. Se dal 1905 abbiamo imparato che né tempo né spazio sono assoluti, la “verità” che ci ha svelato Einstein non fu necessaria per tracciare le prime rotte spaziali.
3. La teoria della relatività è oggi indispensabile per garantire accuratezza e precisione al sistema satellitare di geolocalizzazione. Impossibile trascurare le diversità del fluire del tempo in diversi punti dello spazio, senza incorrere in gravi errori. John Snow, convinto assertore della responsabilità dell’acqua prelevata a valle degli scarichi fognari di Londra per le periodiche epidemie di colera, fondò sui risultati di ricerche epidemiologiche la richiesta di spostare la captazione dell’acquapotabile a monte degli scarichi fognari, e questo ben prima che Robert Koch identificasse nel vibrione del colera il fattore eziologico della malattia. Teorie “vere”, che meglio ci fanno comprendere la realtà, continueranno a sostituirsi nel tempo a quelle “false”, per essere sostituite a loro volta. Ma teorie che verranno superate risulteranno spesso utili per risolvere problemi. Questo vale anche per i modelli, a proposito dei quali George E.P. Box (1919-2013), statistico, pioniere del controllodi qualità, dell’analisi di serie storiche, del disegno degli esperimenti, dell’inferenza bayesiana, sosteneva:
   
   

   «Tutti i modelli sono sbagliati, alcuni sono utili».

   
   In senso generale si può quindi affermare che un modello, lungi dal voler rappresentare la “verità”, sotto condizioni da verificare (assunti) ed entro limiti da rispettare (campo di applicazione), consente previsioni utili relativamente al fenomeno cui bene si adatta. Esso, generalmente matematico, consiste di relazioni funzionali e di parametri, la cui validazione si ottiene dimostrando il buon adattamento dei valori attesi ai dati osservati. In questo consiste l’approccio ipotetico-deduttivo che Karl Popper ha indicato essere fondamento per la “scoperta scientifica”. Il modello deterministico, applicato ai problemi della fisica classica, fornisce una previsione esatta, nei limiti dell’incertezza data da errori casuali di misura. Il modello probabilistico, applicato in biologia, nella ricerca biomedica, nella fisica delle particelle, fornisce previsioni probabilistiche, basate su “valori attesi”, non meno utili e valide delle prime per una decisione a livello individuale, come dimostra la pratica ormai consolidata della Medicina delle Prove di Efficacia (EBM).

Studi epidemiologici di occorrenza: modelli PH e AFT

Per costruire un modello occorre procedere a:

1. analisi del problema, con definizione delle finalità che il ricercatore si propone di ottenere;
2. traduzione del problema nel sistema di relazioni logicofunzionali del modello matematico;
3. stima dei parametri;
4. validazione e verifica degli assunti e della congruenza fra i risultati attesi e i dati osservati.

Nel caso di studio di coorte un gruppo di persone è osservato per un congruo periodo di tempo allo scopo di rilevare il presentarsi di un esito di interesse, malattia o decesso, e il tempo al quale si realizza. Il risultato sarà che alcuni soggetti reclutati nello studio, entro il periodo programmato di osservazione, avranno sperimentato l’esito in studio, altri usciranno indenni dallo studio, o prima della sua conclusione o alla fine. In sintesi, di N reclutati, D avranno manifestato l’esito ai tempi rilevati, (N-D) avranno tempi di osservazione troncati. Per la validità dei risultati occorre che il troncamento sia indipendente dall’esito. Poiché si indica con il termine “failure” il verificarsi e con “survival” il non verificarsi dell’evento, analisi di questo tipo si chiamano “survival analyses”.
Per procedere sono stati definiti due modelli: quello di “Proportional Hazard” (PH) e quello di “Accelerated Failure Time” (AFT). Il modello PH analizza la relazione fra le covariate – le variabili che l’epidemiologo considera importanti per il fenomeno, e fra esse l’esposizione al fattore di rischio in studio – e la variabile risposta; quindi stima l’effetto moltiplicativo di ciascuna variabile sul rischio di base (h₀), quale esso sia. Il modello AFT, invece, analizza la relazione fra le covariate e la variabile  risposta come accelerazione del tempo di occorrenza, assumendo che questa, come generalmente capita in caso di patologia cronico-degenerativa, aumenti  comunque con l’età, in modo indipendente dall’esposizione a un particolare fattore di rischio. Da sottolineare che chi propose per primo il modello PH⁴ suggerì la possibilità di analizzare gli stessi dati in termini di AFT. Il modello PH esprime la forza di associazione fra evento e fattore di rischio ed è più utile alla interpretazione
in senso causale della relazione. Il modello AFT esprime l’impatto dell’esposizione al fattore di rischio sull’esperienza di vita dell’individuo esposto, che ha  sperimentato l’evento. I risultati delle due analisi, la stima di “hazard ratio (HR)”, o RR, nel caso PH, di “time ratio (TR)”, o RAP, nel caso AFT, sono modi per  esprimere l’effetto dell’esposizione secondo metriche differenti. Il modello PH è semi-parametrico e non necessita di modellare il rischio di base (h₀) ma soltanto  assume la costanza del rischio proporzionale:

HR = h₁(t)/h₀(t) = exp(β) = RR

Ma il rischio di base, h₀, può essere espresso in funzione del tempo e presentarsi come invariante (è il caso di failure di molti prodotti industriali), monotonicamente crescente o decrescente (è il caso di molte condizioni di interesse medico), variabile in modo non monotonico. In ogni caso è possibile calcolare la anticipazione temporale dell’evento dovuta alla esposizione calcolando il tempo al quale, nel non esposto, si ottiene la stessa stima di rischio osservata nell’esposto. Questo indipendentemente dalla forma funzionale di relazione fra h e t adattata ai dati osservati, quindi non solo in caso di distribuzione di Weibull.
La distribuzione di densità di probabilità dei tempi all’evento nota come distribuzione di Weibull è la più utilizzata in quanto, per la sua flessibilità, è adattabile a diverse osservazioni empiriche. Per esempio, quando rischio è costante, e cioè h₀(t) = h₀, la distribuzione dei tempi di occorrenza, esponenziale e dipendente solo dalla incidenza istantanea (λ), è rappresentata dal modello di Weibull con parametri γ = 1 e λ. Quando il rischio h₀(t) è monotonicamente crescente o decrescente, la distribuzione dei tempi all’evento è rappresentata dal modello di Weibull con parametro γ rispettivamente maggiore di 1 e minore di 1 (figura 1). Se la relazione è non monotonica si ricorre ad altri modelli di distribuzione dei tempi all’evento e in ogni caso, come detto, è possibile ricavare la stima di anticipazione temporale attribuibile alla esposizione.
Per dare un’idea della flessibilità del modello di Weibull nella descrizione della distribuzione di probabilità dei tempi all’occorrenza e della relazione fra rischio (h) e tempi (t), riporto in figura 2 i grafici relativi a modelli con tempo mediano di 20 anni, valori γ: 0,5; 1,5; 3,0 (λ: 0,15499; 0,00775; 0,0000866). Nella figura 1 riporto i grafici della relazione fra rischio e tempi, per Weibull con valori di γ: 0,8; 1,0; 2,0; 2,5 da cui risulta un andamento monotonicamente decrescente, costante, crescente in modo pressoché lineare e con modalità tendenzialmente esponenziale.

 

La possibilità di tradurre i risultati ottenuti con il modello PH (RR) in quelli del modello AFT (RAP) dipende dalle relazioni matematiche che legano i due modelli.⁵ Dalla formula del modello di Weibull:

f(t)=λγt^γ-1exp(-λt^γ),

dove:
f(t) = funzione di densità di probabilità dei tempi all’occorrenza (t),
λ è parametro di scala e γ parametro di forma della distribuzione (entrambi >0);

h(t) = f(t)/S(t),

dove:
h(t) = funzione di rischio variabile nel tempo (t),
S(t) = 1 – F(t) con F(t) = funzione cumulativa dei tempi all’occorrenza e quindi S(t) funzione di sopravvivenza, proporzione di soggetti a rischio di evento al tempo t

se γ=1:
f(t) = λexp[-λt], S(t) = 1 – F(t) = exp(-λt), h(t) = λ
e i percentili della distribuzione dei tempi saranno:
t(p) = 1/λ Ln(100/(100-p));
mediana per esempio: t(.5) = Ln(2)/λ.

se γ≠1:
f(t) = λγt^γ-1exp(-λt^γ), S(t) = exp[-λt^γ], h(t) = λγt^γ-1
e i percentili della distribuzione dei tempi saranno:
t(p) = [1/λ Ln(100/(100-p))]^1/γ;
mediana t(.5) = (Ln(2)/λ)^1/γ.

Così i dati di uno studio di esposizione a un fattore E possono essere analizzati in termini di accelerazione del tempo stimando a quale valore di età A₀ h₀(t), rischio nel non esposto, eguaglia h₁(t), rischio nell’esposto, all’età A₁:

h_i(t)= h₀ exp(β_EE +β_AA+ Æ©β_jX_j)

con: E = 1, esposti; E = 0, non esposti

dove:
h₁(t) = h₀ exp(β_E+β_AA+Æ©β_jX_j) con β_E, β_A, β_j coefficienti che
esprimono la forza di associazione fra rischio e esposizione,
età in anni, variabili di confondimento;
h₀(t) = h₀ exp(β_ΑA+Æ©β_jX_j)

In questa formulazione, la anticipazione, differenza di età per cui i due valori h(t) si equivalgono, RAP = (A₀ – A₁), risulta eguale al rapporto β_Ε/β_Α.

Modelli di cancerogenesi e andamento dell'incidenza con l'età

Nello studio della relazione fra occorrenza del cancro ed età è fondamentale il lavoro di Armitage e Doll, celebrato a 50 anni di distanza come una pietra miliare con la ripubblicazione integrale su International Journal of Epidemiology.⁶ In tale lavoro si partiva dai tentativi fatti per formulare un modello biologico che rendesse conto della variazione osservata del tasso di mortalità da tumore con l’età, con l’interpolazione ai dati osservati di un modello statistico. In particolare si consideravano le ipotesi formulate da Fisher e Hollomon nel 1951 e da Nordling nel 1953. In entrambi i casi si considerava l’intervallo di età 25-74 anni, trovando una relazione lineare fra il logaritmo dei tassi e il logaritmo dell’età, con coefficiente di regressione lineare stimato di 6.
Fisher e Hollomon interpretarono tale risultato come indicazione del fatto che 6 o 7 fosse il numero di cellule cancerose da raggiungere per innescare la crescita tumorale. In tal caso l’incidenza sarebbe dovuta anche risultare proporzionale alla 5° o 6° potenza della concentrazione del cancerogeno, mentre i dati sperimentali suggerivano una relazione lineare.
Nordling interpretò l’osservazione come dovuta al fatto che una singola cellula diventasse cancerosa alla fine di 7 mutazioni ordinate. Questo faceva anche prevedere la relazione lineare osservata fra logaritmi dell’incidenza e della concentrazione del cancerogeno. Tale modello faceva anche prevedere che un aumento del tasso di mutazione da esposizione in giovane età si sarebbe  manifestato con aumentata incidenza in età avanzata, mentre lo stesso incremento di mutazioni da esposizione in età più avanzata non avrebbe comportato lo stesso aumento di incidenza. Tale fatto previsto dal modello e puntualmente osservato fornì una conferma indipendente alla validità del modello biologico di cancerogenesi multistadio.
In anni più recenti Berry ha affrontato il problema di stabilire con quale probabilità si possa considerare causato da esposizione un tumore che abbia molteplici determinanti.⁷ Per attribuire a un’esposizione particolare, al fumo di sigaretta per esempio, un tumore polmonare insorto in un fumatore si adotta spesso il criterio del “più probabile che non”, che corrisponde a dimostrare l’esistenza di un rischio relativo maggiore di 2 che corrisponde a un rischio attribuibile, negli esposti, >50%. Berry ha argomentato, analizzando il problema in termini di anticipazione dell’evento, che è da attribuire un ruolo causativo dell’esposizione all’evento insorto nell’esposto se solo RR>1 qualora l’incidenza aumentasse con l’età. In tal caso, infatti, per ogni caso esposto si può calcolare un tempo di anticipazione che riguarda il 100% dei casi insorti tra gli esposti. Berry esemplifica il suo argomento considerando la stima di anticipazione dell’evento (RAP) per un fumatore con tumore del polmone diagnosticato ad A anni di età con la formula:

RAP = (A – 20) x (RR1/γ – 1)

ottenuta sostituendo all’età il valore (età – 20) (considerando che i soggetti comincino a fumare in media all’età di 20 anni) e usando per γ il valore di 5,67, stima del parametro di relazione lineare fra Ln(Età – 20 anni) e Ln(Tasso di mortalità per tumore del polmone) ottenuta combinando studi di incidenza nelle popolazioni di UK, US e Australia. I valori A1 e A0 sono quelli per cui risulta: Ln(tasso nei non esposti) = 5,67Ln(A0-20) uguale a Ln(tasso negli esposti) = 5,67Ln(A1-20) x RR.
Limitandosi a considerare tumori insorti tra 40 e 75 anni, Berry riporta la tabella degli anni di anticipazione attesa per valori di RR che vanno da 1,01 a 3,0. Dalla formula, per un fumatore con tumore del polmone diagnosticato a 60 anni e con RR = 2 per il fumo, si ottiene, per esempio, la stima di anticipazione:

(60 – 20) x (2 1/5,67 – 1) = 5,2 anni

Discacciati, Bellavia, Orsini e Greenland,⁸ infine, considerando la relazione fra modello PH e modello AFT, evidenziano come il periodo di anticipazione (RAP):

1. non sia ricavabile dalla differenza fra medie di occorrenza in non esposti ed esposti, perché quest’ultima è influenzata dalla durata della osservazione;
2. non sia ricavabile come spostamento fra curve di sopravvivenza di esposti e non esposti;
3. non dipenda solo da RR ma anche da come l’incidenza aumenta con l’età. Infatti lo stesso valore di RR, con incrementi di incidenza per anno di età diversi, dà una diversa stima di anticipazione, tanto minore quanto maggiore è l’incremento di incidenza con l’età.

Da notare, per il terzo punto, che se il coefficiente associato all’età è zero, ci si trova in una singolarità, poiché il rapporto β_1/β₂ che abbiamo visto misurare la anticipazione nel modello di relazione precedentemente considerato, cresce oltre ogni limite. Ma evidentemente per una patologia la cui occorrenza non aumenti con l’età il concetto di anticipazione non ha nessun contenuto, senza bisogno di fare appello a rischi competitivi, che riguardano una diversa problematica di
stima. Come esempio si può pensare alla relazione fra esposizione a cloruro di vinile monomero e insorgenza di angiosarcoma epatico. La patologia è rarissima nella popolazione non esposta. Il considerarla, qualora si manifesti in un esposto, certamente insorta a causa dell’esposizione è, in questo caso, il senso di lettura di RAP → ∞.

Esposizione ed amianto, occorenza di mesotelioma e RAP

Consonni et al.⁹ hanno pubblicato una comunicazione sui risultati dell’analisi di dati di registri regionali di mesoteliomi pleurici per il periodo 2000-2016 limitandosi alle età 45-74 anni per le regioni Piemonte e Lombardia. Gli Autori stimano in 7 e 10 anni, rispettivamente, la durata del periodo di anticipazione dell’insorgenza di mesotelioma pleurico nei maschi con età alla diagnosi di 45 e di 63 anni. Essi, non disponendo di dati di esposizione ad amianto nella popolazione generale, hanno assunto, ragionevolmente, l’appartenenza al genere maschile come indicatore di esposizione ad amianto stimando i parametri del modello di regressione di Poisson:

Ln(tasso) = a + b₁ genere + b₂ Ln(Età)

codificando il genere, proxy della esposizione, in: 0 = femmine, 1 = maschi.
La stima di RAP che hanno ottenuto (una sottostima, visto che la variabile “genere”, nella sua modalità “femminile”, riflette in modo imperfetto la non esposizione ad amianto) è ottenuta con la formula:

RAP = Età(maschi) x exp[(b₁/b₂)-1]

appropriata alla scala logaritmica che hanno utilizzato sia per l’età sia per la variabile risposta.

Conclusioni

L’uso di modelli probabilistici è, al momento attuale, il nostro modo di modellare eventi che si presentano con modalità intrinsecamente variabili e non deterministicamente prevedibili. Questo non significa che i risultati che un tale modello fornisce non possano essere usati, fatte le appropriate verifiche del rispetto degli assunti di validità e nel campo appropriato di applicazione, per la decisione nel caso individuale, come già avviene nell’uso di evidenze per la decisione clinica nelle applicazioni della EBM. La validità, non la verità, di un modello è valutabile dall’accordo fra valori predetti e valori osservati. Quello che dobbiamo chiedere a un modello è quindi di risultare utile a fare avanzare il nostro percorso verso la verità, nel quale è implicito il riconoscimento delle relazioni matematiche fra modelli per le quali, nel nostro caso, la possibilità di convertire risultati ottenuti in termini di HR o RR, utili allo studio della relazione causale, in risultati espressi in termini di RAP, utili a esprimere l’impatto dell’esposizione sui singoli casi prodotti.

Bibliografia

1. Duca PG, Aurora F, Bai E et al. Mesotelioma: in cerca di posizioni condivise. Epidemiol Prev 2019;43(3-4):112-13.
2. Apostoli P, Boffetta P, Bovenzi M et al. Società Italiana di Medicina del Lavoro. Position Paper Amianto. Med Lav 2019;110: 459-85.
3. Barbieri PG et al. A proposito dell’amianto e del Position Paper della Società Italiana di Medicina del Lavoro. Epidemiol Prev 2020;44(1):89-91.
4. Cox DR. Regression models and life tables (with discussion). J R Statist Soc B 1972;34:187-220.
5. Collett D. Modelling survival data in medical research. London, Chapman & Hall, 1994.
6. Armitage P, Doll R. The age distribution of cancer and multistage theory of carcinogenesis. Br J Cancer 1954;8:1-12. (Comments in: Int J Epidemiol 33(6):1174-96).
7. Berry G. Relative risk and acceleration in lung cancer. Stat Med 2006;26:3511-7.
8. Discacciati A, Bellavia A, Orsini N, Greenland S. On the interpretation of risk and rate advancement periods. Int J Epidemiol 2016;45:278-84.
9. Consonni D, Migliore E, Barone-Adesi F et al. Gender differences in pleural mesothelioma occurrence in Lombardy and Piedmont, Italy. Environ Res 2019;177:108636177.

Conflitti di interesse dichiarati: l'autore è membro del direttivo di Medicina Democratica - Movimento di lotta per la Salute Onlus.