Un quadro metodologico per la sorveglianza dei tumori occupazionali da grandi archivi informativi: controllo della molteplicità e robustezza inferenziale

Introduzione

A partire dalla fine degli anni Novanta, con la costituzione di un registro nazionale dei tumori professionali presso l’INAIL (D.lgs. 81/2008), le banche dati amministrative hanno assunto un ruolo chiave nella sorveglianza dei rischi occupazionali, grazie all’uso di registri di mortalità, ricoveri ospedalieri e archivi occupazionali. In Italia, il record linkage tra flussi ha permesso di condurre studi di mortalità per settore lavorativo e il progetto OCCAM¹ è un esempio di integrazione di registri tumori, schede di dimissione ospedaliera e storie lavorative allo scopo di individuare tumori professionali. Nello stesso periodo, nel Regno Unito e nei Paesi nordici sono state avviate esperienze analoghe, basate su dati di mortalità collegati al censimento e coorti di lungo periodo.^2,3

Questi approcci hanno dimostrato l’utilità di studiare popolazioni ampie, includere esiti rari e monitorare i lavoratori per lunghi periodi, migliorando la qualità dell’informazione sulla loro storia espositiva.^4,5 Rimangono, tuttavia, alcune criticità legate alla frequente mancanza di informazioni su confondenti potenzialmente rilevanti, alla scarsa granularità o alla mancanza di dettagli sulle mansioni e alle possibili distorsioni derivanti da misclassificazione^.6,7 

Il Registro Nazionale dei Mesoteliomi (ReNaM)⁸ e il Registro Nazionale dei Tumori Naso-Sinusali (ReNaTuNS)⁹ rappresentano, nel nostro Paese, strumenti efficaci per la sorveglianza epidemiologica dei tumori professionali caratterizzati da un’elevata frazione eziologica, finalizzati all’individuazione di nuovi casi incidenti e all’identificazione di rischi occupazionali noti o emergenti. Per la maggior parte delle altre sedi tumorali, caratterizzate da una più bassa frazione eziologica (perché associate ad altri determinanti quali gli stili di vita e abitudini al fumo) e più alta incidenza, la metodologia di rilevazione necessita di un approccio sostanzialmente diverso. 

In questo contesto, il progetto OCCAM ha rappresentato un importante passo avanti verso la costruzione di un sistema di sorveglianza epidemiologica potenzialmente estendibile a livello nazionale, grazie all’impiego di fonti informative correnti, quali l’archivio della mortalità ISTAT e l’archivio delle storie contributive dell’INPS. Tuttavia, l’archivio INPS presenta alcune limitazioni rilevanti. In particolare, la copertura del database contributivo è circoscritta ai lavoratori del settore privato, escludendo i lavoratori autonomi, quelli del pubblico impiego, gli artigiani e i lavoratori domestici od occasionali. Inoltre, non sono disponibili informazioni dettagliate sulla mansione specifica svolta, ma soltanto sul settore economico di appartenenza: ciò rende attualmente non praticabili lo sviluppo e l’utilizzo di strumenti analitici come le matrici mansione-esposizione. L’impiego di procedure di record linkage tra archivi sanitari (mortalità per causa) e storie contributive costituisce uno strumento valido per il monitoraggio delle neoplasie a bassa frazione eziologica di sospetta origine professionale. Si tratta di un’opportunità di grande rilevanza, che richiede, tuttavia, un’attenta valutazione dei limiti interpretativi delle stime prodotte, spesso basate su liste di rischi relativamente imperfette, derivate da analisi di mortalità proporzionale per tumori e con la possibilità di controllare solo un numero limitato di confondenti quali età, livello di istruzione e area geografica.

Questi limiti rendono necessaria l’adozione di strategie metodologiche rigorose per evitare fenomeni di data dredging e per garantire la solidità delle inferenze.^10,11

In questo contributo, viene presentata una strategia sistematica per la valutazione dei risultati della sorveglianza dei tumori professionali condotta attraverso record linkage di archivi nazionali di mortalità e di storia contributiva. Nel dettaglio, il presente studio si caratterizza per l’impiego sistematico dei grandi archivi dell’Istituto Nazionale per l’Assicurazione contro gli Infortuni sul Lavoro (INAIL) all’interno del progetto Big Data e Deep Learning nella Sorveglianza dei Tumori Professionali (BEST). INAIL, infatti, grazie a un progetto approvato nel Piano Statistico Nazionale, che prevede l’autorizzazione preventiva con forza di legge dell’Autorità Italiana per la Protezione dei Dati per le opere statistiche di interesse pubblico,¹² ha potuto realizzare tramite record linkage l’unione del Registro Nazionale delle Cause di Morte (ReNCaM) dell’ISTAT 2005-2018 con l’archivio delle storie contributive (INPS) 1974-2018. Da questi grandi archivi, sono stati stimati rischi relativi per sede tumorale e settore occupazionale e adottate metodologie specifiche per affrontare i problemi della molteplicità dei confronti e del confondimento non misurato, che prevedono l’utilizzo di soglie in grafici quantile-quantile per identificare osservazioni aberranti, q-value per la selezione di insiemi di osservazioni da approfondire e modelli gerarchici bayesiani per la costruzione di graduatorie. Sono stati inoltre calcolati gli E-value per stimare l’influenza potenziale di confondenti non misurati.^13,14 

Obiettivi dello studio

Obiettivi di questo studio sono rafforzare la capacità di individuare i rischi oncologici legati al lavoro, descriverne in modo accurato l’incertezza e sviluppare strumenti innovativi per comunicare il rischio alle autorità sanitarie locali, contribuendo così alla prevenzione delle malattie professionali e al miglioramento delle condizioni lavorative. In particolare, si intende affiancare alla stima dei rischi tumorali nei diversi settori occupazionali una gestione adeguata dell’incertezza dovuta alla molteplicità delle analisi e ai potenziali fattori di confondimento, offrendo strumenti critici utili a orientare le iniziative di prevenzione e a definire le priorità di intervento.

Materiali e metodi

Disegno di studio, dati e analisi statistica

Lo studio utilizza un disegno di mortalità proporzionale per cancro (proportional cancer mortality) basato su dati di mortalità. È  più utile considerare gli studi di mortalità proporzionale come studi caso-controllo in cui i casi sono decessi per la causa di interesse e i controlli sono neoplasie selezionate non associate all’esposizione di interesse.^15-17 La misura di associazione utilizzata è il mortality odds ratio (MOR), che rappresenta una stima del rapporto fra tassi di incidenza. La restrizione alle cause neoplastiche per casi e controlli è in genere suggerita per affrontare l’effetto lavoratore sano (healthy worker effect).^15,16

In questo studio, è stata effettuata una serie di analisi caso-controllo su 9 neoplasie di interesse e la loro associazione con l’impiego in 43 settori occupazionali. Per ogni combinazione neoplasia-settore, è stato stimato un cancer mortality odds ratio (CMOR). In questo contesto, sarebbe difficilmente fattibile scegliere, per ogni analisi, una serie di neoplasie di controllo non associate all’esposizione. Pertanto, nell’analizzare una certa neoplasia, sono stati selezionati come controlli tutte le altre neoplasie. Questo approccio è considerato valido, perché ci sono pochi fattori che aumentano sensibilmente il rischio di numerose malattie;¹⁵ in altre parole, se la serie di controllo contiene una o più neoplasie associate all’esposizione, il loro effetto viene “diluito” dalle altre neoplasie non associate all’esposizione. 

Nel presente studio, sono stati analizzati i decessi dal ReNCaM per gli anni 2005-2018. Le informazioni relative al settore lavorativo prevalente derivano dall’Archivio nazionale delle storie contributive dell’INPS per il periodo 1974-2018. Il record linkage tra i dati di mortalità e quelli occupazionali è stato effettuato a livello individuale, mediante codice fiscale, dall’INAIL. Il dataset risultante è stato anonimizzato e trattato nel rispetto delle disposizioni del Regolamento europeo 2016/679 (GDPR) e del D.Lgs. 196/2003, come modificato dal D.Lgs. 101/2018, in materia di protezione dei dati personali.

Poiché gli archivi INPS coprono circa il 55% della forza lavoro italiana, lo studio include esclusivamente i lavoratori dipendenti del settore privato, unici soggetti per i quali sono disponibili informazioni contributive complete. Restano, pertanto, esclusi i lavoratori autonomi, il pubblico impiego, gli artigiani e i lavoratori domestici od occasionali, categorie non presenti negli archivi INPS. Non si dispone dell’informazione della qualifica professionale, ma viene indicato solo se il lavoratore ha svolto la mansione in termini di operaio.

Per ciascun individuo, è stato ricostruito il percorso lavorativo a partire dal 1974, classificando il settore di impiego secondo la codifica ATECO,¹⁸ successivamente armonizzata in 48 categorie economiche sulla base della classificazione statistica delle attività economiche nella Comunità Europea.¹⁹

L’esposizione di interesse è stata definita in relazione al settore occupazionale prevalente in cui ciascun lavoratore ha svolto un’attività, ossia quello in cui ciascun soggetto ha maturato il maggior numero di anni di contribuzione. Ai fini dell’analisi, è stato applicato un periodo di latenza di cinque anni antecedente al decesso.

L’analisi è stata condotta sugli uomini senza distinzione di mansione. Per ciascuna combinazione settore-tumore, sono stati stimati i CMOR con relativi intervalli di confidenza al 90%, mediante modelli di regressione logistica non condizionata. Come gruppo di riferimento, sono stati scelti i lavoratori dei settori dei servizi (alberghi e consumazioni di bevande, commercio, attività connesse ai trasporti, riparazioni non altrove classificate e altre attività di servizi). Le stime sono state corrette per anno del decesso, età al decesso (classi quinquennali), livello di istruzione, classificato in tre categorie (basso: istruzione elementare o media inferiore; medio: istruzione secondaria superiore; alto: istruzione universitaria) e ultima regione di residenza. 

Sono state, inoltre, condotte analisi di sensibilità, rimuovendo il vincolo di latenza o modificando la definizione di non esposto (tutti i settori, anziché il settore servizi), per verificare la robustezza dei risultati (risultati non presentati).

L’analisi si è focalizzata su nove sedi neoplastiche: laringe, naso e seni paranasali, rinofaringe, polmone, mesotelioma (pleura, peritoneo e pericardio), melanoma cutaneo, pelle non melanoma, tessuti linfatici ed ematopoietici, vescica. Per queste sedi, la componente occupazionale risulta chiaramente documentata dalle classificazioni della IARC (gruppi 1 e 2A) e dell’INAIL (Liste I e II). Queste analisi danno luogo a un’ampia matrice di CMOR, che incrocia i comparti occupazionali con le diverse sedi tumorali. In questo contributo, l’attenzione si concentra sugli aspetti metodologici connessi alla valutazione dell’incertezza dei CMOR, affrontando in particolare i problemi legati alla molteplicità dei confronti, alla costruzione di graduatorie e al confondimento non misurato.

Incertezza sulla molteplicità dei rischi stimati

L’incertezza connessa alla molteplicità dei rischi stimati^20-22 ha richiesto l’elaborazione di una strategia specifica; perciò, sono stati adottati tre approcci che rispondono a tre obiettivi inferenziali distinti.

1. Identificazione delle osservazioni aberranti: grafici quantile-quantile con soglie (Q-Q plot con guide rails). Per ogni rischio relativo stimato (CMOR) è possibile calcolare un test statistico (z-score) rapportandolo al suo errore standard. Sotto ipotesi nulla, è lecito confrontare tutti gli z-score così calcolati con la distribuzione normale teorica e riportare un grafico con i quantili della distribuzione osservata degli z-score in ordinata e in ascissa i quantili della distribuzione teorica normale. Per identificare le osservazioni aberranti, si considera la probabilità del valore massimo di n z-score indipendenti, che sarà la probabilità dell’evento congiunto che tutti gli z-score siano inferiori a un determinato valore soglia. Per la regola del prodotto per eventi indipendenti, vale che questa probabilità è data dal prodotto di n identiche distribuzioni cumulative normali. Quindi:

Pr(Z_max ≤ z) = {Φ(z)}ⁿ

da cui la soglia cercata:

1 - α = {Φ(z)}ⁿ

z = Φ-1 {(1 - α)^1/n)}

Si ricorda che questa formula dipende dalla validità dell’assunzione di normalità della distribuzione degli z-score empirici. Inoltre, si assume che la pendenza del grafico quantile-quantile sia pari a 1, ovvero la bisettrice. In altre parole, se ci sono molte osservazioni aberranti, ci si trova in presenza di sovradispersione e questo implicherebbe una pendenza diversa da 1. Law et al. suggeriscono soglie derivate da un’opportuna simulazione, quindi di stimare la sovradispersione;²¹ qui si utilizza il grafico quantile-quantile nella forma più semplice per individuare pochi punti aberranti, cioè le osservazioni che superano soglie critiche definite (α=0,05; α=0,10).²¹ Il grafico, quindi, confronta gli z-score osservati con quelli attesi sotto ipotesi nulla. I punti, relativi a determinati settori lavorativi, che superano le guide rails segnalano la situazione di rischi potenzialmente rilevanti.

2. Selezione del gruppo di osservazioni meritevoli di approfondimento: controllo del false discovery rate (FDR) e q-value. Questo approccio consente di controllare al livello desiderato la proporzione attesa di falsi positivi tra i risultati considerati significativi.^22,23 Il q-value può essere ottenuto come una semplice trasformazione dei p-value. La trasformazione è monotona, quindi l’ordinamento dei p-value dal più piccolo al più grande equivale all’ordinamento dei q-value. In formula: 

dove:

la notazione (i) indica i valori ordinati dal più piccolo al più grande;

m è il numero totale di test;

π₀ è la probabilità a priori dell’ipotesi nulla, che senza perdere in generalità è stata fissata pari a uno.

Avendo un p-value per ogni stima di rischio relativo, si dispone anche del corrispettivo q-value, che è possibile usare euristicamente come era stato suggerito per i p-value da Fisher. L’interpretazione, però, è diversa: se si dichiarano significativi tutti i rischi relativi con q-value inferiore a, ipotizziamo, 0,05, questo implica che il 5% di questi saranno dei falsi positivi. Per esempio, se vengono identificati 20 rischi con q-value <0,05, ci si attende la presenza di un solo falso positivo, anche se non è possibile determinarne l’identità. È stato deciso di riportare in forma grafica i CMOR tramite caterpillar plot, uno strumento grafico che permette di visualizzare, per ciascun settore occupazionale e tumore, le stime di rischio con i relativi intervalli di confidenza. Per coerenza con i q-value, si corregge l’ampiezza dell’intervallo di confidenza in modo da controllare il false coverage rate allo stesso livello scelto per il FDR. In formula, significa calcolare l’intervallo di confidenza al livello:

α' = α × R/m

dove:

R è il numero di rischi dichiarati significativi dall’ispezione dei q-value;

m il numero di rischi stimati;

α il livello di confidenza scelto.

Questa rappresentazione facilita il confronto immediato tra comparti e sedi tumorali, ordinandoli dal livello di rischio più basso a quello più elevato. 

3. Graduatorie dei settori occupazionali e dei tumori professionali attraverso i ranghi e la loro incertezza. Dato un parametro qij di interesse, il rischio relativo di mortalità per tumore nel comparto i e sede tumorale j, è possibile definire il rango R_ij come:

dove I è la funzione indicatrice che assume valore 1 se la condizione è verificata e 0 altrimenti. In altre parole, il rango R_ij rappresenta la posizione relativa del rischio stimato Θ_ij rispetto agli altri rischi della stessa serie (per esempio, le diverse sedi tumorali all’interno di uno stesso comparto). Pertanto, la struttura dei dati è organizzata in una matrice IxJ in cui le righe corrispondono ai “profili” di rischio per comparto e le colonne alle “firme” tumorali per sede. A partire da questa matrice, è possibile derivare matrici di ranghi Rij che consentono di ordinare, in termini probabilistici, i comparti o le sedi tumorali in funzione del rischio stimato e della relativa incertezza.

Per il calcolo dei ranghi e della relativa incertezza, è stato adottato un modello gerarchico bayesiano. In particolare, si assume che, per ciascun settore, il vettore dei log(CMOR) segua una distribuzione normale multivariata MVN(µ_i, ∑), dove µi rappresenta il vettore medio dei rischi attesi per il settore i-esimo e ∑ la matrice di covarianza comune tra le diverse sedi tumorali. Sui parametri di queste distribuzioni, sono state adottate a priori non informative: distribuzioni normali a media zero e precisione 0,0001 per ciascun elemento di mi e distribuzione di Wishart W(Ω,ν) per ∑. Maggiori dettagli in Catelan et al.²⁴

Le distribuzioni a posteriori dei ranghi Rij sono state ottenute mediante simulazione Markov Chain Monte Carlo (MCMC) utilizzando WinBUGS 1.4.²⁵

La media a posteriori dei ranghi fornisce una misura sintetica dell’ordinamento dei comparti o delle sedi tumorali per rischio. L’incertezza viene riportata con intervalli di credibilità all’80%.^24,26 

Sensibilità delle associazioni rispetto ai confondenti non misurati

Per valutare, invece, la sensibilità delle associazioni rispetto a confondenti non misurati (per esempio, l’abitudine al fumo o il consumo di alcol), è stato utilizzato l’E-value, una metrica proposta da VanderWeele e Ding.14 L’E-value quantifica la forza che un confondente non misurato dovrebbe avere sia con l’esposizione (settore lavorativo) sia con l’esito (tumore) per spiegare interamente l’associazione osservata. Per rischi relativi (RR), l’E-value si calcola come: 

E - value = RR + √ RR x (RR – 1)

Questo approccio può essere esteso, con opportune modifiche, anche ad altre misure di effetto. Nel caso di odds ratio (OR) o hazard ratio (HR), se l’esito è raro (per esempio, < 15%), è possibile applicare la stessa formula. Tuttavia, quando l’esito è più frequente, nella formula sopra indicata, il RR deve essere sostituito con la radice quadrata dell’OR o dell’HR.¹⁴ Inoltre, qualora l’effetto osservato sia protettivo, nella formula si utilizza l’inverso della misura di effetto. 

Un E-value elevato è tipicamente indicazione di robustezza dell’associazione osservata, quindi difficilmente attribuibile solo a fattori non misurati; un E-value basso segnala, invece, maggiore vulnerabilità del risultato. Nello studio, gli E-values sono stati ottenuti a partire dai CMOR, limitatamente per i comparti con un eccesso di rischio (CMOR ≥1), insieme ai valori corrispondenti (associazione con l’esposizione RR_EU e associazione con l’esito RR_UD) stimati per alcune delle covariate osservate (istruzione e anno del decesso). Ciò ha reso possibile confrontare la solidità delle associazioni osservate rispetto a fattori noti, considerandoli, di volta in volta, come non misurati.²⁷ Le analisi sono state condotte mediante il software R usando il package “Evalue”.²⁸ 

Risultati

In tabella 1 sono riportati numero e percentuale di lavoratori deceduti per sede primaria del tumore, secondo la classificazione ICD-10.

La quota più elevata riguarda il tumore del polmone, che rappresenta oltre un quarto del totale (25,8%). Alcune sedi rare, ma di forte interesse occupazionale, come i tumori della pleura (1,4%) e quelli naso-sinusali (0,1%), pur rappresentando una piccola frazione dei casi, rivestono particolare rilevanza per la sorveglianza epidemiologica in ragione delle esposizioni professionali ben documentate (per esempio, amianto, polveri di legno e di cuoio). Più della metà dei decessi (56,9%) ricade nella categoria “altre sedi tumorali”.

Nella tabella S1 (vedi materiali supplementari online) sono riportati i CMOR per i settori e le neoplasie scelte. Gli eccessi più consistenti sono stati osservati per il tumore del polmone, il mesotelioma e i tumori naso-sinusali. L’analisi ha messo in evidenza differenze marcate tra settori occupazionali nella mortalità per tumore polmonare. In particolare, il settore dell’edilizia ha mostrato un rischio aumentato (CMOR 1,14; IC90% 1,13-1,16), così come il settore della pesca (CMOR 1,26; IC90% 1,19-1,33). Rischi inferiori al valore nullo sono stati osservati per il settore dell’agricoltura (CMOR 0,88; IC90% 0,87-0,89). Per i tumori della pleura e del peritoneo, l’associazione più forte è stata rilevata tra i settori dei cantieri navali (CMOR 6,36; IC90% 5,72-7,07), delle costruzioni meccaniche (CMOR 1,76; IC90% 1,67-1,86) e della produzione di materiali per l’edilizia (CMOR 1,82; IC90% 1,63-2,02).  Sebbene rari, anche i tumori delle cavità nasali hanno mostrato CMOR elevati in specifici settori occupazionali. In particolare, il settore dell’industria della pelle e delle calzature (CMOR 4,13; IC90%: 3,05-5,59), così come il settore del legno (CMOR 3,98; IC90% 3,15-5,02).

In figura 1 sono mostrati i due ranking marginali.

La figura 1A rende evidente che i mesoteliomi risultano al primo posto tra tutte le categorie occupazionali, seguiti dai tumori cutanei. Le neoplasie del rinofaringe, della laringe, naso-sinusali, polmonari ed ematopoietiche mostrano una distribuzione più sovrapposta. Meno rilevanti, considerando l’insieme delle categorie occupazionali, risultano i tumori della vescica e il melanoma cutaneo. 

La figura 1B presenta il ranking marginale delle categorie occupazionali rispetto a tutte le neoplasie esaminate. È stata osservata una notevole sovrapposizione, come indicato dall’ampiezza degli intervalli di credibilità, con l’eccezione dell’agricoltura, stabilmente collocata all’ultima posizione. Questo pattern riflette la presenza di fattori di rischio occupazionali specifici per ciascuna categoria, non condivisi tra i diversi settori. Ogni categoria occupazionale presenta, infatti, un proprio rischio tumorale peculiare e, pertanto, nessuna emerge in modo netto quando si marginalizza su tutte le sedi tumorali.

Il QQ-plot (figura 2A), basato sugli z-score, combina la stima puntuale dell’OR con il relativo errore standard e ha individuato come categorie occupazionali fuori soglia: costruzioni, trasporti, pesca e ristorazione. L’inferenza selettiva basata sul controllo della false coverage (figura 2B) ha posto l’attenzione su: pesca, rifiuti, ceramica, trasporti marittimi, pulizie, trasporti, ristorazione e costruzioni. Questi schemi, pur largamente coerenti tra loro, dipendono da un diverso peso attribuito alla stima puntuale dell’OR, all’errore standard e all’eterogeneità tra categorie occupazionali (ovvero il numero di OR divergenti).

I ranking rappresentano uno strumento utile per descrivere la graduatoria complessiva dei rischi, esplorare profili specifici o sintetizzare un quadro generale. Le prime posizioni forniscono indicazioni immediate su quali possono essere le priorità di intervento. La matrice degli OR può essere analizzata in due modi: considerando una sede tumorale specifica (profilando le categorie occupazionali) oppure una categoria occupazionale (analizzandone la “firma” per sede tumorale). Per esempio, il profilo delle categorie occupazionali per il tumore del polmone (figura 3A) mette in luce il rischio più elevato nel settore della pesca, seguiti da quelli dei rifiuti, ceramica, parrucchieri, costruzioni, trasporti, ristorazione e pulizie. Un’ulteriore applicazione è rappresentata dalla “firma” per sede tumorale della categoria delle costruzioni (figura 3 B), in cui i rischi per tumore della laringe e tumore del polmone emergono come i più rilevanti.

La tabella 2 riporta i CMOR aggiustati e i corrispondenti valori di E-value per il tumore polmonare, escludendo i settori occupazionali per i quali è stato osservato un CMOR inferiore a 1.

Il valore di E-value più elevato si osserva nel comparto della pesca (1,49), seguito da parrucchieri (1,44) e raccolta e smaltimento rifiuti (1,39). Oltre a questi settori, i risultati di questo studio hanno permesso di osservare eccessi di rischio più modesti nei settori dell’edilizia (1,34) e dei trasporti (1,29). Queste evidenze risultano particolarmente rilevanti, poiché entrambi i comparti sono storicamente associati a cancerogeni occupazionali ben noti. Analogamente, CMOR elevati sono stati ottenuti nei settori di ristorazione e servizi di pulizia e disinfestazione (entrambi con E-value 1,34), sebbene le possibili fonti di esposizione siano meno chiare e possano riflettere in parte un confondimento legato a fattori socioeconomici.

Valori di E-value inferiori, invece, sono stati ottenuti in relazione ai settori occupazionali di gomma, chimica, industria del tabacco e petrolio (E-values compresi tra 1 e 1,08).  La tabella 2 riporta il valore dell’E-value sia per la stima del CMOR sia per l’estremo inferiore del relativo intervallo di confidenza. È possibile, quindi, dedurre che l’intervallo di confidenza per molti dei comparti lavorativi riportati comprende il valore 1 e, di conseguenza, l’E-value per il limite inferiore dell’intervallo è pari a 1, indicando la non necessità di avere confondenti non misurati per portare a nulla l’associazione osservata. Per esempio, per il tumore al polmone, il settore dei laboratori fotografici presenta un CMOR pari a 1,13 (IC90% 0,98-1,29); il corrispondente valore di E-value per la stima puntuale è 1,32, ma per l’estremo inferiore dell’IC si ottiene il valore 1.

In tabella 2 sono presentati anche i valori delle associazioni di istruzione e anno del decesso sia con il comparto lavorativo sia con il tumore polmonare (RREU e RRUD). Questi contribuiscono a valutare la robustezza delle associazioni osservate. Per i settori considerati, in linea con il tipo di studio e con la definizione generica di esposizione adottata, le stime delle associazioni e, di conseguenza, i valori di E-value, sono risultate modeste; tuttavia, l’entità dei valori dei confondenti misurati è simile o inferiore, suggerendo una certa robustezza per le stime ottenute.

Discussione

L’integrazione degli archivi nazionali di mortalità e di contribuzione previdenziale ha garantito una numerosità elevata, che ha reso possibile la valutazione dei rischi associati a tumori rari, ma di particolare interesse per l’epidemiologia occupazionale. L’impiego di banche dati amministrative, caratterizzate da elevata copertura, rappresenta pertanto un pilastro fondamentale per la sorveglianza dei tumori di origine professionale.⁷ 

I risultati qui presentati, in accordo con l’evidenza epidemiologica, confermano un aumento del rischio di tumore polmonare tra i lavoratori edili, di mesotelioma tra gli addetti ai cantieri navali e di tumori naso-sinusali tra i lavoratori dell’industria del cuoio e del legno. Queste associazioni risultano coerenti con le esposizioni documentate a cancerogeni noti, tra cui amianto, silice cristallina e polveri di cuoio e di legno.^29,30

Oltre a questi rischi, per larga parte attesi, sono stati osservati rischi aumentati in settori non identificati come prioritari per la sorveglianza, come il settore dei servizi di pulizia e di ristorazione. Questi risultati richiedono ulteriori approfondimenti, poiché potrebbero riflettere sia esposizioni occupazionali emergenti sia distorsioni da selezione (legate al disegno con controlli deceduti per neoplasia) o da informazione (legate alla definizione dell’esposizione basata sul settore economico) oppure da imperfetta misura dei confondenti o, come in tutti gli studi osservazionali, da fattori confondenti non misurati.³¹

Tra i lavoratori agricoli, si osserva un decremento del rischio che potrebbe essere in parte attribuito all’healthy worker effect, ma anche a fattori protettivi legati allo stile di vita o all’ambiente. Questa interpretazione, tuttavia, deve essere considerata con la dovuta cautela, alla luce dei limiti dello studio e anche della natura eterogenea delle esposizioni in agricoltura, considerando che questo settore risulta comunque esposto ad altri agenti (per esempio, pesticidi) con effetti sulla salute riportati in un’estesa letteratura.^32,33

Accanto ai punti di forza, vanno riconosciute alcune limitazioni. Gli archivi amministrativi non consentono una ricostruzione dettagliata delle esposizioni: pur disponendo di informazioni sul settore industriale, mancano dettagli sulle mansioni specifiche e sull’intensità dell’esposizione, con il conseguente rischio di misclassificazione. Inoltre, non sono disponibili informazioni individuali su fattori come il fumo, il consumo di alcol, le comorbidità e la familiarità, che rappresentano potenziali confondenti per molte neoplasie. Infine, non è stato possibile considerare i tempi di latenza specifici delle diverse sedi tumorali.

Dal punto di vista metodologico, il nostro approccio si colloca in un filone innovativo volto a rafforzare la sorveglianza e a definire priorità preventive basate sull’evidenza. La questione dei test multipli, raramente affrontata nell’epidemiologia occupazionale, assume un ruolo centrale quando si analizzano numerosi settori e sedi tumorali. I metodi tradizionali di controllo del Family Wise Error Rate (FWER), come la correzione di Bonferroni,³⁴ risultano troppo conservativi e riducono drasticamente la potenza statistica. Un’alternativa è rappresentata dal FDR13 e dal relativo q-value,²² ormai strumenti consolidati per l’interpretazione dei risultati.

In tale contesto, i metodi di ranking basati su modelli bayesiani gerarchici si sono dimostrati particolarmente efficaci nel sintetizzare ampie serie di associazioni occupazione-patologia, consentendo di identificare combinazioni settore-tumore a rischio aumentato e, al contempo, di quantificare l’incertezza. Questo approccio, già applicato nell’epidemiologia ambientale,²⁴ è stato qui esteso ai tumori occupazionali, in cui i problemi di test multipli e di inferenza selettiva sono particolarmente rilevanti.²¹ L’integrazione di ranking bayesiano, Q-Q plot e controllo del FDR ha rafforzato la robustezza del processo di prioritizzazione, migliorandone l’interpretabilità.

Parallelamente a ciò, il potenziale impatto sulle stime del confondimento non misurato è stato valutato attraverso il calcolo degli E-values. Questo approccio, introdotto inizialmente per il rischio relativo, è stato esteso a diverse misure di effetto¹⁴ e ha permesso una più accurata interpretazione delle associazioni osservate. Tuttavia, è necessario mettere in luce alcune limitazioni che si accompagnano all’utilizzo dell’E-value. In primo luogo, l’interpretazione risulta meno immediata per esposizioni continue, in cui contrasti più ampi generano valori di E-value più elevati, senza necessariamente implicare una maggiore robustezza dell’associazione.³⁵ Oltre a ciò, in presenza di molteplici confondenti non misurati, alcune associazioni possono superare il valore calcolato, richiedendo l’interpretazione con cautela anche degli E-values elevati.³⁶ Inoltre, il metodo si fonda su assunzioni semplificative, come l’uguaglianza e la concordanza delle associazioni confondente-esposizione e confondente-esito.³⁵ Nei casi di effetti apparentemente protettivi, sono necessari ulteriori adattamenti: VanderWeele³⁷ suggerisce di affrontare questa questione considerando come esposizione di interesse la categoria caratterizzata da un rischio maggiore. Infine, una revisione degli studi applicativi³⁸ mette in guardia contro un uso meramente meccanico dell’E-value, sottolineando la necessità di contestualizzare i risultati rispetto ai possibili confondenti non misurati e, quando appropriato, di integrarli con analisi di sensibilità alternative.³⁷ Nel presente studio, l’applicazione dell’E-value ha permesso di quantificare la robustezza delle associazioni rispetto ai potenziali confondenti non misurati. I valori ottenuti sono risultati modesti, in linea con il disegno dello studio e con le limitazioni legate alla definizione dell’esposizione. Inoltre, l’aggiustamento non ha potuto tenere conto di confondenti rilevanti legati allo stile di vita individuale come fumo e consumo di alcol. Tuttavia, nei casi in cui l’entità delle associazioni dei confondenti misurati con esposizione ed esito è risultata paragonabile o inferiore agli E-values, vi è un supporto della robustezza delle stime ottenute.

Conclusioni

L’integrazione di strumenti grafici di diagnostica, del controllo del FDR, del ranking bayesiano ed E-values delinea un quadro metodologico innovativo nell’ambito dell’epidemiologia dei tumori occupazionali, estendendo approcci già applicati in ambito di salute ambientale.^14,22,24 

Sono stati impiegati strumenti grafici e statistici per esaminare un esteso insieme di associazioni tra sedi tumorali e settori occupazionali, con l’obiettivo di identificare e prioritizzare i segnali emergenti, valutando al contempo l’influenza di potenziali confondenti non misurati. L’applicazione integrata di questi approcci ha offerto un esempio concreto di implementazione e interpretazione metodologica, mettendone in evidenza la rilevanza in contesti applicativi reali.

In definitiva, questo quadro integrato contribuisce a rafforzare la solidità e l’interpretabilità dei risultati, supportando la definizione di priorità basate sull’evidenza per la sorveglianza della salute occupazionale.

Conflitti di interesse dichiarati: Dario Consonni è stato consulente del Giudice in procedimenti penali riguardanti malattie da amianto.

Bibliografia

1. Crosignani P, Massari S, Audisio R et al. The Italian surveillance system for occupational cancers: characteristics, initial results, and future prospects. Am J Ind Med 2006;49(9):791-98.
2. Harrington J. Reviews: Occupational Health Decennial Supplement. Health Education J 1995;54:381-82.
3. Kauppinen T, Heikkilä P, Plato N et al. Construction of job-exposure matrices for the Nordic Occupational Cancer Study (NOCCA). Acta Oncologica 2009;48(5):791-800. doi: 10.1080/02841860902718747
4. Registries for Evaluating Patient Outcomes: A User’s Guide. Rockville (MD), Agency for Healthcare Research and Quality, 2020.
5. Thygesen LC, Ersbøll AK. When the entire population is the sample: strengths and limitations in register-based epidemiology. Eur J Epidemiol 2014;29(8):551-58. doi: 10.1007/s10654-013-9873-0
6. Blair A, Stewart P, Lubin JH, Forastiere F. Methodological issues regarding confounding and exposure misclassification in epidemiological studies of occupational exposures. Am J Ind Med 2007;50(3):199-207. doi: 10.1002/ajim.20281
7. Kauppinen T. Exposure assessment – a challenge for occupational epidemiology. Scand J Work Environ Health 1996;22(6):401-3. doi: 10.5271/sjweh.160
8. INAIL. Registro Nazionale dei Mesoteliomi: ottavo rapporto. Roma, INAIL, 2025. Disponibile all’indirizzo: https://www.inail.it/portale/it/inail-comunica/pubblicazioni/catalogo-generale/catalogo-generale-dettaglio.2025.02.registro-nazionale-dei-mesoteliomi--ottavo-rapporto.html
9. INAIL. ReNaTuNS sorveglianza epidemiologica dei tumori naso-sinusali. Manuale operativo. Roma, INAIL, 2020. Disponibile all’indirizzo: https://www.inail.it/portale/it/inail-comunica/pubblicazioni/catalogo-generale/catalogo-generale-dettaglio.2020.11.renatuns-sorveglianza-epidemiologica-dei-tumori-naso-sinusali-manuale-operativo.html
10. Mehlum IS, Turner, MC. Challenges of Large Cohort and Massive Data in Occupational Health. In: Wahrendorf M, Chandola T, Descatha A (eds). Handbook of Life Course Occupational Health. Cham, Springer International Publishing, 2023; pp. 95-120. doi:10.1007/978-3-031-30492-7_3
11. Benjamini Y. Selective Inference: The Silent Killer of Replicability. Harvard Data Science Review 2. 2020. doi: https://doi.org/10.1162/99608f92.fc62b261
12. Sistema Statistico Nazionale (SISTAN). Programma statistico nazionale 2023-2025. 2024. Disponibile all’indirizzo: https://www.sistan.it/index.php?id=688 (2024).
13. Benjamini Y, Hochberg Y. Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. J R Stat Soc Series B Methodol 1995;57:289-300.
14. VanderWeele TJ, Ding P. Sensitivity Analysis in Observational Research: Introducing the E-Value. Ann Intern Med 2017;167(4):268-74. doi: 10.7326/M16-2607
15. Checkoway H, Pearce N, Kriebel D. Research Methods in Occupational Epidemiology. Oxford University Press 2004).
16. Miettinen OS, Wang JD. An alternative to the proportionate mortality ratio. Am J Epidemiol 1981;114(1):144-48. doi: 10.1093/oxfordjournals.aje.a113161
17. Richardson DB, Keil AP, Tchetgen Tchetgen E, Cooper G. Negative Control Outcomes and the Analysis of Standardized Mortality Ratios. Epidemiology 2015;26(5):727-32. doi: 10.1097/EDE.0000000000000353
18. ISTAT. Classificazione delle attività economiche ATECO. Roma, ISTAT, 2025. Disponibile all’indirizzo: https://www.istat.it/classificazione/classificazione-delle-attivita-economiche-ateco/
19. Eurostat. NACE Rev. 2 - Statistical Classification of Economic Activities. Lussemburgo, Office for Official Publications of the European Communities, 2008.
20. Carpenter LM, Swerdlow AJ, Fear NT. Mortality of doctors in different specialties: findings from a cohort of 20000 NHS hospital consultants. Occup Environ Med 1997;54(6):388-95. doi: 10.1136/oem.54.6.388
21. Law GR, Cox DR, Machonochie NE, Simpson J, Roman E, Carpenter LM. Large tables. Biostatistics 2001;2(2):163-71. doi: 10.1093/biostatistics/2.2.163
22. Storey JD. The positive false discovery rate: a Bayesian interpretation and the q-value. Ann Stat 2003;31:2013-35.
23. Benjamini Y, Yekutieli D. False Discovery Rate – Adjusted Multiple Confidence Intervals for Selected Parameters. J Am Stat Assoc 2005;100(469):71-81. doi: 10.1198/016214504000001907
24. Catelan D, Biggeri A. A statistical approach to rank multiple priorities in environmental epidemiology: an example from high-risk areas in Sardinia, Italy. Geospat Health 2008;3(1):81-89. doi: 10.4081/gh.2008.234
25. Spiegelhalter DJ, Thomas A, Best N, Lunn D. WinBUGS User Manual. Cambridge, MRC Biostatistics Unit, 2003.
26. Shen W, Louis TA. Triple-goal Estimates in Two-stage Hierarchical Models. J R Stat Soc Ser B Stat Methodol 1998;60:455-71.
27. Ding P, VanderWeele TJ. Sensitivity Analysis Without Assumptions. Epidemiology 2016;27(3):368-77. doi: 10.1097/EDE.0000000000000457
28. R core Team. R: A Language and Environment for Statistical Computing. Vienna (Austria), R Foundation for Statistical Computing, 2023.
29. International Agency for Research on Cancer. Agents Classified by the IARC Monographs, Volumes 1-133. Lione, IARC, 2023. Disponibile all’indirizzo: https://monographs.iarc.who.int/agents-classified-by-the-iarc
30. Mirabelli D, Kauppinen T. Occupational exposures to carcinogens in Italy: an update of CAREX database. Int J Occup Environ Health 2005;11(1):53-63. doi: 10.1179/oeh.2005.11.1.53
31. Gariazzo C, Massari S, Consonni D, Marchetti MR, Marinaccio A. Cancer-Specific Mortality Odds Ratios in the Food, Accommodation, and Beverage Activities in Italy. J Occup Environ Med 2024;66(7):572-79. doi: 10.1097/JOM.0000000000003114
32. Blair A, Ritz B, Wesseling C, Freeman LB. Pesticides and human health. Occup Environ Med 2015;72(2):81-82. doi: 10.1136/oemed-2014-102454
33. Togawa K, Leon ME, Lebailly P et al. Cancer incidence in agricultural workers: Findings from an international consortium of agricultural cohort studies (AGRICOH). Environ Int 2021;157:106825. doi: 10.1016/j.envint.2021.106825
34. Sedgwick P. Multiple significance tests: the Bonferroni correction. BMJ 2012;344:e509. doi:10.1136/bmj.e509
35. Ioannidis JPA, Tan YJ, Blum MR. Limitations and Misinterpretations of E-Values for Sensitivity Analyses of Observational Studies. Ann Intern Med 2019;170(2):108-11. doi: 10.7326/M18-2159
36. VanderWeele TJ, Ding P, Mathur M. Technical Considerations in the Use of the E-Value. J Causal Infer 2019;7. doi: 10.1515/jci-2018-0007
37. VanderWeele TJ. Are Greenland, Ioannidis and Poole opposed to the Cornfield conditions? A defence of the E-value. Int J Epidemiol 2022;51(2):364-71. doi: 10.1093/ije/dyab218
38. Blum MR, Tan YJ, Ioannidis JPA. Use of E-values for addressing confounding in observational studies – an empirical assessment of the literature. Int J Epidemiol 2020;49:1482-94. doi: 10.1093/ije/dyz261