n1

POLITECNICO DI MILANOFacoltà di Ingegneria dei Sistemi

Corso di Laurea Magistrale inIngegneria Gestionale

Un modello matematico per la predizionedel costo di cura dello scompenso

cardiaco sulla base dei dati contenuti neidatabase amministrativi della Regione

Lombardia

Relatore:Prof.ssa Cristina Masella

Co-Relatore: Ing. Giulia Garavaglia

Tesi di Laurea di:Lucas Bertolotti

Matr. 783061

Anno Accademico 2013-2014

n2

Ringrazio l'Ing. Giulia Garavaglia per la guida fornitanell'impostazione del lavoro, delle discussioni sui risultati dell'analisi

e della lettura delle bozze

n3

Indice1. Introduzione

1.1 Obiettivi del lavoro1.2 Lo scompenso cardiaco

2. Il Sistema sanitario2.1 Servizio Sanitario Nazionale2.2 Database amministrativi sanitari2.3 Sistema di Rimborso DRG

3. Descrizione del progetto HF DATA3.1 Introduzione (obiettivi e unità operative)3.2 Criteri di estrazione3.3 Database degli eventi

4. Analisi della letteratura4.1 Ricerca Strutturata4.2 Componenti del costo4.3 Dati mancanti e censored4.4 Modelli matematici di previsione4.4.1 Estrapolazione4.4.2 Modelli di regressione4.4.3 Modelli joint4.4.4 Statistica non parametrica4.5 Articoli di riferimento

5. Analisi dei dati5.1 Analisi esplorativa5.1.1 Metodologia per l'analisi esplorativa5.1.2 Risultati dell'analisi esplorativa5.2 Costruzione del modello5.2.1 Regressione con responso modellizzato tramite distribuzione gamma5.2.2 Regressione con responso modellizzato tramite distribuzione NormaleInversa Gaussiana (NIG)5.2.3 Regressione (NIG) con insieme di training mobile5.2.4 Calcolo dei p-value5.2.5 Analisi di rischio5.3 Analisi di segmentazione

6. Conclusione e svillupi futuri7. Appendice8. Bibliografia

n4

Introduzione

1.1 Obbiettivi del lavoro

Lo scompenso cardiaco è una sindrome eterogenea che porta spesso alle persone chene sono affette limitazioni funzionali e spesso da prognosi infausta e costituisce unadelle cause di ospedalizzazione più frequenti. Il miglioramento delle modalità di curadella patologia potrebbero avere un impatto sugli outcome del paziente scompensato,sia in termini di mortalità che in termini di riospedalizzazioni. Pertanto comprenderel'evoluzione epidemiologica dello scompenso cardiaco, quella dei profili di cura e delcorrelato consumo di risorse economiche risulta molto importante per i policy makersanitari. Lo studio dello scompenso cardiaco si è spesso basata sulle informazioni raccolte inregistri e database all'interno di studi clinici che hanno portato a enormi passi avantinella definizione, ad esempio, dei trattamenti farmacologici più efficaci e sullascoperta di quali siano i fattori prognostici dello scompenso. Tuttavia, come spessoaccade, la popolazione arruolata nei trial clinici non può essere consideratarappresentativa di quella reale, che spesso, dopo la dimissione ospedaliera, vieneseguita dal medico di medicina generale senza il supporto di un centro specializzatonella cura dello scompenso. I database amministrativi sanitari rappresentano le unichefonti di dati raccolti sugli eventi di pazienti non selezionati per studi clinici. Il progetto HF DATA, finanziato dal Ministero della Salute e Regione Lombardia,nell'ambito della quale il presente lavoro di tesi si svolge, nasce con l'obiettivo diprodurre risultati che possano supportare il policy maker regionale e quello nazionalenella comprensione dell'evoluzione delle modalità di cura dello scompenso cardiaco edei risultati fino ad ora ottenuti ed, eventualmente, a implementare azioni correttivedelle attuali politiche sanitarie.Lo scopo di questo lavoro è sviluppare un modello matematico perl'analisi/predizione del costo di cura dello scompenso cardiaco sulla base dei daticontenuti nei database amministrativi della Regione Lombardia e fornire una solidabase di evidenza a supporto del processo decisionale del policy maker regionale.

1.2 Lo scompenso cardiaco

Nella medicina c'è una differenza concettuale tra diagnosi di una malattia edefinizione di una malattia:

• Una malattia è una condizione biologica di funzionamento del corpo umano incondizioni ritenute non salutari, che genera effetti ritenuti non desiderabili:tipicamente morte e/o sofferenza.

n5

• La diagnosi invece è un procedimento medico atto a identificare nel pazienteuna o più malattie conosciute dai professionisti sanitari.

Da questa prima considerazione segue la necessità di classificare le malattie e lerelative metodologie diagnostiche. La branca della medicina che si occupa diclassificare le malattie è la nosologia. Storicamente ci sono stati diversi sistemiproposti per classificare le malattie, tra cui[56]:

• il più vecchio è quello di John Graunt (1663) con il titolo Natural andPolitical Observations Made upon the Bills of Mortality

• nell' Ottocento William Farr propose un sistema di classificazione per finalitàstatistiche

• Jacques Bertillon propose un sistema internazionale di classificazione nel1893, che si diffuse con il nome International List of Causes of Death

• nella quinta conferenza internazionale per la revisione della International Listof Causes of Death (1938) fu riconosciuta la necessità di una classificazioneche fosse valida per piani assicurativi, ospedali, servizi medici militari eorganizzazioni simili

• nella sesta revisione (1946), l'Organizzazione mondiale della sanità (OMSitaliano/WHO inglese) ha ricevuto la proposta del Committee on Joint Causesof Death degli USA per unificare la classificazione per fini di mortalità conquella per fini di morbidità

• la nona revisione (1975) introdusse la possibilità di classificare le malattie inbase alla parte del corpo affetta

Per questo lavoro useremo la classificazione ICD9-CM, che è una modificaamericana della nona revisione, elaborata dal National Center for Health Statistics edal Center for Medicare and Medicaid Services[64]. In particolare ci focalizzeremo sui pazienti che sono stati diagnosticati con la malattiache viene denominata scompenso cardiaco, di cui segue la definizione medica:

“Guasto della pompa cardiaca, in condizioni di carico fisiologiche, nell'imprimereuna sufficiente energia idraulica per mantenere la circolazione fisiologica”[1][57]

n6

Il Sistema Sanitario

2.1 Servizio Sanitario Nazionale

In Italia, il Servizio Sanitario Nazionale (SSN) è responsabile della “tutela dellasalute come diritto fondamentale dell'individuo e interesse della collettività”[51]. Ladefinizione del Piano Sanitario Nazionale viene fatta contestualmenteall'individuazione delle risorse finanziare destinate al SSN nel Documento diprogrammazione economico-finanziaria. Le Regioni contribuiscono alla definizionedel piano sanitario nazionale con proposte, la relazione annuale sullo stato diattuazione del Piano Sanitario Regionale, i risultati di gestione e la spesa prevista perl'anno successivo.Sono posti a carico del SSN le tipologie di assistenza, i servizi e le prestazionisanitarie che presentano, per specifiche condizioni cliniche o di rischio, evidenzescientifiche di un significativo beneficio in termini di salute (efficacia), a livelloindividuale o collettivo, a fronte delle risorse impiegate (efficienza). Le prestazioniinnovative per le quali non sono disponibili evidenze scientifiche di efficacia possonoessere erogate solo con programmi di sperimentazione autorizzati dal Ministero dellasalute.Il Piano sanitario nazionale indica:

a) le aree prioritarie di intervento, anche ai fini di una progressiva riduzione dellediseguaglianze sociali e territoriali nei confronti della Salute; b) i livelli essenziali di assistenza sanitaria da assicurare per il triennio di validitàdel Piano; c) la quota capitaria di finanziamento per ciascun anno di validità del Piano e lasua disaggregazione per livelli di assistenza;d) gli indirizzi finalizzati a orientare il Servizio sanitario nazionale verso ilmiglioramento continuo della qualità dell’assistenza, anche attraverso larealizzazione di progetti di interesse sovra-regionale;e) i progetti-obiettivo, da realizzare anche mediante l'integrazione funzionale eoperativa dei servizi sanitari e dei servizi socio-assistenziali degli enti locali; f) le finalità generali e i settori principali della ricerca biomedica e sanitaria,prevedendo altresì il relativo programma di ricerca; g) le esigenze relative alla formazione di base e gli indirizzi relativi allaformazione continua del personale, nonché al fabbisogno e alla valorizzazione dellerisorse umane; h) le linee guida e i relativi percorsi diagnostico-terapeutici allo scopo difavorire, all’interno di ciascuna struttura sanitaria, lo sviluppo di modalitàsistematiche di revisione e valutazione della pratica clinica e assistenziale e diassicurare l’applicazione dei livelli essenziali di assistenza;

n7

i) i criteri e gli indicatori per la verifica dei livelli di assistenza assicurati inrapporto a quelli previsti.

La Relazione sullo stato sanitario del Paese, predisposta annualmente dal Ministro della Salute:

a) illustra le condizioni di salute della popolazione presente sul territorio nazionale;b) descrive le risorse impiegate e le attività svolte dal Servizio sanitarionazionale;c) espone i risultati conseguiti rispetto agli obiettivi fissati dal Piano sanitarionazionale;d) riferisce sui risultati conseguiti dalle regioni in riferimento all'attuazione deipiani sanitari regionali;e) fornisce indicazioni per l’elaborazione delle politiche sanitarie e laprogrammazione degli interventi.

Le Regioni si occupano della determinazione dei principi sull'organizzazione deiservizi e sull'attività destinata alla tutela della salute e dei criteri di finanziamentodelle unità sanitarie locali e delle aziende ospedaliere, le attività di indirizzo tecnico,promozione e supporto nei confronti delle suddette unità sanitarie locali ed aziende,anche in relazione al controllo di gestione e alla valutazione della qualità delleprestazioni sanitarie. In particolare in base alla legge 724 del 23/10/1994, le Regionipossono emettere prestiti obbligazionari esclusivamente se destinati al finanziamentodegli investimenti, cioè non viene permessa dalla legge questa pratica per finanziarespese di parte corrente.La Regione disciplina anche:

a) l’articolazione del territorio regionale in unità sanitarie localib) i principi e criteri per l’adozione dell’atto aziendale c) la definizione dei criteri per l’articolazione delle unità sanitarie locali indistretti, tenendo conto delle peculiarità delle zone montane e a bassa densità dipopolazione;d) il finanziamento delle unità sanitarie locali, sulla base di una quota capitariacorretta in relazione alle caratteristiche della popolazione residente e) le modalità di vigilanza e di controllo, da parte della regione medesima, sulleunità sanitarie locali, nonché di valutazione dei risultati delle stesse f) l’organizzazione e il funzionamento delle attivitàg) fermo restando il generale divieto di indebitamento, la possibilità per le unitàsanitarie locali di: h) anticipazione nella misura massima di un dodicesimo dell’ammontare annuo

n8

del valore dei ricavi, inclusi i trasferimenti, iscritti nel bilancio preventivoannuale; contrazione di mutui e accensione di altre forme di credito, di duratanon superiore a dieci anni, per il finanziamento di spese di investimento e previaautorizzazione regionale, fino a un ammontare complessivo delle relative rate,per capitale e interessi, non superiore al quindici per cento delle entrate propriecorrenti, a esclusione della quota di fondo sanitario nazionale di parte correnteattribuita alla regione; le modalità con cui le unità sanitarie locali e le aziendeospedaliere assicurano le prestazioni e i servizi contemplati dai livelli aggiuntividi assistenza finanziati dai comuni

Le unità sanitarie locali si costituiscono in aziende con personalità giuridica pubblicae autonomia imprenditoriale, la loro organizzazione e funzionamento sonodisciplinati con atto aziendale di diritto privato, nel rispetto dei principi e criteriprevisti da disposizioni regionali. Le aziende ospedaliere devono chiudere il propriobilancio in pareggio. L'eventuale avanzo di amministrazione è utilizzato per gliinvestimenti in conto capitale, per oneri di parte corrente e per eventuali forme diincentivazione al personale da definire in sede di contrattazione. Il verificarsi di nongiustificati disavanzi di gestione o la perdita delle caratteristiche strutturali e diattività prescritte comportano rispettivamente la sospensione degli organi direttivi daparte della regione e la revoca dell'autonomia aziendale. I contratti di fornitura dibeni e servizi, il cui valore sia inferiore a quello stabilito dalla normativa comunitariain materia, sono appaltati o contrattati direttamente secondo le norme di dirittoprivato indicate nell’atto aziendale.Gli ospedali che non siano costituiti in aziende ospedaliere conservano la natura dipresidi dell'unità sanitaria locale, a cui è attribuita autonomia economico-finanziariacon contabilità separata all'interno del bilancio dell'unità sanitaria locale, con le stessedisposizioni previste per le aziende ospedaliere. Nelle unità sanitarie locali nelle qualisono presenti più ospedali, questi possono essere accorpati ai fini funzionali. Le unitàsanitarie locali e le aziende ospedaliere hanno disponibilità del patrimonio secondo ilregime della proprietà privata.Le regioni emanano norme per la gestione economico finanziaria e patrimoniale delleunità sanitarie locali e delle aziende ospedaliere:

a) la tenuta del libro delle deliberazioni del direttore generale;b) l’adozione del bilancio economico pluriennale di previsione nonché delbilancio preventivo economico annuale relativo all’esercizio successivo;c) la destinazione dell’eventuale avanzo e le modalità di copertura deglieventuali disavanzi di esercizio;d) la tenuta di una contabilità analitica per centri di costo e responsabilità checonsenta analisi comparative dei costi, dei rendimenti e dei risultati;e) l’obbligo delle unità sanitarie locali e delle aziende ospedaliere di rendere

n9

pubblici, annualmente, i risultati delle proprie analisi dei costi, dei rendimenti edei risultati per centri di costo e responsabilità; f) il piano di valorizzazione del patrimonio immobiliare anche attraversoeventuali dismissioni e conferimenti.

2.2 Database amministrativi sanitari

I database amministrativi sanitari sono basi di dati alimentate dai flussi informativigenerati dai provider di prestazioni sanitarie ai fini del rimborso e inviati all'ente chesi occupa del finanziamento dell'assistenza sanitaria. In Italia, dove il ServizioSanitario è pubblico, i provider inviano alle Regioni i dati necessari al finanziamentodelle prestazioni erogate.In Lombardia ogni provider del sistema sociosanitario che eroga una prestazione inregime di accreditamento con il Servizio Sanitario Regionale, per poter ottenere ilrimborso della prestazione, è tenuto a inviare un flusso informativo standard aRegione Lombardia. Per ciascun tipo di servizio esiste un dataset specifico. Ilcomplesso dei database alimentati dai diversi flussi informativi relativi ai diversiservizi forma il Datawarehouse Sanitario di Regione Lombardia. Ogni record cheviene inviato dal provider è associato a un paziente, che può essere identificatoattraverso il codice fiscale. Pertanto è possibile, tramite il codice fiscale, individuaretutti gli accessi effettuati dall'assistito al Servizio Sanitario Regionale consultando idiversi flussi.I database amministrativi sanitari possono essere usati per diversi scopi di ricerca,che contribuiscono al buon funzionamento del sistema sanitario nazionale[116][101]:

1. studi di farmaco-utilizzazione, distinguibili in due classi: statistiche sull’usodei farmaci e gli studi sull’appropriatezza d’uso di specifici gruppi di farmaci.Le statistiche sull’uso dei farmaci permettono di elaborare un rank dei farmacipiù prescritti e venduti, trend dei consumi di farmaci in funzione del tempo estudiare le differenze tra le prescrizioni in base alle aree geografiche. Gli studidi appropriatezza permettono un'analisi più accurata, come studiare l'utilizzo difarmaci in base all'età e al sesso.

2. Studi sulle ricette mediche3. Sviluppo di sistemi di sorveglianza e farmacosorveglianza: consiste nella

capacità di monitorare il rapporto benificio-rischio dei farmaci venduti, questoè particolare aiuto nel soddisfare l'esigenza di efficacia basata su evidenzescientifiche richiesta dalla legge

4. epidemiologia, tra cui farmacoepidemiologia5. analisi economiche, tra cui quelle strumentali al controllo di gestione regionale6. analisi di costo efficacia di trattamenti7. analisi sui percorsi di cura

n10

Se da un lato la disponibilità di questi database rende possibile studiare l'efficacia el'efficienza del sistema sanitario, dall'altro è importante tenere conto anche dellelimitazioni dell'utilizzo di questi strumenti, che saranno descritte in seguito.

A differenza di quanto accade per la costruzione dei database degli studi clinici, nelcaso dei database amministrativi sanitari il processi di acquisizione dei dati non èsotto il controllo diretto del ricercatore, ma il record viene generato solo se un eventoviene generato in concomitanza all'accesso al sistema sanitario da parte del cittadino(ricovero, accesso ambulatoriale, prescrizione di un farmaco, ecc.)Alcune fonti di possibili di errori nel database sono[101]:

• codifica dei medicinali (poco probabile)• assenza di diagnosi secondarie (comune)• cambiamenti nel hardware, software o nei metodi di codifica • fusione di diversi piani di salute possono portare al raddoppiamento dei numeri

di identificazione del paziente (nel caso di sistemi sanitari di altre areegeografiche1)

La presenza di regolamentazioni può essere un ostacolo alla realizzazione dellostudio, per esempio, negli Stati Uniti il Health Insurance Portability andAccountability Act (HIPAA) ha messo forti vincoli sull'utilizzo di dati per scopi nondirettamente legati alla cura dei pazienti. Il motivo di ciò è l'attenzione legata allaprotezione della privacy del paziente.[101]

L'utilizzo dei database amministrativi è perciò condizionato da una serie di criticità.Nel seguito si riportano alcuni contributi derivanti dalla letteratura scientifica cheindicano in che modo questi problemi possono essere affrontati e superati. Lericerche proposte sono state selezionate anche in virtù del focus, poiché infattiriguardano patologie cardiache o cardiovascolari. Lusignan et al. 2007[92]descrivono una procedura chiamata Total data quality management (TQDM) cheviene applicata a database di pazienti affetti da corono-patia, la principale causa discompenso cardiaco, che consiste nella seguente procedura:

1) Fase di design: perfezionamento della domanda di ricerca, identificazione del dataset idoneo, governance dell'informazione e protezione dati, feedback delle fasisuccessive2) Inserimento dati3) Estrazione dei dati: nell'esempio citato si è usato MIQUEST (Morbidity

1In Lombardia il linkage tra le diverse basi di dati si fa usando il codice fiscale, che èunico

n11

Information Query and Export Syntax) per prelevare i dati4) Migrazione: necessita di un mezzo di trasferimento, per esempio un floppy-disk, edi una conversione del formato file che permetta di alimentare il repository di dati5) Integrazione: tabelle dati, dati su sottoinsiemi di pazienti e dati seriali devonoessere collegati6) Pulizia: correzione di dati sbagliati o mancanti7) Processamento: conversione dei dati in inglese8) Analisi: analizzare la qualità dei dati ottenuti (completezza, accuratezza,consistenza, conversione valuta)

Walraven et al. 2012[93] indicano come sia importante considerare il frame dicampionamento e i metodi di campionamento e considerare che molti di questi data-set vengono creati per scopi non esclusivamente legati alla ricerca, tutto ciò può inlinea di principio compromettere la completezza e l' accuratezza. In generale il data-set tende ad essere più completo quando i responsabili della suacompilazione hanno un incentivo nel fornire correttamente i dati, per esempio se imedici vengono pagati in un sistema “fee-for-service”, allora saranno più incentivatinella compilazione dei dati. Altri aspetti importanti sono identificare la popolazione che è stata esclusa dal data-set e fornire una descrizione di tutte le variabili usate nel data set. In base a ciò èpossibile valutare l'impatto dei dati mancati e la possibilità di distorsione sistematichenel calcolo delle statistiche.Specificare i quality check (logici o random) usati e ricerche effettuare usando lostesso data-set è un altro aspetto importante, la review effettuata dagli autoristabilisce che 76% degli studi che usano database amministrativi usano codicidiagnostici e codici procedurali per definire la coorte di pazienti, successivamenteviene attribuito un codice in base a una sistema di classificazione standard la cuiaccuratezza dipende dalla validità della diagnosi e dalla procedura di assegnazionedei codici:

• L'accuratezza della diagnostica potrebbe variare per esempio a seconda dellaspecializzazione del medico che effettua la diagnostica, oppure differire tragli ospedali per quanto riguarda le diagnostiche che necessitano di tecnologiemeno comuni.

• Per quanto riguarda la procedura di assegnazione del codice, il rischio dimisclassificazione e la presenza di bias dovrebbero essere calcolati in base aun'analisi statistica

Tre design di studio sono tipicamente usati:

1. Studi ecologici: comparano statistiche calcolate sulla base di differenti

n12

popolazioni2. Studi di re-astrazione: ripetono la procedura di assegnazione dei codici 3. Studi gold standard: comparazione del codice con un gold standard (variabili

cliniche standard, panello oppure un database più accurato)

Kuwabara et al. 2010[94] costruiscono una database amministrativo, incollaborazione con il ministero della salute giapponese, il quale contiene dati didimissione e dati claim in formato elettronico. Questo database è stato utilizzato dalministero per valutare la performance di 731 ospedali. I pazienti scompensati sonostati identificati usando la classificazione ICD10, sulla base dei giudizi dei clinici si èvalutato se si trattasse di un caso di malattia acuta (ospedalizzazione urgente e nonpianificata) o cronica (impossibilità di eseguire attività giornaliere), i casi di trapiantosono stati esclusi. I costi ospedalieri vengono calcolati in base a un sistema dipagamento “fee-for-service”, nello studio sono stati considerati costi di laboratorio,costi di strumenti, costi del personale medico e costi amministrativi. Per valutare lagravità della cronicità è stato usato il Charlson Comorbidity Index[95]. Alcuni limitidello studio sono il periodo di rilevazione delle informazioni per pazienti dimessi (6mesi nel 2006), assenza di dati ambulatoriali, assenza di dati sulle condizionifisiologiche di ingresso dei pazienti.

Chen et al. 2011[96] prendono in considerazione il database National HealthInsurance Research Database (NHIRD) di Taiwan. Il data set NHIRD contiene recordriguardanti visite mediche (sia con medicina occidentale che con medicina orientale),visite dentistiche, ospedalizzazione, farmaci somministrati, esami diagnostici dilaboratorio e codici procedurali. Per proteggere la privacy dei pazienti i codiciidentificativi sono unici al data-set. Una delle critiche mosse a questo data-set era lamancanza di note dei medici, report di laboratorio e immagini diagnostiche

Saczynski et al. 2012[97] conducono una review della letteratura dove analizzano glistudi epidemiologici in base a: sistema codifica, database usato, popolazione studiata,periodo temporale, distinzione tra prevalenza e incidenza, algoritmo usato peridentificare i casi di scompenso cardiaco, criterio di aggiudicazione, processo divalidazione e statistiche di validazione. Le codifiche usate erano ICD-8, ICD-9 eICD-10 e DRG. Sono stati inseriti nella review solo articoli con “κ di Cohen”2

superiore a 0.65.

Quach et al. 2010[98], in una review, analizzano 25 studi condotti tra il 1992 e 2008 ehanno rilevato che:

2 , dove Pr(a) è il livello di accordo tra i reviewer, Pr(e) è la probabilità di

concordanza casuale

n13

• 9 sono stati realizzati in Canada, 11 negli USA, uno nella Svezia, Portogallo,Regno Unito, Australia e in Danimarca rispettivamente.

• La dimensione del campione varia tra 126 e 58816 pazienti.• Più di 70 codici ICD sono stati usati per definire casi di scompenso tra cui:

ICD-8, ICD-9 e ICD-10.• 19 studi hanno usato la revisione di cartelle cliniche come benchmark e gli

altri 6 hanno usato sondaggi, registri o record ricodificati• Diagnostica: 7 studi hanno accettato le diagnostiche nelle cartelle, 4 hanno

usato casi documentati di scompenso, 8 hanno usato criteri diagnostici comeFramingham, Carlson, Boston, the National Health and Nutrition ExaminationSurvey (NHANES), the New York Heart Association (NYHA) or theEuropean Society of Cardiology (ESC).

• 21 studi hanno effettuato una validazione tramite dati di dimissioneospedaliera

• 1 studio ha integrato i dati con altre fonti come facility, phyisician claim efarmaci

• 4 studi hanno esaminato la differenza tra l'uso di codifiche principali versuscodifiche secondarie

Avendo quindi inquadrato il tema dei database in relazione al sistema sanitario e airisultati della letteratura, descriviamo più in dettaglio il Sistema Informativo Socio-Sanitario (SISS) di Regione Lombardia, il quale sarà fondamentale per l'analisi dellapopolazione di scompensati. Questo sistema informativo nasce come riconfigurazionedel vecchio Sistema Sanitario Lombardo in base alla Legge 31/97. Precedentementeil sistema era costituito da un insieme di applicazioni ad uso di Regione Lombardiache evolvevano in modo disgiunto sul territorio. Una pietra miliare fu l’introduzionedella piattaforma CRS-SISS (Carta Regionale dei servizi - Sistema Informativo Socio-Sanitario), questo sistema “racconta” la storia dei cittadini e non quella delle strutturedell’offerta (flussi tradizionali di rendicontazione), consentendo di effettuare inmaniera più veloce ed efficace attività come i controlli di appropriatezza degliinterventi (analisi dei percorsi) o la programmazione regionale, che può esseresupportata da analisi ex-post, fondate sui criteri di occupazione dell’offerta, ed analisiprevisionali ex-ante, fondate sulla unione di analisi epidemiologiche e percorsiassistenziali di riferimento.

Il sistema è stato concepito per raccogliere i dati generati durante i percorsi delcittadino all'interno del sistema sanitario e permetterne l'analisi in base a diversi puntidi vista, consentendo l’attuazione di politiche mirate da parte dell’Amministrazioneregionale. Ogni dato raccolto mediante il SISS durante lo svolgimento di tali percorsiha una doppia valenza e viene qualificato sia come dato clinico sia come dato

n14

amministrativo. In tal modo si vengono a “costituire” due grandi basi informative organizzate pergenerare le conoscenze necessarie all’esercizio della professione clinica e al governodel Sistema Socio-Sanitario Regionale, ovvero: il Fascicolo Sanitario Elettronico(FSE) e il data warehouse Sanità (DWH). Il Fascicolo sanitario elettronico (FSE) èuna raccolta di dati e informazioni sanitarie che costituiscono la storia clinica e disalute di una persona. La consultazione del FSE avviene in forma protetta e riservataattraverso l’utilizzo di credenziali personali. Il DWH consente di accedere a report erealizzare analisi sofisticate, inoltre la disponibilità in tempo reale su DWH delprescritto sia farmaceutico, sia specialistico ambulatoriale, è di supporto per ulteriorianalisi di trend di eventi sanitari, importanti per anticipare e correggere l’emergere dianomalie di sistema o per il processo di pianificazione e controllo. Esiste anche laBanca Dati Assistiti (BDA), la quale è una “vista” del DWH ove sono aggregati tutti iconsumi sanitari, strumento essenziale sia per analisi a supporto dellaprogrammazione della spesa sanitaria, sia per analisi tipicamente epidemiologichestatistiche, sia per valutare l'appropriatezza e la qualità delle prestazioni, pereffettuare analisi trasversali incentrate sulle diverse patologie e sui percorsi di cura. Il SISS è alimentato da tutti i flussi amministrativi generati dai provider accreditatidel Servizio Sanitario regionale, ad esempio il flusso delle Schede di DimissioneOspedaliera (SDO), il flusso delle prestazioni ambulatoriali specialistiche e il flussodella farmaceutica (FUR)[52].Lo studio oggetto di questa tesi si basa sull'impiego del database delle SDO.Nel seguito si riporta una descrizione della struttura del database e delle informazioniin esso contenute.La Scheda di Dimissione Ospedaliera (SDO) è stata istituita dal D.M. 28/12/1991quale “strumento ordinario per la raccolta delle informazioni relative ad ognipaziente dimesso dagli istituti di ricovero pubblici e privati in tutto il territorionazionale”.Lo stesso decreto stabilisce anche i contenuti della scheda, successivamente spiegatinel D.M. 26/07/1993 e poi aggiornati nel D.M. 27/10/2000:

a) la sezione prima, che contiene le informazioni anagrafiche di seguito riportate: 1) denominazione dell'ospedale di ricovero;2) numero della scheda: primi 2 caratteri indicano l'anno e i successivi 6 unanumerazione progressiva dentro l'anno3) cognome e nome del paziente; 4) sesso; 5) data di nascita; 6) comune di nascita;

n15

7) stato civile; 8) comune di residenza; 9) cittadinanza; 10) codice sanitario individuale; 11) regione di residenza; 12) azienda unità sanitaria locale di residenza; b) la sezione seconda, che contiene almeno le informazioni del seguente elenco, la cui numerazione riprende e prosegue la numerazione dell'elenco di cui alla precedente lettera a): 1) denominazione dell'ospedale di ricovero; 2) numero della scheda; 13) regime di ricovero: ricovero ordinario oppure ricovero diurno (<24h)14) data di ricovero; 15) unità operativa di ammissione: primi 2 caratteri indicano specialità clinica esuccessivi e 2 successivi indicano l'unita operativa16) onere della degenza: soggetto su cui ricade l'onere di rimborsare le spese relativeal ricovero17) provenienza del paziente:

• 1) paziente che acceda all'istituto di cura senza proposta di ricovero • formulata da un medico;• 2) paziente inviato all'istituto di cura con proposta del medico di base; • 3) ricovero precedentemente programmato dallo stesso istituto di • cura; • 4) paziente trasferito da un istituto di cura pubblico; • 5) paziente trasferito da un istituto di cura privato accreditato; • 6) paziente trasferito da istituto di cura privato non accreditato; • 7) paziente trasferito da altro tipo di attività di ricovero (acuti, • riabilitazione, lungodegenza) o da altro regime di ricovero (ricovero • diurno o ordinario) nello stesso istituto; • 9) altro.

18) tipo di ricovero: ricovero programmato, ricovero urgente, ricovero pertrattamento sanitario obbligatorio, ricovero programmato con preospedalizzazione 19) traumatismi o intossicazioni; 20) trasferimenti interni; 21) unità operativa di dimissione; 22) data di dimissione o morte; 23) modalità di dimissione:

n16

• 1) nel caso in cui il paziente sia deceduto; • 2) dimissione ordinaria al domicilio del paziente; • 3) dimissione ordinaria presso una residenza sanitaria assistenziale (RSA); • 4) dimissione al domicilio del paziente con attivazione di ospedalizzazione

domiciliare; • 5) dimissione volontaria (da utilizzare anche nei casi in cui il paziente in ciclo

di trattamento diurno non si sia ripresentato durante il ciclo programmato);• 6) trasferimento ad un altro istituto di ricovero e cura, pubblico o privato, per

acuti; • 7) trasferimento ad altro regime di ricovero o ad altro tipo di attività di ricovero

nell'ambito dello stesso istituto; • 8) trasferimento ad un istituto pubblico o privato di riabilitazione; • 9) dimissione ordinaria con attivazione di assistenza domiciliare integrata.

24) riscontro autoptico; 25) motivo del ricovero in regime diurno; 26) numero di giornate di presenza in ricovero diurno; 27) peso alla nascita; 28) diagnosi principale di dimissione: codifica ICD9-CM29) diagnosi secondarie; 30) intervento chirurgico principale o parto; 31) altri interventi terapeutiche.

Vale anche il seguente caveat [5]:

“Il confronto tra dati relativi a diversi anni può risentire del diverso grado dicompletezza ottenuto nelle Regioni oppure di modifiche organizzative intervenute odi cambiamenti nelle definizioni o nelle codifiche adottate.Da alcune elaborazioni specifiche sono esclusi i casi in cui si rilevano errori nelleinformazioni analizzate, che non ne hanno consentito l'utilizzazione. Nonostanteun’intensa attività di controllo e di miglioramento della qualità dei dati, in alcunicasi non è stato possibile eliminare alcune rare incongruenze tra informazionilogicamente correlate.”

2.3 Sistema di Rimborso DRG

Il sistema di rimborso e quindi di finanziamento degli ospedali accreditati inLombardia è regolato dal sistema dei Diagnosis Related Group (DRG). Il sistema deiDRG classifica alla dimissione i pazienti di un ospedale in classi omogenee per

n17

diagnosi e consumo di risorse., più precisamente si può quindi attribuire ad ogni SDOun codice numerico DRG (Diagnosis Related Group) e ad ogni DRG è possibileassociare un valore monetario. Un' ospedale provvede ai suoi pazienti una serie diservizi: se i pazienti possono essere aggruppati tramite attributi clinici e somiglianzenel processo di erogazione dei servizi, allora abbiamo a che fare con un problema dicase mix, ovvero un mix di casi medici. Per quanto riguarda la Lombardia i DRGvengono attribuiti con il software Grouper della 3M Company, e i valori monetari(tariffe dei DRG) vengono fissati a livello regionale e poi pubblicati nella GazzettaUfficiale.

Concomitante allo sviluppo dei DRG, c'è stato quello delle Major DiagnosticCategories (MDC) [55] che consisteva nel suddividere tutte le possibili diagnosiprincipali in 23 categorie mutuamente esclusive. Le diagnosi del MDC corrispondonoa diverse parti del corpo umano e per motivi di interpretabilità da parte del personalemedico fu stabilito che uno stesso DRG non poteva contenere dentro di se diversicodici di MDC.

L'idea dei DRG risale principalmente agli studi fatti nell'Università di Yale dairicercatori Robert Fetter e John D.Thompson e poi successivamente introdotto dalgoverno americano nel 1983 nel sistema Medicare, sostanzialmente un piano diassicurazione per persone con età superiore a 65 anni. [6] Questo metodo dei DRG hasubito una lunga evoluzione nel lungo degli anni ed è anche stato adattato alle diverserealtà ospedaliere.

Possiamo comunque capire l'intuizione dell' approccio dall'articolo più citato di Fetterdove esamina un database di record ospedalieri del New Jersey [7]. Sia il responsoindicato con Y e la covariata con X, la quale assume N valori distinti. Il sottoinsiemedelle osservazioni con valore Xi (dove 1 ≤ i ≤ N) viene chiamato categoria. Se M i è ilnumero di osservazioni nella categoria i, allora la somma totale dei quadrati TSSQ(total sum of squares) è data da:

Dove la media del responso sul dataset è:

Vogliamo quindi partizionare il dataset in base ai valori {1,2,....,N} assunti dallacovariata in G insiemi disgiunti Rk :

n18

Questi insiemi Fetter et al. chiamano (impropriamente nel senso matematico) gruppi,si può quindi calcolare il TWGSSQ(G) (Total Within Groups Sum of Squares):

Questo è la somma dei SSQ dei responsi effettuata per categoria e per gruppo.L'algoritmo di classificazione del software AUTOGRP [8] procedeva quindi aminimizzare TWGSSQ per una covariata assegnata: in questo modo ogniosservazione appartiene ad un solo gruppo e il valore del responso è simile alla mediadel gruppo. I ricercatori hanno rimosso alcune osservazioni: quelle con pazientimorti, errori di codifica, dati mancanti e osservazioni con un periodo di degenzaparticolarmente lungo (outlier dal loro punto di vista). Le covariate date in inputall'algoritmo erano le diagnosi, procedure chirurgiche, età, sesso e servizi cliniciricevuti. Il responso scelto è stato la durata della degenza, usato come proxy per ilcosto.É stato importante complementare il modo di procedere iterativo dell'algoritmo con igiudizi dei medici. Lo schema finale che risultò doveva e deve ancora soddisfarealcuni requisiti:

• Deve essere interpretabile dai medici sia dal punto di vista della diagnosi, siadal punto di vista del management dei pazienti

• I gruppi devono essere definiti usando variabili tipicamente disponibili pressogli ospedali

• Il numero finale di gruppi deve essere sulla scala delle centinaia• I gruppi devono contenere pazienti con responsi simili• La definizione dei gruppi deve essere comparabile attraverso diversi schemi di

codifica

Questa analisi a suo tempo risultò in 383 DRG.Esistono oggi nel mondo diverse tipologie di DRG, tra cui[61]:

• DRG Medicare• DRG Refined (R-DRGs)• Severity DRG (S-DRGs)

n19

• All patients DRGs (AP-DRGs)• All Patient Refined DRGs (APR-DRGs)

Alcuni paesi, oltre a recepire i DRG americani, hanno elaborato DRG nazionali, peresempio in Germania esistono i G-DRG. In Italia sono usati allo stato attuale solo iDRG Medicare[62], i quali vengono elaborati e descritti dal Center for Medicare &Medicaid Services (CMS) nei DRG Definitions Manual.[53] L'emissione dei DRG, oltre ad essere fondata su analisi di tipo statistico e medico, èanche regolata dalla legge, in particolare il D.M. 15/04/1994 individua (con una seriedi allegati) le prestazioni da remunerare, i criteri di determinazione delle tariffe, leattività di controllo da esercitare. Per quanto riguarda le tariffe:

1. Le regioni e le province autonome determinano le tariffe delle prestazioni daapplicare nel proprio ambito territoriale. Le tariffe sono fissate sulla base delcosto standard di produzione e dei costi generali, in quota percentuale rispetto aicosti standard di produzione.2. Il costo standard di produzione per prestazione è calcolato in via preventivadalle regioni e dalle province autonome, sulla base dei costi rilevati presso uncampione di soggetti erogatori, pubblici e privati, operanti rispettivamentenell'ambito del servizio sanitario nazionale del territorio regionale e provinciale,preventivamente individuato secondo criteri di efficienza ed efficacia. Tale costofa riferimento alla composizione ed alla qualità di fattori produttivi utilizzati per laproduzione della prestazione, valorizzati sulla base dei prezzi unitari medi di acquistoriferiti all'ultimo anno e delle relative eventuali variazioni attese in ragione del tassodi inflazione programmato. Le componenti di costo da considerare per il calcolo delcosto standard di produzione della prestazione sono le seguenti:

• il costo del personale direttamente impiegato• il costo dei materiali consumati • il costo delle apparecchiature utilizzate (manutenzione, ammortamento),

proporzionato ad un tasso di utilizzo predeterminato a livello regionale;• i costi generali della unità produttiva della prestazione, ossia il costo dei fattori

di produzione attribuiti alla unità produttiva ma non direttamente utilizzatinella produzione della singola prestazione, distribuiti proporzionalmente tratutte le prestazioni da questa prodotte.

3. Per la determinazione delle tariffe, il costo standard di produzione perprestazione viene incrementato di una quota percentuale corrispondente alvalore medio rilevato del rapporto tra tali costi generali di struttura e la sommadei costi precedenti4. Nel determinare le tariffe per le prestazioni rese dai soggetti erogatori per i

n20

quali siano documentati lo svolgimento di attività istituzionali di ricerca edidattica e/o la presenza di servizi obbligatori individuati dalla programmazionenazionale e regionale, le regioni e le province autonome possono incrementare ilcosto standard di produzione di una quota percentuale proporzionale ai costiattribuiti a tali attività.5. In fase di prima applicazione del decreto, ai fini del calcolo del costo standarddi produzione delle prestazioni le regioni e le province autonome possonolimitare la rilevazione dei costi ad un campione di prestazioni. Al fine di stimare icosti standard di produzione delle prestazioni non comprese nel campione le regioni ele province autonome possono utilizzare il sistema di pesi relativi, con gliaggiustamenti eventualmente conseguenti alle rilevazioni campionarie dei costi diproduzione per prestazione. Le regioni e le province autonome provvedonoannualmente a verificare, ed eventualmente a rettificare, il sistema di pesi relativisulla base dei costi di produzione.6. Le regioni e le province autonome, con periodicità almeno triennale,provvedono all'aggiornamento delle tariffe, tenendo conto delle innovazionitecnologiche e delle variazioni dei costi delle prestazioni rilevate.7. Tutti i soggetti erogatori, pubblici e privati, che operano nell'ambito delServizio sanitario nazionale sono tenuti a trasmettere alle rispettive regioni eprovince autonome le necessarie informazioni sui propri costi di produzione,nonché ad attestarne la veridicità e la corrispondenza alle proprie scritture contabili,secondo le modalità e la periodicità definite dalle regioni e province autonome diappartenenza.

Le regioni e le province autonome vigilano sulla corretta applicazione da parte delleunità sanitarie locali del sistema di remunerazione mediante tariffe definite ai sensidel suddetto decreto, avvalendosi anche delle commissioni regionali per lapromozione della qualità delle attività sanitarie. Al fine di consentire l'acquisizionedelle informazioni necessarie alla programmazione sanitaria nazionale, le regioni e leprovince autonome inviano al Ministero della Sanità i provvedimenti regionali eprovinciali di determinazione delle tariffe delle prestazioni, corredati dei relativi datidi riferimento sui costi, entro sessanta giorni dalla loro approvazione.I Medicare DRG hanno subito una serie di revisioni nel tempo, in Lombardia inparticolare sono state usate le seguenti versioni:

• La versione 14.0 dei CMS DRG è stata usata fino al dicembre 20023

• La versione 19.0 dei CMS DRG è stata introdotta il 1° gennaio 2003[58][60]• La versione 24.0 dei CMS DRG è utilizzata in fase di test durante il 2008, e poi

introdotta ufficialmente il 1° gennaio 2010[59][60]

3 In Italia sono state usate solo le versioni 10.0, 14.0, 19.0 e 24.0 (vedi http://www.drg.it/ oppure Nonis e Rosatti[63])

n21

Descrizione del progetto HF DATA

3.1 Introduzione (obiettivi e unità operative)

Questo lavoro si inserisce all'interno del progetto HF DATA[50], un progetto diricerca finanziato dal Ministero della Salute dal 2009 e da Regione Lombardia, il cuiresponsabile di progetto è l'Azienda Ospedaliera Niguarda Cà Granda, con i seguentiobiettivi:

1. Descrivere l'epidemiologia e la storia naturale dei pazienti affetti da scompensocardiaco

2. Descrivere il consumo di risorse per la cura dello scompenso cardiaco3. Valutare le analogie e le differenze nella cura dei pazienti in rapporto all'area,

all'età, al sesso, e ad altre variabili4. Valutare la prescrizione di farmaci e i fattori che la condizionano

Lo studio si basa sui dati contenuti nei database amministrativi di RegioneLombardia, in particolare i pazienti ospedalizzati almeno una volta per scompensocardiaco e residenti in Lombardia tra il 2005 e 2012. L'attività di progetto vieneripartita in 8 unità operative:

• UO1 Regione Lombardia - DG Sanità – Unità organizzativa programmazione esviluppo piani: fornisce l'expertise epidemiologico e di metodologia di analisidei database istituzionali, sulla base dell'esperienza di progetti precedenti sonostate proposte diverse metodologie di filtro per l'estrazione di dati

• UO2 Politecnico di Milano – MOX Laboratorio di Modellistica e CalcoloScientifico-Dipartimento di Matematica: si occupa dello studio dei modellistatistici multi-stadio per l'analisi della progressione di una malattia, sia dalpunto di vista della fruizione delle risorse cliniche, sia dal punto di vista delpeggioramento effettivo della malattia

• UO3 Politecnico di Milano – DIG Dipartimento di Ingegneria Gestionale: sioccupa della definizione dei criteri di estrazione del campione di malati sullabase del flusso SDO, selezione di altre fonti informative, definizione dellametodologia per il trattamento dei dati finalizzata alla strutturazione di unmodello di malattia

• UO4 CEFRIEL -Digital Enterprise & Information Systems: contribuisce allagestione del progetto con la realizzazione e formalizzazione degli incontri, ildisegno del sito di progetto http://hfdata.cefriel.it/ e la sua messa in linea,analisi di fattibilità delle tecniche di text mining per analizzare le letter didimissione dei pazienti

n22

• UO5 Lombardia Informatica S.p.A. -Funzione Service Management Regionale– Sistemi clinici: fornisce attraverso una piattaforma di DataWareHouse unacesso integrato ai dati amministrativi tradizionali e ai dati SISS

• UO6 IRCSS Fondazione S. Maugeri Lumezzane, UO7 AO Niguarda Ca'Granda – SC Cardiologia2, UO8 AO San Carlo Milano – SC Cardiologia:contribuiscono nella definizione dei dati da estrarre dai database, modalità diassociazione dei dati per valutare la gravità del profilo di malattia e distingueretra le diagnosi e le procedure eseguite in ricovero, quelle che costruisconofattori prognostici favorevoli e sfavorevoli

3.2 Criteri di estrazione

L'estrazione dei dati avviene in tre fasi:

1. Estrazione dati delle SDO relative ai ricoveri effettuati in degenza ordinaria(2000-2012), limitatamente alle MDC 1,4,5 e 11, integrati con dati dagliaccessi di pronto soccorso (2005-2012) e con i dati anagrafici dell'eventualedata di decesso.

2. Estrazione di tutti i ricoveri dal 2000-2012 e alcune prestazioni ambulatoriali,dei pazienti individuati nel “gruppo SDO scompenso”

3. Vengono effettuate richieste mirate di dati relativi alla protesica e alleprescrizioni farmaceutiche

I dati estratti nella prima fase sono stati sottoposti a un accurato controllo di qualità esono stati manipolati per facilitare le analisi effettuate dalle diverse unità operative.I pazienti che entrano nello studio a partire dal 2005 non hanno eventi di scompensonei 5 anni precedenti e possono perciò considerarsi come nuovi casi incidenti discompenso. Sono stati identificati due criteri per l'estrazione dei casi di scompenso dal flussoSDO: Inpatient Quality Indicators (IQI), proposti dalla Agency for HealthcareResearch & Quality (AHRQ), e Hierarchical Condition Category (HCC) proposti dalCenter for Medicare and Medicaid Services (CMS). Più in particolare è stato usato l'IQI-16 (versione 4.3 agosto 2011) e l'HCC 80 (versione 12), applicati alle MDC 1, 4, 5 e 11, con l'obbiettivo di rendere più omogenea la casistica degli scompensati.Questi due criteri sono stati individuati dalla UO3 DIG attraverso una review diletteratura usando come chiave di ricerca: “administrative data” OR “claims data”in All Fields AND “Heart Failure” in Mesh. Il criterio HCC è costruito a partiredai profili di rischio dei pazienti: i codici ICD9-CM vengono raggruppati in 805gruppi diagnostici (DXG) in modo univoco ed esaustivo, i DXG vengono poiraggruppati in 189 categorie di condizioni (CC) e infine viene imposta una strutturagerarchica alle categorie di condizioni arrivando cosi alle HCC.

n23

I due criteri non sono completamente sovrapponibili, per i nostri scopi consideriamol'intersezione AHRQ∩HCC e le differenze AHRQ\HCC=AHRQ\ eHCC\AHRQ=HCC\:

Partendo da questa logica, gli eventi di ricovero relativi a ogni paziente sono staticlassificati sulla base dei codici di diagnosi contenute nella SDO del primo ricoveroin 4 gruppi sulla base delle indicazioni dei clinici coinvolti nel progetto:

1. Gruppo 1: criterio HCC in diagnosi primaria e criterio HCC in diagnosisecondaria

2. Gruppo 2: criterio HCC\ come diagnosi primaria e HCC\ come diagnosisecondaria

3. Gruppo 3: nessuno evento di scompenso nella diagnosi primaria, criterioAHRQ∩HCC nella diagnosi secondaria

4. Gruppo 4: criterio AHRQ\ come diagnosi primaria o AHRQ\ come diagnosisecondaria, nel caso di nessuno evento di scompenso nella diagnosi primaria

3.3 Database degli eventi

Le analisi sono state effettuate sul database degli eventi di ricovero per scompenso

n24

cardiaco. Un evento può essere costituito da un ricovero o da più ricoveri consecutivi.Due ricoveri sono stati considerati come consecutivi se il secondo inizia nello stessogiorno o nel giorno successivo alla fine del primo, un insieme di ricoveri consecutiviindividua una unità di osservazione o evento. Una variabile di flag, è una variabilebinaria che assume valore 1 nel caso sia soddisfatta la condizione di flag.Il database degli eventi, alimentato dai campi del flusso SDO, è costituito dalleseguenti variabili4:

V1: COD_REG codice identificativo del paziente

V2: evento_id numero che identifica quell’evento all’interno di tutti glieventi del paziente, non necessariamente parte da 1 oincrementa sempre di 1

V3: età età del paziente a inizio evento

V4: data_inizio_ev data inizio del primo ricovero dell’evento

V5: data_rif_ev data di fine del primo ricovero dell’evento

V6: data_dec_tronc data decesso del paziente dalla Nuova Anagrafe Regionale,non da SDO

V7: data_mpr_Min data della main procedure per il calcolo del DRG

V8: RES_LOMB_Min da usare insieme a “RES_LOMB_Max” (V17):• se res_lomb_min = 0 e res_lomb_max = 1, allora il

paziente è stato sia residente che non residente durantel’evento

• se res_lomb_min = 1 e res_lomb_max = 1, allora il paziente è sempre stato residente durante l’evento

V9: flag_riab_Min da usare insieme a “flag_riab_Max” (V18):• se res_ riab _Min = 0 e res_ riab _Max = 1, allora

nell’evento ci sono stati ricoveri di riabilitazione e anche non di riabilitazione

• se res_ riab _Min = 1 e res_ riab _Max = 1, allora l’evento è fatto solo da ricoveri di riabilitazione

• se res_ riab _Min = 0 e res_ riab _Max = 0, allora nell’evento non ci sono ricoveri di riabilitazione

V10: data_fine_ev data di fine dell’evento

V11: dec_intra flag di decesso intraospedaliero

V12: AHRQ_1dx_Max flag che indica se uno dei ricoveri dell’evento ha un codiceAHRQ in prima diagnosi

V13: AHRQ_nddx_Max flag che indica se uno dei ricoveri dell’evento ha un codiceAHRQ in una delle diagnosi secondarie

V14: HCC80_1dx_Max flag che indica se uno dei ricoveri dell’evento ha un codiceHCC in prima diagnosi

V15: HCC80_nddx_Max flag che indica se uno dei ricoveri dell’evento ha un codiceHCC in una delle diagnosi secondarie

4 Per sinteticità e scopi di programmazione si indicano le variabili del data-set con V1, V2,..., VN

n25

V16: n_ric_inev numero di ricoveri che compongono l’evento

V17: RES_LOMB_Max vedi “RES_LOMB_Min” (V8)

V18: flag_riab_Max vedi “flag_riab_Min” (V9)

V19: flag_ti_Max flag che indica se in uno dei ricoveri dell’evento il paziente èpassato in terapia intensiva, una terapia necessaria per ilsupporto delle funzioni vitali del paziente dopo un interventochirurgico

V20: flag_cardiochir_Max flag che indica se in uno dei ricoveri dell’evento il paziente èpassato in un reparto di cardiochirurgia

V21: flag_ICD_Max flag che indica se in uno dei ricoveri dell’evento il paziente hasubito un intervento per impiantare un defibrillatore, unatecnologia che consiste nell'uso di una corrente elettrica nelcuore del paziente per combattere l'aritmia

V22: flag_drgchir_Max flag che indica se uno dei ricoveri dell’evento ha un DRGchirurgico

V23: AHRQ_escl_1dx_Max flag che indica se uno dei ricoveri dell’evento ha un codiceesclusivo AHRQ in prima diagnosi

V24: HCC80_escl_1dx_Max flag che indica se uno dei ricoveri dell’evento ha un codiceesclusivo HCC in prima diagnosi

V25: AHRQ_e_HCC80_1dx_Max flag che indica se uno dei ricoveri dell’evento ha un codicedell’intersezione tra AHRQ e HCC in prima diagnosi

V26: AHRQ_escl_nddx_Max flag che indica se uno dei ricoveri dell’evento ha un codiceesclusivo AHRQ in una delle diagnosi secondarie

V27: HCC80_escl_nddx_Max flag che indica se uno dei ricoveri dell’evento ha un codiceesclusivo HCC in una delle diagnosi secondarie

V28: AHRQ_e_HCC80_nddx_Max flag che indica se uno dei ricoveri dell’evento ha un codicedell’intersezione tra AHRQ e HCC in una delle diagnosisecondarie

V29: VAL_TOT_EURO_Sum somma delle valorizzazioni dei ricoveri che compongonol’evento

V30: VAL_PROTESI_Sum somma delle valorizzazioni delle protesi dei ricoveri checompongono l’evento

V31: VAL_TOT_EURO_Nmiss numero di ricoveri dell’evento per cui non è stata data lavalorizzazione

V32: VAL_PROTESI_NMiss variabile da non usare

V33: SESSO sesso del paziente al primo ricovero dell’evento

V34: ASL_RESIDENZA ASL residenza del paziente al primo ricovero dell’evento

V35: STRUTT_RIC_SUB_ID variabile da non usare

V36: DISCIP_ENT_ID disciplina del reparto medico di entrata del primo ricoverodell’evento

V37: oss variabile da non usare

V38: data_studio_out data di uscita dallo studio

n26

V39: desc_studio_out descrizione dell’uscita dalla studio, può assumere i seguentivalori:

• NON APPLICABILE: il paziente non è residente enon ha data di decesso, il paziente è residente ma nonha codice identificativo, la data di decesso è compilatama precedente al 2000

• DECEDUTO: il paziente ha una data di decesso• TRONCATO: il paziente è arrivato vivo al 31-12-2012• PERSO: il paziente non è deceduto, ma è uscito

dall’anagrafica prima del 31-12-2013

V40: anno_rif_ev data di fine del primo ricovero dell’evento; usata perassegnare gli anni agli eventi

V41: durata_ev differenza tra data_fine_ev e data_inizio_ev, cioè la durata diospedalizzazione; i valori negativi pari a -1 o -2 sono statiposti a zero, gli altri valori negativi sono stati posti a missing

V42: DISCIP_USC_ID disciplina del reparto medico di uscita dell’ultimo ricoverodell’evento

V43: tempo_da_evprec variabile da non usare

V44: metastatic_romano_Max flag che indica la presenza di cancro metastatico

V45: chf_romano_Ma flag che indica la presenza di scompenso cardiaco congestizio

V46: dementia_romano_Max flag che indica la presenza di demenza

V47: renal_elixhauser_Max flag che indica la presenza di insufficienza renale

V48: wtloss_elixhauser_Max flag che indica la presenza di perdita di peso

V49: hemiplegia_romano_Max flag che indica la presenza di emiplegia

V50: alcohol_elixhauser_Max flag che indica la presenza di abuso di alcool

V51: tumor_romano_Max flag che indica la presenza di tumore

V52: arrhythmia_elixhauser_Max flag che indica la presenza di aritmia cardiaca

V53: pulmonarydz_romano_Max flag che indica la presenza di malattie polmonare croniche

V54: coagulopathy_elixhauser_Max flag che indica la presenza di coagulo-patia

V55: compdiabetes_elixhauser_Max flag che indica la presenza di complicazioni del diabete

V56: anemia_elixhauser_Max flag che indica la presenza di deficienze legate all'anemia

V57: electrolytes_elixhauser_Max flag che indica la presenza di disturbi legati ai fluidi eelettroliti

V58: liver_elixhauser_Max flag che indica la presenza di malattia del fegato

V59: pvd_elixhauser_Max flag che indica la presenza di arteriopatia obliterante periferica

V60: psychosis_elixhauser_Max flag che indica la presenza di psicosi

V61: pulmcirc_elixhauser_Max flag che indica la presenza di disturbi della circolazionepolmonare

V62: hivaids_romano_Max flag che indica la presenza di HIV/AIDS

V63: hypertension_elixhauser_Max flag che indica la presenza di ipertensione

V64: weight peso delle comorbidità

V65: tipo di struttura del primo ricovero dell’evento (Aziende

n27

TIPO_ENTE_APP_AMM_DESC Ospedaliera, Aziende Sanitarie Locali, etc.)

V66: ENTE_APP_TERRIT_ID Codice dell'ASL della struttura del primo ricovero dell’evento,il codice indica la posizione geografica dell'ASL, in generalela provincia e l'indirizzo di riferimento

V67: ENTE_APP_TERRIT_DESC provincia dell'ASL della struttura del primo ricoverodell’evento

V68: str_lomb flag che indica se la struttura del primo ricovero dell’evento èin Lombardia

V69: primo_ev flag che indica che si tratta del primo evento del paziente; dausare insieme a “ultimo_ev” (V69)

V70: ultimo_ev vedi “primo_ev” (V70)

V71: evento_numero progressivo di evento del paziente; va da 1 al numero totale dieventi del pz con passo 1

V72: gg_inizio_out giornate da inizio evento a uscita dallo studio; i valori negativipari a -1 o -2 sono stati posti a zero, gli altri valori negativisono posti a missing

V73: gg_fine_out giornate da fine evento a uscita dallo studio; i valori negativipari a -1 o -2 sono stati posti a zero, gli altri valori negativisono posti a missing

V74: gg_inizio_inizio giornate da inizio evento precedente a inizio evento attuale; ivalori negativi pari a -1 o -2 sono stati posti a zero, gli altrivalori negativi sono posti a missing

V75: gg_fine_inizio giornate da fine evento precedente a inizio evento attuale; ivalori negativi pari a -1 o -2 sono stati posti a zero, gli altrivalori negativi sono posti a missing

V76: gg_inizio_dec questa variabile è uguale alla variabile “gg_inizio_out” (V72)se il paziente è deceduto; i valori negativi pari a -1 o -2 sonostati posti a zero, gli altri valori negativi sono posti a missing

V77: gg_fine_dec questa variabile è uguale alla variabile “gg_fine_out” (V73) seil paziente è deceduto; i valori negativi pari a -1 o -2 sono statiposti a zero, gli altri valori negativi sono posti a missing

V78: flag_pediatrico_ref flag che indica se il paziente aveva età inferiore a 2 anni alprimo evento

V79: flag_pediatrico flag che indica se il paziente aveva età inferiore a 18 anni alprimo evento

V80: gruppo1 Questa variabile contiene i dati divisi in base ai quattro gruppiidentificati nella sezione 4.2

V81: strutt_ric_id variabile da non usare

V82: strutt_ric_ricod variabile da non usare

V83: flag_shock sintomo/sindrome di shock circolatorio, una condizione dovec'è insufficienza di ossigeno per l'attività metabolica celulare

V84: flag_PTCA_Max angioplastica coronarica percutanea transluminale, una

n28

intervenzione di consente che di riprendere la circolazionesanguigna per vasi sanguini ostruiti, tipicamente dallapresenza di grasso

V85: flab_cabg intervento chirurgico con bypass aortocoronarico, unaintervenzione che consente di deviare il sangue da una aortaostruita verso il cuore

V86: flag_MDC5_Max evento classificabile come MDC5, ovvero diagnosi relativa alsistema circolatorio

Il peso delle comorbidità (variabile V64) è calcolato usando la tabella di Gagne et al.(2010)[102]:

Comorbidità Peso

Cancro metastatico 5

Scompenso cardiaco congestizio 2

Demenza 2

Insufficienza renale 2

Perdita di peso 2

Emiplegia 2

Abuso di alcool 1

Tumore 1

Aritmia cardiaca 1

Malattia polmonare cronica 1

Coagulopatia 1

Complicazioni del diabete 1

Deficienze legate all'anemia 1

Disturbi legati ai fluidi e elettroliti 1

Malattia del fegato 1

Arteriopatia obliterante periferica 1

Psicosi 1

Disturbi della circolazione polmonare 1

HIV/AIDS -1

Ipertensione -1

n29

Analisi della letteratura

4.1 Ricerca Strutturata

Per quanto riguarda il tema di questo lavoro, è stata fatta una ricerca strutturata suldatabase Scopus di articoli originali o review della letteratura:

• Usando come chiave di ricerca: “Heart failure” AND “cost” AND “data” inArticle Title, Abstract, Keywords

• pubblicati dal 2013 al 2003• risultati sono stati integrati pubblicati prima del 2003

Gli articoli sono stati scelti in base alla lettura dell'abstract e classificati in base alcontenuto in 4 gruppi:

• Articoli riguardanti l'utilizzo di risorse ospedaliere/componenti di costo dellacura dello scompenso cardiaco, importanti per capire quali sono i servizi checoncorrono in larga parte alla determinazione del costo della cura delloscompenso

• Articoli riguardanti il tema del disease management/ home care del paziente eriguardante utilizzo di risorse/costo della cura per le suddette problematiche,importante perché c'è una tendenza al trattamento preventivo dei pazienti fuoridall'ospedale per evitare riospedalizzazioni

• Articoli riguardanti metodi matematici/previsionali usati per predire ilcosto/utilizzo futuro di risorse ospedaliere

• Articoli riguardanti l'utilizzo di database amministrativi per lo studio delloscompenso cardiaco, già usati nella sezione 2.2

La maggior parte degli articoli sono strutturati come report di statistica descrittivaoppure con la struttura dello studio statistico. Uno studio statistico è solitamentestrutturato con il seguente schema:

1. Disegno dello studio: descrizione di come sono stati ottenuti i dati e dellapopolazione di pazienti

2. Metodi di analisi usati: descrizione sintetica del modello matematico usato esoftware usato per effettuare l'analisi

3. Risultati ottenuti: statistiche descrittive, analisi di trend, risultati di survery etabelle di covariate se si tratta di un'analisi di regressione

4. Discussione dei risultati

n30

Nella ricerca statistica, gli studi si dividono in 2 grandi gruppi: negli studisperimentali il ricercatore interferisce nell'oggetto dello studio, mentre negli studiosservazionali il ricercatore si limita ad osservare delle variabili, questi ultimi sonoquelli prevalenti nell' epidemiologia e nella biostatistica. Per capire i metodi dianalisi usati dagli autori è opportuno chiarire come vengono analizzati dati mancanti,censored e le tipologie di studi presenti in letteratura:

Studi sperimentali:

• Controllo per placebo: in medicina usato per valutare l'efficacia di un farmacosomministrato a un gruppo di pazienti rispetto a un gruppo a cui vienesomministrato un placebo (sostanza presunta inerte)

• Controllo per randomizzazione: pazienti assegnati casualmente a diversi gruppiper tenere il bias del ricercatore sotto controllo

• Esperimento cieco: i risultati vengono analizzati da individui che non sono aconoscenza di come è stato realizzato l'esperimento

Studi osservazionali:

• Case control[37]: molto usato in epidemiologia, i dati vengono divisi in 2gruppi diversi e esaminati retrospettivamente, per esempio: utilizzo di tabaccoin pazienti con cancro (case) e senza cancro (control)

• Cross-sectional[37]: dati ottenuti da una popolazione in cui non si considera ladimensione temporale

• Longitudinale[36]: variabili osservate ripetutamente nel lungo del tempo. Se idati sono scelti da un sottoinsieme specifico, allora si dice che è uno studio dicoorte o panello.

In questi studi possono essere usate diverse metodologie di campionamento, tra cui:

• Campionamento casuale semplice: osservazioni prelevate casualmente con lastessa probabilità

• Campionamento casuale sistematico[35]: se una popolazione è costituita da Nindividui, e il campione che si vuole prelevare è di dimensione n, allora k=n/N(eventualmente arrotondato) è l'intervallo di campionamento (skip) con cuibisogna prelevare le osservazioni dalla popolazione

n31

4.2 Componenti di costo e utilizzo di risorse ospedaliere

Prima di esaminare il costo di cura dello scompenso, cerchiamo di capire ladimensione del fenomeno scompenso cardiaco. Annualmente la American HeartAssociation(AHA) pubblica nel journal Circulation un report con il titolo “HeartDisease and Stroke Statistics”, il quale è citato con molta frequenza in letteratura (piùdi 200 citazioni5). La Association for the Prevention and Relief of Heart Disease èstata fondata a New York nel 1915, per colmare l'assenza di informazioni sul temadelle malattie cardiovascolari: il primo risultato furono una serie di studi per scoprirese i pazienti affetti da malattie cardiovascolari potevano ritornare a lavorare. Icardiologisti Lewis A. Conner, Robert H. Halsey, Paul D. White, Joseph Sailer,Robert B. Preble, Hugh D. McCulloch fondarono la American Heart Association nel1924, riconoscendo la necessità di approfondire gli studi su questo tema[114]. Trai ivari dati divulgati in questo report, ci atteniamo a quelli riguardanti lo scompensocardiaco (heart failure), il quale è una delle 5 principali cause di morte per malattiacardiovascolare negli Stati Uniti, insieme alla coronopatia, infarto miocardico,ipertensione e malattie delle arterie (I10-I15 nella codifica ICD-10).

Percentuale di morti con malattia cardiovascolare negli USA (2009) Fonte: National Heart, Lung,and Blood Institute from National Center for Health Statistics reports and data sets [104]

La mortalità è ovviamente legata all'età del paziente, in particolare è una funzioneconvessa in relazione alla variabile età:

5 Il report del 2010 ha ottenuto 2315 citazioni su Scopus

n32

Mortalità ospedaliera: maschi (blu) e femmine (rosso). Fonte: Dati del NationwideInpatient Sample6 [111]

Osserviamo come l'attributo maschio e femmina non sembra essere particolarmentecaratterizzante in relazione alla mortalità ospedaliera, tranne per le femmine con etàinferiore a 29 anni che hanno una mortalità significativamente più bassa. Però, comevedremo in seguito, i pazienti con età inferiore a 40 anni sono una piccola parte dellapopolazione scompensata americana.

Il report della AHA contiene anche dei dati riguardanti l'incidenza e la prevalenzadello scompenso cardiaco: in epidemiologia la prevalenza indica la proporzione diuna popolazione soggetta a una certa condizione, invece l'incidenza è il tasso dicrescita della condizione nella popolazione. Notiamo, dai grafici della prossimapagina, che gli uomini sono maggiormente affetti dallo scompenso, tranne perpazienti con età superiore a 80 anni, e che per i pazienti con età superiore ai 65 anni iltasso di crescita della malattia è più grande per gli uomini.

6 Il Nationwide Inpatient Sample fa parte della famiglia di database denominata Healthcare Cost and Utilization Project, sponsorizzata dalla AHRQ

n33

Prevalenza dello scompenso cardiaco in base al sesso e all'età (2007-2010). Fonte: National Centerfor Health Statistics. CVD includes International Classification of Diseases[104]

Incidenza dello scompenso cardiaco in base al sesso e all'età (1982-2003). Fonte: National Centerfor Health Statistics. CVD includes International Classification of Diseases [104]

n34

Lo scompenso cardiaco può manifestarsi attraverso due patogenesi principali: ladisfunzione sistolica, un'incapacità del cuore di contrarsi efficamente, e ladisfunzione diastolica, difficoltà dei ventricoli a riempirsi adeguatamente. Hogg et al.(2004)[108] hanno osservato, come vediamo dal prossimo grafico, che ai casi gravi dientrambi le patologie sono associati a una percentuale significativa di pazienti affettidallo scompenso cardiaco.

Prevalenza dello scompenso cardiaco nella popolazione dell'Omlsted County (Minnesota) in base a2 patologie: disfunzione sistolica e diastolica. Fonte: Olmsted County Study, Minnesota[108]

Come abbiamo accennato prima nel punto 2.3, la durata di ospedalizzazione èstoricamente associata al costo di cura dello scompenso cardiaco. Per capire comepossa variare questa variabile, facciamo un confronto tra diversi paesi:

Tempo medio di ospedalizzazione (giorni) nei diversi paesi. Fonte: Acute Study of ClinicalEffectiveness of Nesiritide in Decompensated Heart Failure [109]

n35

Osserviamo che mentre gli Stati Uniti sono tra i paesi con più basso tempo diospedalizzazione medio, l'Italia è tra i paesi con tempo di ospedalizzazione più alto. Èopportuno notare che se da un lato è opportuno ridurre il tempo di ospedalizzazione,dall'altro lato è importare assicurarsi che ciò non impatti negativamente nella qualitàdel ricovero del paziente. Oltre alla localizzazione geografica, esistono altre variabiliin base alle quali è possibile discriminare, come quelle viste nel grafico sottostante,ottenute dal Kids’ Inpatient Database (progetto HCUP):

Durata media di ospedalizzazione per bambini scompensati. Fonte: Kids’ InpatientDatabase[112]

Le 4 variabili di segmentazione con tempo medio di ospedalizzazione più alto sono:dispositivo di assistenza ventricolare e ossigenazione extracorporea a membrana (2tecnologie), setticemia e insufficienza renale acuta (2 malattie). Il NationwideInpatient Sample ci fornisce anche dati sulla durata di ospedalizzazione in base all'etàdel paziente, vediamo dal grafico B della prossima pagina che i pazienti più giovanidi 29 anni hanno una durata di ospedalizzazione più alta.

n36

Durata di ospedalizzazione: maschi (blu) e femmine (rosso). Fonte: Dati del NationwideInpatient Sample [111]

Come vediamo dal grafico sottostante, la grande parte dei costi di cura sono quelliospedalieri anche se è presente una parte non trascurabile di costi legati al homecare/disease management, sostenuti con l'obbiettivo di diminuire la ri-ospedalizzazione del paziente, e quindi auspicabilmente di contenere i costi cumulatidi cura nella vita del paziente.

Ripartizione di costi per la cura dello scompenso in USA (2008). Fonte: European Society ofCardiology [105]

n37

Concentrandoci sulla componente ospedaliera del costo per la cura dello scompenso,che sarà oggetto delle analisi del presente lavoro di tesi, e utilizzando come fonte ildatabase PREMIER vediamo la distribuzione (asimmetrica) del costo dei ricoveri, ealcune delle variabili che hanno un impatto significativo su di esso:

Distribuzione di costo totale (A) per la unità di “acute care” (2004-2005) con codifica ICD-9 per scompenso e età superiore a 18 anni. I dati sono estratti dal database americano PREMIER. Fonte: PREMIER's Perspective Comparative Database [107]

n38

Distribuzione di costo (popolazione precedente) in base al reparto. VasoA=medicinali vasoattivi.Fonte: PREMIER's Perspective Comparative Database[107]

Per “Room & Board” si intendono quei costi necessari per l'utilizzo di una camera daletto e il cibo servito al paziente, questi costi sono più contenuti per quei pazientisoggetti a medicinali vasoattivi, che dall'altro lato sono responsabili per costi più altinel' unità di cura acuta. Per i pazienti ospedalizzati con più di 5 giorni abbiamo unconsumo più alto nello stesso reparto.I costi di cura dello scompenso cardiaco continuano ad aumentare negli ultimi anni(uno dei motivi di interesse su questo tema a livello americano e internazionale)anche se tale costo in proporzione al costo totale per la cura di malattiecardiovascolari sembra rimanere stabile:

Costi per lo scompenso cardiaco negli USA (2003-2010). La linea piena indica il valore in miliardidi dollari, mentre la linea tratteggiata indicata il valore sulla percentuale dei costi totali per lemalattie cardiovascolari. Fonte: European Society of Cardiology [105]

n39

L'AHA proietta per il futuro un aumento dei costi di cura per le malattiecardiovascolari, tra cui è presente anche lo scompenso cardiaco:

Costi proiettati per malattie cardiache per il 2015-2030 dalla AHA. Fonte: Unpublished datatabulated by the American Heart Association [104]

Consideriamo che il PIL americano del 2013 è compreso tra i 16.000 e i 17.000miliardi di dollari, in particolare la stimativa della banca mondiale è 16.800 miliardidi dollari[115], quindi per effetto di comparazione abbiamo che i costi totali proiettatiper il 2015 sono pari al 2.13% del PIL del 2013, e i costi dello scompenso il 9% deicosti totali per la cura delle malattie cardiovascolari.Sotto è riportato la ripartizione dei costi di cura in base all'etiologia per un'ospedale-università della Repubblica Ceca:

n40

Ripartizione dei costi in base alla etiologia per la University Hospital Brno. Fonte: AHEAD Coreregistry[110]

La legenda è la seguente:

• new-onset:scompenso cardiaco di nuova insorgenza• ADCHF: acute decompensation of chronic heart failure• ACS: acute coronary syndrome• CIHD: chronic ischemic heart disease• valvular dysfunction: valvulopatia• arrhythmia: artmia• HT crisis: hypertensive crisis• other: altri• ICU-Intensive Cardiology Unit: Unità di cura cardiologica intensiva• SCU- Standard Cardiology Unit: Unità di cura cardiologica standard• CAG- Coronary AngioGraphy: angiografia coronarica• PCI- Percutaneous Coronary Intervention: intervento coronarico percutaneo• antiarrhytmic interventions: interventi anti-aritmia

È da notare che la fetta di costo maggiore è sempre quella dell'angiografia o degliinterventi anti-aritmia, eccezione fatta per la valvulopatia dove prevalgono i costidell'unità di cura cardiologica intensiva.

n41

4.3 Dati mancanti e censored

Prima di passare alla descrizione dei modelli matematici è opportuno tenere presenteuna problematica comune negli studi epidemiologici, ovvero la presenza di datimancanti e/o censored. Little e Rubin (2002) [44] offrono una guida ai possibilimeccanismi che possono generare dati mancanti, sia Y=(yij) un generico dataset,M=(Mij) una matrice binaria che indica la presenza di osservazioni mancanti e f(M|Y;φ) una densità di probabilità condizionata con uno(o più) parametro(i) incognito(i) φ:

• Missing completely at random (MCAR): , ovvero sela struttura M dei valori mancanti non dipende da Y

• Missing at random (MAR): , ovvero se lastruttura M dei valori mancanti dipende solo dai valori osservati Yobs

• Not missing at random (NMAR): , ovvero se lastruttura M dei valori mancanti dipende anche dai valori mancanti del dataset Y

I metodi proposti per queste problematiche sono classificabili nelle seguenti categorienon mutuamente esclusive:

• Procedimenti basati su unità statistiche completamente osservate: scartare leosservazioni mancanti e analizzare solo i dati completi.

• Procedimenti di pesatura: i parametri o i momenti incogniti vengono stimaticon dei pesi che specificano la probabilità di una unità statistica essere o nonessere inclusa nel campione.

• Procedimenti di imputazione: i valori mancanti vengono sostituiti e ilcampione viene analizzato con metodologie standard. L' imputazione puòessere: “hot deck”, quando i valori mancanti vengono sostituiti usandoosservazioni del campione, “mean imputation”, quando i valori mancantivengono sostituiti usando le medie di sottoinsiemi del campione, “regressionimputation”, quando i valori mancanti (di un responso) vengono stimati conuna funzione di regressione.

• Procedimenti basati sulla specificazione di un modello: tipicamente si specificauna funzione di densità o di verosimiglianza per i dati, e poi si stimano iparametri incogniti con un metodo appropriato.

Il problema del censoring è un caso particolare di dati NMAR, che presenta lacaratteristica peculiare che i dati vengono tagliati dal campione in funzione di unavariabile temporale, quindi limitandoci allo scopo di questo lavoro, la mancanza difollow-up del paziente fino al momento della morte rappresenta un problema di rightcensoring. In un lavoro di review su questo tema, Wijeysundera et al.(2012)[38]individuano 4 proprietà importanti dei dati di costo clinici:

n42

• Grande parte della popolazione generale è sana e ha costo pari a zero• Distribuzione asimmetrica con pochi individui ad alto costo nella coda destra• Violazione della proprietà di omoschedasticità• Dati incompleti

Le funzioni di costo incidenti per un singolo paziente sono convesse (grafico A nellaprossima pagina), quindi se il follow-up del paziente non comprende la morte, i costirischiano di venire sottostimati anche di parecchio visto che una parte rilevante dicosti è quella della fase terminale (grafico B). Se il periodo di osservazione termina inmodo casuale allora abbiamo a che fare con un problema di random right censoring.

Fonte: Wijeysundera et al. (2012)[113]

Gli autori rassegnano un insieme di tecniche sviluppate per il problema del rightcensoring che verranno descritte in seguito.Sia S(t)=Pr(tempo di vita del paziente<t) la funzione di sopravvivenza e0≤t1'≤t2'≤...≤tfinale' la sequenza di osservazioni dei tempi di vita, allora nel metodo diKaplan-Meier[39] lo stimatore è:

Dove tr'≤t, si chiama funzione dell' hazard, che èl'espressione matematica del tasso di mortalità. Se Sc(t)=Pr(tempo che il pazienterimane uncensored>t) è la probabilità di essere uncensored dopo il tempo t eanalogamente:

n43

Un requisito fondamentale di questo metodo è un censoring indipendente[38], ovvero:

• il “tempo di censoring” deve essere stocasticamente indipendente dal “tempo fino alla morte”.

• per quanto riguarda i costi cumulati, i “costi fino al momento del censoring” difficilmente sono stocasticamente indipendenti dai “costi fino alla morte”.

Per ogni paziente i la funzione di costo cumulata è Citotal=Ai(ti), esiste una famiglia di

metodi che ripesano le osservazioni in modo tale che ogni caso completo rappresenti anche i casi censored. Lo stimatore di Lin et al. (1997)[40] si basa sulla suddivisione dell' intervallo di osservazione in K intervalli ugualmente grandi, se solo i costi C=Ci

total sono disponibili allora:

Dove T è il tempo latente di sopravvivenza e Ti è il tempo latente di sopravvivenza per il paziente i.Bang e Tsiatis (2000)[41] descrivono uno stimatore a pesatura con probabilità inversa, in inglese inverse probability weighting (IPW), 1/Sc(Ti) che ripesa ogni osservazione uncensored in modo tale da ottenere:

Δi è una variabile binaria che assume valore 1 se il paziente è uncensored, altrimenti 0se è censored e N è il numero di pazienti presenti nello studio. Un problema dello stimatore IPW è l'inefficienza statistica, Raikou e McGuire (2004)[42] hanno effettuato un'analisi di simulazione nella quale hanno determinato che nella presenza di un fattore di censoring superiore al 50% lo stimatore IPW diventa instabile.Una variante “partioned” di questo ultimo stimatore consiste nel dividere l'intervallo di osservazione in K partizioni e Mi

j indica la porzione di costo per ogni partizione j, allora:

n44

Quest'ultimo è il metodo consigliato da Wijeysundera et al. (2012) per la suaefficienza. Un altro metodo è quello del phase based costing che consiste nei seguentipassi:

• Definire le fasi della malattia. Per esempio, analogamente a quello che si fa nella letteratura di ingegneria dell'affidabilità[43], si può usare una curva da vasca da bagno per descrivere l'andamento dei costi del paziente: una fase iniziale dove i costi decrescono, una fase intermedia con costi pressoché stabili, e una fase terminale dove i costi crescono rapidamente fino al momentodella morte.

• Siccome la curva da vasca da bagno è convessa per costruzione si possono usare i punti di flesso della cumulata C=Ci

total per determinare la durata delle diverse fasi

• Allocare i costi dei pazienti ad ogni fase• Determinare il costo medio per ogni fase• Usando le medie dei dati di costo e la funzione di sopravvivenza S(t) costruire

la curva cumulata dei costi

4.4 Modelli matematici previsionali

Un modello matematico è una astrazione selettiva di un sistema reale[3], un buonmodello è quindi quello che riesce a rappresentare le principali caratteristiche di unsistema reale senza l'uso di una eccessiva complicazione e senza l'uso di ipotesi nonrealistiche. Essendo la matematica una disciplina deduttiva tutte le proprietà delmodello sono una conseguenza diretta delle ipotesi di modellazione. Alcuni modellisemplicemente indicano un trend futuro, per esempio di costo, senza specificare se èpossibile interferire sul fenomeno sotto studio, altri prevedono la possibilità che cisiano delle variabili decisionali su cui agire (leve). I modelli statistici contengonosempre una componente aleatoria che a volte è quella predominante.La seguente tabella è un riassunto schematico della review di letteratura sui modellimatematici:

n45

Nota: i modelli indicati con il colore verde sono tutti classificabili come modelli di regressione.

Goehler et al. 2011[116], nella review che hanno condotto su questo tema neidatabase MEDLINE/PreMEDLINE/EMBASE hanno identificato anche comemodelli:

• Modelli Markoviani• Sistemi di equazioni differenziali• Simulazione a eventi discreti

Noi ci concentriamo sugli articoli individuati nella tabella precedente, i qualiapplicano modelli matematici alla problematica dell'utilizzo (presente e futuro) dirisorse ospedaliere e/o costi clinici legati allo scompenso cardiaco.

4.4.1 Estrapolazione

Un metodo abbastanza semplice per ottenere una previsione dell'andamento dei costiclinici è quello di identificare un trend (tendenza) nei dati storici ed estrapolarequesto trend nel futuro. Fang et al. (2008) [9] prendono in analisi un campione casuale sistematico di ospedalinegli USA (dati NHDS) e calcolano il tasso di ospedalizzazione dal 1979 al 2004come rapporto tra numero stimato di ospedalizzazioni e popolazione civile

n46

americana, e poi congetturano qualitativamente che con l'invecchiamento dellapopolazione e interventi terapeutici il trend di aumento delle ospedalizzazioni dovràcontinuare nel futuro. Stewart et al. (2002) [10] usano 2 campioni di studi cross-sectional nella GranBretagna affetti da disfunzione sistolica sintomatica del ventricolo sinistro (1995) perstimare la popolazione di scompensati cardiaci nel UK. Successivamente stimano ilcosto per l'anno 2000 sulla base di estrapolazione di trend di costo (dal 1990 al 1995)di diversa natura: demografici, consultazioni di medicina generale, ospedalizzazioni epazienti discharged. Garrett et al. (2007) [67] esaminano dati provenienti dal piano assicurativoHealthPartners tra il 2002 e 2003, dati i quali vengono aggregati unità cliniche cherappresentano interi episodi di ricovero. Le malattie dei pazienti includono, oltre alloscompenso cardiaco, coronaropatia, diabete, condizioni psichiatriche, asma,dipendenza chimica e infertilità. I costi comprendono pagamenti dei pianiassicurativi, pagamenti di benefici e debiti dei membri, corretti con il tasso diinflazione Consumer Price Index (CPI).Gli autori applicano le proporzioni di costo per età, sesso e malattie alle proiezioni dicrescita della popolazione americana tra il 2000 e 2050 del US Census Bureau. [68]Quindi anche se i suddetti paper adottano l'idea del' estrapolazione di trend, notiamoalcune differenze significative:

• Fang et. al estrapolano il tasso di ospedalizzazione, Stewart et al. estrapolanodati di costo e Garrett et al. estrapolano il tasso di crescita della popolazioneamericana

• Fang et. al e Garrett et al. fanno una previsione su eventi che avvengono dopola pubblicazione del paper, mentre Stewart et al. fanno una previsione sulpassato, proiettando dati del 1995 per l'anno 2000

• Fang et. al fanno una previsione di natura puramente qualitativa (il trend diaumento potrà continuare nel futuro), mentre Stewart et al. e Garrett et al.fanno una stima quantitativa

• Tutti gli approcci ricadono comunque nella casistica dei metodi diestrapolazione di serie storiche [3]

4.4.2 Modelli di regressione

I modelli di regressione sono ampiamente usati nell'analisi di costi clinici e inepidemiologia, quindi prima di procedere è opportuno formalizzare la varietà diquesta classe di modelli, visto che questi modelli sono anche abbastanza variegati perquanto riguarda le ipotesi di modellazione. Segue per primo la trattazione classica che prende il nome di regressione linearemulti-variata. [2][3]

n47

Sia un dataset composto da m osservazioni e n+1 variabili, di cui una è il responso(variabile target o anche risposta) e le altre n variabili sono le covariate (predittori). Ilresponso verrà indicato con Y e le covariate con X1, X2, …, Xn.. Le covariate formanoquindi un vettore indicato con . Come si è soliti indicare in statistica, la lettera maiuscola è la variabile aleatoria e lalettera minuscola è il valore assunto da questa variabile sul campione: quindi ilresponso Y assume sul campione i valori yi , mentre la matrice delle osservazionidelle covariate X1, X2, …, Xn. viene indicata con X=[xij].Si ipotizza l'esistenza di una funzione che lega il responso con lecovariate: + ε , dove ε viene chiamato rumore (errore):

• Se , allora si

dice che il modello è una regressione lineare• In un campione il modello assumerà dei valori puntuali yi = ∑bi xij+ εj e sarà

indicato con • Nel linguaggio della probabilità si dice che il modello che stiamo studiando è il

valore atteso condizionato di Y dato un certo numero di osservazioni X1=x1,X2=x2,, …, Xn=xn , cioè E[Y| X1=x1, X2=x2,, …, Xn=xn]. Quindi formalmente:Y=E[Y| X1=x1, X2=x2,, …, Xn=xn]+ ε

• L'approccio classico è quello di assumere che il rumore sia normalmentedistribuito ε ~N(0,σ2)

Questo modello di base, che è la sintesi dei risultati classici di Gauss, Legendre,Galton e Pearson, è stato soggetto ad alcune generalizzazioni. La distribuzionenormale fa parte di una famiglia più ampia di distribuzioni chiamata famigliaesponenziale (classe esponenziale) [4]. Sia Y una variabile aleatoria con fdp(funzione di densità di probabilità) descritta da un solo parametro θ con forma:

Dove s, t, a, b sono funzioni note. É comune anche scrivere la fdp sostituendos(y)=exp(d(y)) e t(θ)=exp(c( θ)):

Possiamo anche considerare una fdp della famiglia esponenziale formata da piùparametri, ovvero da un vettore di parametri:

n48

Quest'ultima viene più comunemente presentata nella forma:

È da notare che per la definizione di fdp (∫f(x;θ)dx=1), che t(θ) viene logicamentedeterminato una volta scelti s(y), a(y) e b(θ) (discorso analogo per c(θ) e ).La famiglia esponenziale contiene un grande numero di distribuzioni: normale,esponenziale, gamma, chi-quadro, Poisson solo per citarne alcune.

Consideriamo, per semplicità, distribuzioni della famiglia esponenziale con un soloparametro. Siano Y1, Y2, … ,Yq q responsi appartenenti alla famiglia esponenziale :

• ogni Yi ha fdp • le fdp di ogni Yi hanno tutte la stessa forma, e con l'ipotesi di indipendenza

stocastica delle Yi la fdp congiunta è data:

• Indichiamo, come è comune, E(Yi)=μi, μi sarà una funzione di θ

Un modello lineare generalizzato è quindi dato da:

E assume valori:

g si chiama la funzione di link. L'estensione alla famiglia esponenziale con piùparametri è analoga.

Se le fdp delle Y1, Y2, … ,Yq sono delle e la funzione di link è g(x)=x allorail modello diventa:

n49

Da notare che è un modello lineare generalizzato con le Yi prese da una fdp dellafamiglia esponenziale con 2 parametri. Oppure riscritto:

Questo modello viene chiamato modello generale lineare, perché generalizza laregressione lineare multipla prima esposta al caso di multipli responsi. Si noti che ininglese il modello generale lineare si dice General Linear Model7 (GLM) e il modellolineare generalizzato si chiama Generalized Linear Model (GLM), quindi solo usarel'acronimo può essere fuorviante senza far vedere anche la forma funzionale delmodello. Il modello generale lineare in R è dato dalla funzione lm(), mentre ilmodello lineare generalizzato è dato dalla funzione glm().

Esistono diverse tecniche per ottenere il fit di un modello di regressione. Se m è ladimensione del campione:

• Se si minimizza , somma dei quadrati degli errori (Sum of

Squared Errors) , allora si dice che è una regressione ottenuta con il metodo deiminimi quadrati (least squares)

• Se si minimizza , allora la regressione ottenuta con in metodo

delle deviazioni minime (least absolute deviation)

• Se si massimizza la funzione , allora la

regressione è ottenuta con il metodo della massima verosimiglianza (maximum

likelihood). É comune cercare il massimo della , la

log-likelihood. Entrambe hanno il massimo nello stesso valore di θ.• Tutte le metodologie precedenti, quando applicate generano sempre delle

equazioni di stima (dei parametri incogniti). Esistono molte equazioni di stimaproposte in letteratura, però ci soffermiamo brevemente sulle generalestimating equations (GEE). Se N sono i responsi stocasticamente

indipendenti, le GEE vengono usate per stimare i

parametri dei modelli lineari generalizzati: Di=∂μi/∂βj , e

Ai=diag(var(yit)) sono matrici. R è la matrice di correlazione per y e φ è unacostante.

7 La Dobson nel suo Introduction to Generalized Linear Models lo chiama anche Normal Linear Model

n50

• Se i responsi e allora abbiamo laregressione di Poisson come caso particolare del modello lineare generalizzato.Passando ai logaritmi abbiamo la forma classica della regressione di Poisson

, dove il termine lnni viene chiamato termine di offset[4].

I modelli di regressione, come prima descritto, hanno una natura deduttiva: con ilverificarsi dei valori delle covariate, si deduce il valore del responso. Esempi dicovariate trovate in questa letteratura sono: età, durata di ospedalizzazione, numero diinvertenti medici sul paziente, intensità della malattia, presenza di comorbidità. Cosìcome l'estrapolazione di serie storica, i modelli di regressione consentono anche dieffettuare previsioni sul responso, purché si riescano a prevedere anche i valori dellecovariate. Da questo punto di vista si parla di interpolazione quando si prediconovalori che stanno dentro l'insieme di osservazione del responso e di estrapolazionequando si predicono valori che stanno fuori dell'insieme di osservazione del responso.La grande differenza concettuale tra i modelli di regressione e l'estrapolazione diserie storica è che i modelli di regressione suggeriscono anche su quali variabilibisogna agire per cercare di controllare il valore del responso. Quindi, per esempio,se si ritiene che la durata dell'ospedalizzazione abbia un impatto rilevante sul costo,questo suggerisce agli addetti alla gestione del sistema sanitario che per ridurre ilcosto è opportuno ridurre la durata di ospedalizzazione, con però alcuni caveat:

• Tenere conto della possibilità che riducendo la durata di ospedalizzazione delpaziente possa ridursi anche la qualità del ricovero (quindi impatto su unresponso non modelizzato)

• Tenere conto che tanto più lontana è l'estrapolazione dei dati del dataset, tantopiù rischio c'è di commettere errore: in questo caso per durate ospedalieremolto lunghe la stima di costo può risultare fuorviante

• Tenere conto che negli studi di regressione esiste il rischio di sovrastimare ilcontributo delle covariate per spiegare le variazioni del responso rispetto allavariabilità naturale del responso, fenomeno indicato in letteratura comeregressione verso la media[66]

• Tenere conto che in un modello con molte variabili, aumenta la probabilità direlazioni spurie tra le variabili, ovvero relazioni funzionali tra le variabili chenon sono replicabili nel tempo ma dovute alla pura casualità, problema ormaiinevitabile alla luce del fenomeno dei Big Data.

• Idealmente in un modello di regressione le covariate sono scorrelate tra di loro,anche perché se una covariata Xi dipende funzionalmente da una covariata Xj,ha poco senso tenere dentro il modello entrambe le variabili. Oltre a questodiscorso concettuale c'è anche il problema del calcolo della matrice inversa(X'X)-1 che diventa tanto più grave quanto quest'ultima è prossima a unamatrice singolare. Tuttavia si può cercare di rimediare questo problema

n51

numerico usando metodi di regressione ridge.

Come prima annunciato, descriveremo i modelli di regressioni proposti, i quali sonoabbastanza variegati e i contesti di utilizzo sono i più diversi.Dunlay et al. (2010) [11] realizzano uno studio di tipo coorte in Olmsted County(USA) con pazienti dal 1987 al 2006: il modello di regressione è quello di Tian eHuang (2007).[12] Il modello a due parti di Tian e Huang è stato sviluppatooriginariamente per l'analisi di dati di costo Medicare in studi di coorte. Il modello sipropone di catturare 3 elementi importanti:

• Presenza di asimmetria dei dati di costo con coda destra (right-skew)• Presenza di molti zeri per i pazienti con casi poco gravi• Presenza di follow-up incompleto (censoring)

Le variabili del modello sono:

• T indica l'istante della morte• C indica l'instante del censoring• Z=min( T, C )• M(u) indica il costo cumulato all'instante u, 0 ≤ u ≤ Z, M(u) indipendentemente

e identicamente distribuito• è il vettore di covariate• i è l'indice di paziente• Yi è il costo cumulato per il paziente i fino al tempo Ti=Ti

L =min( Ti , L ), doveL è una costante di upper bound per il periodo di analisi

Il modello a 2 parti è dato da:

g() e m() sono funzioni di link (inverse secondo il formalismo prima usato), laseconda parte è un modello lineare generalizzato. I parametri β'0 e γ'0 sono calibraticon il metodo delle equazioni di stima, che Tian e Huang costruiscono basandosi sulfilone di letteratura sulla pesatura con probabilità inversa prima descritto.Dunlay et al. hanno ottenuto il fit del modello separatamente per ricoveri ordinari ericoveri diurni:

• Per la prima parte del modello (probabilità di avere costi positivi) è stata usata

n52

una regressione logistica• Per la seconda parte (GLM) è stata usata come funzione di link un logaritmo• Per il responso si è scelto di usare la distribuzione gamma, motivata dal test di

Park modificato[14]• Solo costi diretti

I dati mancanti in questo studio erano pochi (<2%) tranne per i dati sulla frazione dielezione (ejection fraction) (26% dati mancanti), questi ultimi sono stati corretti con ilmetodo dell'imputazione multipla [45]. In particolare è stato creato un vettore di 5elementi per ogni osservazione mancante (quindi per un totale di 5 dataset) e poianalizzati con il metodo suggerito da Donald Rubin, ovvero: applicare metodistandard separatamente ad ogni dataset per stimare i parametri incogniti, e poiottenere la stima per il parametro incognito del campione originale come media delle5 stime fatte per ogni dataset. Per tutti gli articoli di questa sezione verrà riportat unatabella come la seguente, dove sono evidenziate le covariate testate dagli autori e, sedisponibile, i p-value calcolati:

Autori Dunlay et al.Data pubblicazione 2010Età 0.0001Sesso 0.375Massa corporea 0.642Creatina pre operatoria 0.314Ipertensione 0.862Malattia cerebrovascolare 0.003Diabete mellito 0.003Malata polmonare cronica 0.576Infarto miocardico 0.411Frazione di eiezione 0.041

Esposito et al. (2009)[15] prendono in esame dati Medicare e Medicaid di 4 statiamericani: Arksansas, California, Indiana e New Jersey. La regressione è effettuatacon un modello a 2 parti[16]:

Come al solito, la seconda parte del modello è un modello lineare generalizzato. Ilresponso è il costo medico.

n53

Autori Esposito et al.Data pubblicazione 2009aderenza ai medicinali 0.01

Titler et al. (2008)[17] studiano dati di un centro del Midwest americano (1998-2002)usando i metodi di stima di Liang e Zegler[18]:

• Yi indica il costo per il paziente i, che appartiene alla classe esponenziale e che assume valori yit

• indica il vettore di covariate per il paziente i e che assume valori •

Se le Yi presentano indipendenza stocastica per ogni paziente, allora si può usare ilmetodo di stima delle GEE. Gli autori usano anche un test di correlazione ρXY diordine zero8 tra il costo ospedaliero e ognuna delle covariate: se il p-value per lacandidata a covariata era maggiore di 0.15, questa è stata scartata.

Autori Titler et al.Data pubblicazione 2008Anemia 0.483Numero di ospedalizzazioni 0.0001Gravità malattia 0.1777Rapporto infermiera professionale per paziente 0.0001numero di procedure 0.0001diagnosi cardiaca non invasiva 0.002diagnosi cardiaca invasiva 0.067terapia cardiaca invasiva 0.001numero di medicine 0.001Benzodiazepine 0.0032medicinali per occhi, orecchie, naso e gola 0.0001antiobici per occhi, orecchie, naso e gola 0.0001Antiemetico 0.382medicinali gastrointestinali 0.0018Cefalosporine 0.0004Diuretici risparmiatori di potassio 0.0022antibiotici cutanei 0.0001antifungino 0.0012irrigazione 0.0002vitamina k 0.148agenti del sistema nervoso centrale 0.0036

Jayadevappa et al. (2006)[21] analizzano retrospettivamente uno studio case-controldi una unità di cura acuta per anziani (età≥65) (in inglese ACE) dal 1999 al 2000, le

8 La correlazione (parziale) di ordine uno ρXY|Z è la correlazione dei residui di X e Y su Z, quindi

una correlazione che tiene conto di una variabile di controllo, su R si può usare la funzione cor.test()per la correlazione di ordine zero e pcor.test() per la correlazione di ordine uno.

n54

codifiche primarie sono di scompenso cardiaco, pneumonia e infezione delle vieurinarie. I dati sono stati ottenuti dal Pennsylvania Integrated Clinical andAdministrative Research Database (PICARD), i costi sono stati stimati con unrapporto cost-to-charge pari a 0.8 e alcuni test sono stati applicati sui dati di costo:

• Test t di Student[2][22] sul logaritmo dei dati di costo• Test bootstrap per il costo medio[23] (ACE vs non ACE)• Test dei ranghi con segno di Wilcoxon per il costo mediano[24] (ACE vs non

ACE)

Una regressione loglineare è stata usata per analizzare il costo aggiuntivodell'unità ACE. Una tecnica di smearing[25] è stata usata per correggere il biaspresente nella ritrasformazione, cioè per ritornare dalla scala logaritmica lnY allascala originale Y: sia g la trasformazione, h=g-1 la ritrasformazione, allora si puòstimare usando:

Gli autori hanno anche usato una regressione di Poisson per analizzare i numeri diriospedalizzazioni annue.

Autori Jayadevappa et al.Data pubblicazione 2006Dialisi post operatoria 0.95Ecocardiografia 0.43medicinali gastrointestinali 0.53laboratorio 0.4terapia respiratoria 0.23farmaci 0.62radiologia 0.94costo della camera 0.003rifornimento medico 0.04Scintigrafia polmonare 0.7Laboratorio prove di funzionalità respiratoria 0.47terapia fisica 0.005telemetria 0.12nutrizione 0.86terapia vocale 0.36

Dall et al. (2010)[65] esaminano un campione di pazienti provenienti dal programmaTRICARE, un programma di assistenza medica per i militari americani. Questocampione esclude i pazienti che hanno abbandonato il programma TRICARE, i

n55

pazienti con meno di 6 mesi di assistenza, i pazienti che sono passati al programmaMedicare, infetti da HIV/AIDS, oppure morti.Lo studio è di tipo case-control, che presenta un problema non banale nella scelta delgruppo control, visto che usare non-partecipanti al programma TRICARE comecontrol comporta il rischio di efetto di selezione (“selection bias”), ovvero creare unadistorsione nella stima degli indicatori, mentre una comparazione tra un prima e undopo (pre-post evaluation) comporta il rischio di sovrastimare variazioni strutturaliper l'effetto di regressione verso la media.[66] Gli autori hanno optato per creare ungruppo denominato come Historical Control Group (HCG) fatto di pazienti affetti daasma, scompenso cardiaco congestizia o diabete.L'analisi di regressione è stata condotta con dati tra l'ottobre 2004 e settembre 2005:una regressione di Poisson per ricoveri ambulatoriali e non-ambulatoriali, un modelolineare generalizzato con distribuzione gamma per le spese mediche e una regressionelogistica per le ricette e i farmaci.

Autori Dall et al.Data pubblicazione 2010Scompenso cardiaco congestizio *Durata ospedalizzazione *visite di emergenza *ambulatorio *prescrizioni 0.05

Kaul et al. (2011)[26] prendono in esame una coorte di pazienti scompensati residentinell'Alberta (Canada) con almeno una ospedalizzazione oppure almeno 3 diagnosi discompenso ambulatoriali in un periodo consecutivo di 3 mesi, i costi usati nellaregressione sono: costi del personale medico, costi dei farmaci, costi diospedalizzazione, costi di diagnosi, costi di procedimenti chirurgici, costi diinterventi, costi di protesi. Gli autori hanno usato i test di Cochran-Mantel-Haenzel[27][28][29] per l'analisi di trend temporali e un modello linearegeneralizzato con funzione di link logaritmica e responso con distribuzione diPoisson. Gli autori hanno effettuato la regressione solo sugli ultimi 6 mesi di vita deipazienti, in ciò che loro chiamano costi “End-of-Life” (EOL).

n56

Autori Kaul et al.Data pubblicazione 2011Età *Sesso *Ipertensione *Arteriopatia obliterante periferica *Malattia cerebrovascolare *Diabete mellito *Durata ospedalizzazione *Insufficienza renale *cancro *Broncopneumopatia cronica ostruttiva *Demenza *anno morte *

Setoguchi et al. (2010)[19] simultaneamente analizzano dati di scompenso cardiaco ecancro, usando come fonte di dati il programma Medicare, registro statale di cancro edati farmaceutici di erogazione del programma di assistenza della Pennsylvania (PA).La coorte per lo scompenso cardiaco è stata definita usando solo pazienti conospedalizzazioni≥2 nel PA (1997-2004), solo pazienti la cui è morte è stata permalattie cardiache e solo pazienti che hanno assunto diuretici dell'ansa 365 giorniprima della morte. L'analisi dei dati è stata effettuata con una regressione di Poissonmodificata (da Zou 2004[20]): se xi è una esposizione a un presunto fattore di rischio(covariata), yi una variabile binaria che rappresenta l'accadimento dell'eventorischioso e π(xi) è il rischio per il paziente i, possiamo usare una funzione di linklogaritmica nella regressione:

Questo tipo di problematica può essere riassunto in una tabella di contingenza:yi=1 yi=0 Totale

xi=1 a b n1=a+b

xi=0 c d n0=c+d

n=n1+n0

E di conseguenza gli autori stimano il rischio relativo:

Whitfield et al. (2006)[32] usano dati di 5 primary care trust (PCT) del Regno Unitocome input dell'equazione di Framingham[33]. Siano x1,...,xk k fattori di rischio perun individuo, se T è il tempo necessario per il verificarsi dell'evento rischioso

n57

(infarto, scompenso cardiaco, ischemia,...) che segue una distribuzione di Weibull e laregressione :

Quest'ultima equazione stima il rischio dell'evento cardiaco usando la distribuzione diGumbel. Whitfield et al. affermano inoltre che siccome i pazienti con il diabete hannoun rischio cardiaco più elevato è più opportuno usare il modello UKPDS[34] (cheperò non riguarda lo scompenso cardiaco, ma altre malattie cardiache). Gli eventivengono valorizzati con il sistema inglese Health Resource Grouping (HRG). Lecovariate usate non sono disponibili nel paper.

Howel et al. (2009)[30] studiano dati del Queensland Hospital Admitted PatientsData Collection (QHAPDC) con pazienti ospedalizzati in emergenza, aventi comegruppo 28 malattie diverse (tra cui scompenso cardiaco) dal 2005 al 2006. I pazientisono stati identificati con il sistema AR-DRG e gli autori usano una regressionelogistica[3] per stimare la probabilità di riammissione del paziente:

Gli autori usano 75% del dataset come training per il modello e l'altro 25% per lavalidazione del modello. Gli autori affermano che i pazienti a maggior rischio diriammissione sono anche quelli più costosi. Le covariate sono state selezionate con un metodo di selezione intenzionale(purposeful selection) per ridurre la colinearità[31]:

1. Analisi univariata di ogni variabile candidata al modello: per le variabilicategoriche si può usare un test chi-quadro per il calcolo del rapporto diverosimiglianza, per le variabili continue si può usare il test di Wald. Se il p-value del test è <0.25, allora la variabile viene usata nel prossimo passo.

2. Usando le variabili selezionate prima si ottiene un fit con un modello logisticomultivariato, usando ancora una volta il test di Wald si calcola il p-value perogni covariata. Le variabili che non contribuiscono, con livelli tradizionali disignificativa statistica, vengono eliminate generando un nuovo modello piùpiccolo. Si esegue un test sul rapporto di verosimiglianza tra modello nuovo emodello vecchio.

n58

3. Se Δβ>20% la covariata viene re-inclusa nel modello e si ritorna al passo 2. SeΔβ <20% si va al passo 4.

4. Se è rimasta qualche covariata con Δβ>20%, si ritorna al passo 3. Altrimenti siottiene un modello preliminare per gli effetti.

5. Per ogni variabile continua verificare che il modello cresce linearmente con lacovariata, se questo procede, allora abbiamo il modello principale per glieffetti.

6. Si esegue un check sulla possibilità che ci siano variabili che interagiscono tradi loro, se si aggiunge un termine di interazione si ritorna al passo 2. Altrimentiabbiamo il modello finale.

7. Verificare la goodness-of-fit del modello finale.

Per lo scompenso cardiaco è stato calcolato un p-value <0.0001.

Zhongmin et al. (2012)[70] prendono in esame dati dal California CABG9 OutcomesReporting Program (CCORP) collegati (linkage) ai dati discharge degli ospedali dellaCalifornia dal 2009 al 2010. Il linkage è di natura probabilistica[77] con un tasso disuccesso del 99.5%. Gli autori studiano le riospedalizzazioni dopo CABG solo se ipazienti sono stati riammessi a una unità di cura acuta con un codice principale didiagnosi ICD9-CM per malattie del cuore. Per pazienti riospedalizzati multiple voltenell'arco di 30 giorni è stato considerato solo il primo evento. Siano i responsi binariYij=pij+eij che indicano se il paziente i è stato riammesso all'ospedale j, con i l'indice dipaziente e j l'indice di ospedale e i X1,...,Xk k fattori di rischio, la regressione logisticaè data da:

Si ha che Y~Bin(1, πijX), eij indipendenti, E(eij)=0, Var(eij)=σ2e=pij(1-pij) .

Gli autori hanno ottenuto un fit anche con una versione a 2 livelli del modello logitprima esposto (regressione logistica gerarchica):

Con uj~N(0,σ2u) iid che sono variabili caratteristiche dell'ospedale. Gli autori

calcolano per questo modello un coefficiente di correlazione intraclasse (ICC), cioè

9 CABG=Coronary Artery Bypass Grafting, un procedimento chirurgico usato per superare condottivascolari ostruiti

n59

un coefficiente che stima quanto gli elementi di uno stesso insieme si assomiglianotra di loro: in questo caso la percentuale di varianza residua dovuta alle differenze tragli ospedali:

Autori Zhongmin et al.Data pubblicazione 2012Scompenso cardiaco congestizio 0.0001Età 0.0001Sesso 0.0001Etnia 0.0001Massa corporea 0.002Status Procedura (Elettiva, Urgente o Emergente) 0.0001Creatina pre operatoria 0.003Ipertensione 0.0001Arteriopatia obliterante periferica 0.0001Malattia cerebrovascolare 0.0001Ictus 0.0001Diabete mellito 0.0001Malata polmonare cronica 0.0001Trattamento immunosoppressivo 0.023Dialisi 0.0001Flutter atriale 0.0001Blocco cardiaco 0.235Tachicardia 0.096Infarto miocardico 0.0002Shock cardiogeno 0.01Classe NYHA 0.0001Chirurgia cardiaca 0.137Angioplastica 0.0572Frazione di eiezione 0.0001Stenosi sinistra 0.165Numero arterie malate 0.551Insufficienza mitrale 0.0004Endocardite infettiva 0.439Rianimazione 0.439Durata ospedalizzazione 0.0001Assicurazione 0.0001Rendita mediana 0.0001Sternale profondo 0.0001Coma continuo 0.314Ictus post operatorio 0.003Tamponamento cardiaco 0.0001Occlusione del trapianto 0.426Fibrillazione post operatoria 0.0001Ventilazione prolungata 0.0001Dialisi post operatoria 0.0001Insufficienza renale 0.0001

n60

4.4.3 Modelli joint

Questa sezione descrive brevemente i modelli di tipo joint, i quali per essere calibratirichiedono tecniche di integrazione numerica: in particolare tecniche di quadraturagaussiana. Entrambi i modelli sono applicati ai pazienti scompensati del databaseclinico del Uva Health System (Virginia).Liu et al.(2008)[71] costruiscono un modello usando le seguenti variabili:

• tij è il tempo dell'evento j del paziente i• tempo della morte del paziente i• ri(t) rischio di riospedalizzazione• wi

R è vettore di covariate per rischio di riospedalizzazione• Yij costo paziente (oppure tempo cura)• dNij(t) intervallo di rilevamento• zij è un vettore di covariate per il costo del paziente• λi(t) rischio di morte del paziente i• wi

D è vettore di covariate per rischio di morte

I termini di rumore sono: ui, vi, eij~N(0,σ2e).

Successivamente Liu (2009)[72] seguendo lo stesso approccio propone:

Uij=I(Yij>0) una variabile indicatore per il costo, Zij il vettore di covariate, i termini dirumore sono: ai~N(0,σ2

a), bi~N(0,σ2b), eij~N(0,σ2

e).

Esistono anche modelli joint che non necessitano di quadrature numeriche per esserecalibrati: Sun et al. (2012)[49] sviluppano un modello matematico che applicano a un

n61

repository di dati clinici del University of Virginia Health System: 1475 pazienti conetà tra I 60 e 89 anni che sono stati diagnosticati con scompenso cardiaco e trattati perla prima volta nel 2004 e seguiti fino al 31 luglio 2006. Le variabili del modello sono:

• responso Y(t)• un vettore di covariate • υ1 , υ2 sono 2 variabili latenti che sono stocasticamente indipendenti da • D il momento della morte• C il momento in cui avviene il censoring• T=min(C,D)• δ=I(D≤C)• una variabile N(t) che conta il numero di osservazioni fino al tempo t

Il modello joint proposto è:

Dove Λ0 è la funzione del hazard, α0(), μ0(), β0, γ0 sono 2 funzioni e 2 parametri da stimare con metodi di regressione. Y(t) in questo caso era il logaritmo dei dati di costo medici.

Autori Sun et al.Data pubblicazione 2012Età 0.0215Sesso 0.4663Etnia 0.0002

4.4.4 Statistica non parametrica

Gregory et al. (2005)[46] prendono in esame un campione di 82 pazienti valutati pertrapianto cardiaco del Tufts-New England Medical Center, un programmaaccademico per scompensati, tra il 2000 e 2001 che vengono seguiti per 2 annisuccessivi. Le osservazioni sono state ottenute integrando dati dal sistema dicontabilità interna, cartelle cliniche e database di mortalità del Social Security edomande fatte al gruppo di lavoro del programma.É stata fatta una analisi di statistica descrittiva per separare i pazienti che hannosubito il trapianto da quelli che non lo hanno subito, alle variabili discrete sono stati

n62

applicati il test chi-quadro di Pearson e il test esatto di Fisher, alle variabile continueè stato applicato il test di Student per 2 campioni (test paired), ai dati di mortalità èstato fatto un fit con la curva di Kaplan-Meier S(t). I tassi medi cumulati perospedalizzazioni inpatient (prescrizione medico) e per incontri outpatient (senzaprescrizione) sono stati stimati sulla base di tutto il campione, gli autori hanno anchestimato i tassi medi cumulati per pazienti che hanno subito trapianto e per quelli chenon lo hanno subito. I tassi medi cumulati per questi eventi sono stati stimati usandomodelli non parametrici di affidabilità right censored[47], sia n la dimensione delcampione, di il numero di pazienti morti nell' i-esimo intervallo, allora si stima lafunzione di ripartizione F(t)=Pr(T<t) con:

Tassi medi ottenuti con il modello matematico. Fonte: Gregory et al. (2005)[106]

Il fatturato per eventi inpatient è stato stimato usando i DRG Medicare e le spese perammissione sono state stimate con la metodologia Medicare 2003, questo fatturatocomprende sia gli incassi per l'educazione indiretta dei medici (IME) e outlier dicosto. I tassi basici di remunerazione sono stati stimati in base a quelli dei grandiospedali, corretti con un tasso salariale del 1.15 e un tasso geografico del 1.1. I costidiretti (inpatient e outpatient), che vengono usati nella economia sanitaria comeproxy per i costi marginali[48], sono stati calcolati usando i costi diretti dell'unitàinfermieristica, costi diretti clinici e costi diretti secondari. I costi sono statiattualizzati con il tasso di inflazione del Consumer Price Index. Sono stati usatimodelli affidabilistici per il fatturato e i costi diretti, corretti con le probabilità disopravvivenza per ogni evento temporale.

n63

Fatturato e costi cumulati. Fonte: Gregory et al. (2005)[106]

4.5 Articoli di riferimento

Siccome questo lavoro utilizza una logica di controllo di gestione, quindi per ilcontrollo della spesa10, i modelli di regressione appaiono quelli più adatti in un primomomento, infatti gli altri modelli non consentono di ipotizzare delle leve percontrollare la spesa, eccezione fatta per i modelli joint, che però rischiano di portare auna complicazione eccessiva per gli scopi di questa analisi. Sarà costruito un modelloin cui il responso sarà il costo, le covariate verranno individuate attraverso un'opportuna analisi.Gli articoli più rilevanti per la nostra analisi sono Dunlay et al. (2010) e Esposito etal. (2009): gli autori considerano esplicitamente il censoring, l'asimmetria statisticadei dati costo e l'utilizzo dei modelli lineari generalizzati, non sarà necessario peròusare il formalismo dei modelli a 2 parti, poichè che non avremo osservazioni nulle dicosto.

10 Il dipartimento PAC è il responsabile lombardo delle funzioni legate alla programmazione, all'acquisto e al controllo di prestazioni necessarie a soddisfare i bisogni di salute espressi dagli assisti dell’ASL. Le analisi di costo vengono effettuate dalla struttura S.s.d. Contratti Sanitari, Valutazione Economica e Statistica.

n64

Analisi dei dati

L'analisi dei dati è stata fatta utilizzando il linguaggio R.[54] Nella costruzione delmodello collegheremo gli aspetti matematici con quelli gestionali.

5.1 Analisi esplorativa

L'analisi esplorativa ha due principali obiettivi:

1. Descrivere sinteticamente la popolazione: questo è particolarmentesignificativo per i data-set attuali, i quali possono contenere anche migliaia diosservazioni e centinaia di variabili, e non sono quindi interpretabiliintuitivamente da un essere umano

2. Verificare la plausibilità delle ipotesi di modellazione, per esempio sullapesantezza della coda della distribuzione ipotizzata

Sono state necessarie circa 4867 linee di codice R e un'intensa attività di debuggingper stampare l'output dell'analisi: tabelle, grafici e indici.

5.1.1 Metodologia per l'analisi esplorativa

L'analisi esplorativa è stata realizzata secondo lo schema presente in Vercellis (2006)[3] e arricchita con alcuni risultati della letteratura sulle code pesanti/code spesse(“heavy tail”/ “fat tail”)[69]. Indichiamo, per sinteticità, il momento (campionario) di

ordine n con , gli indici calcolati dall'algoritmo per le variabili

continue sono:

• Media: , l'indice di posizionamento più usato, corrisponde al

baricentro della distribuzione di probabilità• Media troncata: la media calcolata rimuovendo il 5% degli eventi estremi, aiuta

a dare un 'idea di quanto impattano gli eventi estremi sui momenti• Mediana = , un indice di posizionamento che divide la distribuzione in

2 parti equiprobabili• Range = xmax – xmin,

, un indice di dispersione

• Deviazione Media Assoluta = , un indice di dispersione

• Deviazione Standard: , un altro indice di dispersione

n65

• Coefficiente variazione =

• Asimmetria: , se positiva la distribuzione è asimmetrica, altrimenti vale

l'inverso• Curtosi: , viene tipicamente usata per misurare quanto i dati siano

distanti dalla distribuzione normale, però siccome questa interpretazione puòrisultare fuorviante per distribuzioni asimmetriche[117], dato che non ci sonovariabili continue simmetriche usabili come covariate, questo indice verràtrascurato in questa analisi. Ciononostante, l'analisi dello spessore delle codenon verrà trascurata.

La letteratura sulle code spesse (pesanti), importanti in ambito economico, ha ormaiuna longa storia, possiamo citare come pionieri in questo argomento Cauchy, Pareto eLevy. Analizzare le code della distribuzione di probabilità è importante perché aiuta acapire l'impatto degli eventi estremi sul fenomeno sotto studio. Qui riportiamo, perscopi di sintesi, due classificazioni proposte per le distribuzioni di probabilità in basealle caratteristiche delle code. Entrambe le classificazioni tengono conto del seguente risultato: se una funzione ètale per cui f(cx)=Cf(x), dove c e C sono due costanti, allora si dice che f gode dellaproprietà di invarianza di scala. In particolare se Pr(X>x)=1-Pr(X<x)=Cx-α, questalegge di potenza (o anche per certi valori di C, distribuzione di Pareto) gode dellaproprietà di invarianza per scala. É noto che per α≤2 la varianza diverge, e per α≤1 lamedia diverge. Più in generale possiamo dire che hanno code paretiane funzioni deltipo Pr(X>x)=L(x)x-α, dove L(x) soddisfa:

La classificazione più vecchia è quella del matematico Benoit Mandelbrot[74]. SianoX' e X'' variabili aleatorie iid con fdp f(x), X=X'+X'', f2(x)=∫f(u)f(x-u)du è la fdp di Xottenuta come convoluzione11. Se X è nota allora, la distribuzione condizionata,denominata rapporto di porzionatura di breve periodo (short-run portioning ratio) è:

Intuitivamente se X' e X'' sono dello stesso ordine di grandezza, allora si dice che ilrapporto di porzionatura è equilibrato, altrimenti concentrato. Più tecnicamente si può

11 La convoluzione di due funzioni f e g è:

n66

studiare la moda di min(X',X'') usando log(f(x)): se log(f(x)) è concava allora ilrapporto di porzionatura è equilibrato, se log(f(x)) è convessa allora il rapporto diporzionatura è concentrato. La porzionatura di lungo periodo è data dai classiciteoremi di convergenza della statistica. Oltre a questi criteri Mandelbrot suggerisceanche lo studio della funzione di ripartizione inversa Pr(x)-1 e del fattore di scalaσ(q):

• Aleatorietà corretta lieve (Proper mild randomness): rapporto di porzionaturadi breve periodo equilibrato per N=2

• Aleatorietà borderline lieve (Borderline mild randomness): rapporto diporzionatura di breve periodo concentrato per N=2, però diventa equilibrato seN è più grande di un valore finito.

• Aleatorietà lenta con momenti delocalizzati finiti (Slow randomness with finitedelocalized moments): Pr(X>x)-1 cresce più rapidamente di |logx| però piùlentamente di |logx|1/w con w<1, oppure σ(q) cresce più velocemente di q peròpiù lentamente di q1/w.

• Aleatorietà lenta con momenti localizzati finiti (Slow randomness with finiteand localized moments): Pr(X>x)-1 cresce più rapidamente di ogni potenza |logx|1/w però più lentamente di exp(|logx|γ) con γ<1, oppure σ(q) cresce piùvelocemente di q però rimane finito.

• Aleatorietà pre-selvaggia (Pre-wild randomness): Pr(X>x)-1 cresce piùrapidamente exp(|logx|γ) con γ<1 però più lentamente di x-1/2, oppure σ(q)diventa infinito per q>2.

• Aleatorietà selvaggia (Wild randomness): E(X2) infinito, però E(Xq) finito perun valore di q>0

• Aleatorietà estrema (Extreme randomness): tutti i momenti infiniti

Vediamo di rappresentare le Pr(x)-1 proposte da Mandelbrot in un grafico:

n67

Pr(X>x)-1 per la classificazione di Mandelbrot

Evidentemente non possiamo rappresentare l'infinito intorno all'origine: per alcuni valori di w e γsarebbe necessario ridimensionare il grafico

Prima di continuare facciamo la seguente osservazione: parecchi concetti della teoriadella probabilità si possono riassumere usando il formalismo dello spazio Lp, lanorma-p è tradizionalmente definita sullo spazio vettoriale come:

Se X è un campione, ||X-μ||1 /n è la deviazione media assoluta campionaria e ||X-μ||2/√n è la deviazione standard campionaria non corretta.Mentre sullo spazio misurabile Y=(Ω, A, μ) (in questo caso spazio di probabilità) consideriamo l'insieme di funzioni misurabili f tale che:

Questa è appunto la definizione di spazio Lp.. La norma definita in questo spazio si

n68

indica con o con ||f||p. Notiamo subito che il fattore di scala di Mandelbrot altro non è che:

Finché σ(q) è finito e fX è definita su [0,+∞) la famiglia di funzioni xq forma unospazio Lp. Più in generale se:

Allora σ(q)=E[Xq]1/q < ∞.Un' altra classificazione più recente è una elaborazione dell'autore Nassim NicholasTaleb[69] del testo su eventi estremi di Embrechts et al.[73] usando le derivate diRadon-Nikodym ∂ν/∂μ12:

12 Alcuni teoremi della teoria della probabilità si possono ottenere usando la teoria della misura.Uno spazio di probabilità è la terna (Ω, A ,P[]), dove Ω è lo spazio campionario, A lo spazio deglieventi e P[] la misura di probabilità. Il teorema di Radon-Nikodym garantisce che sotto opportunecondizioni, data una f: Ω→[0,+∞), si possono assegnare 2 misure ν e μ tale per cui per ogni

sottoinsieme misurabile X di Ω: , dove f=dν/dμ è la derivata di Radon-Nikodym

n69

Fonte: Taleb (2014) [69]

Intuitivamente se otteniamo un buon fit di distribuzione da un campione, questo nonimpedisce che il campione osservato sia stato generato da una distribuzioneappartenente a un livello più alto del tableau.

Esempio: la lognormale è classificata nella sezione subesponenziale, però unadistribuzione a un livello più alto potrebbe aver generato una coda che sembra esserelognormale e ciò non è visibile dal campione per insufficienza di sampling nelle code.D'altra parte una coda spessa come quella della lognormale non può essere statagenerata da una distribuzione a un livello più basso, per esempio da una distribuzionegamma (gaussiana per approssimazione nel tableau).

Taleb riassume questa asimmetria nel principio logico: l'assenza dell'evidenza nonequivale all'evidenza dell'assenza, questa considerazione sarà utile più avanti nellacostruzione del modello. Riassumiamo brevemente le considerazioni di Taleb.Muovendoci sulla diagonale del tableau:

• Legge debole dei grandi numeri: siano X1,...,Xn n variabili aleatorie iid e sia lamedia campionaria , allora se la media μ è finita

n70

• Teorema centrale del limite (TCL): siano X1,...,Xn n variabili aleatorie iid conmedia μ e varianza σ2≠0 finite, allora

Teorema di Berry-Eseen: valgono le stesse ipotesi del TCL, in più se E|X|3 èfinito, allora

Dove Φ(x) è la funzione di ripartizione della N(0,1) e c una costante positiva.• Supporto compatto: la funzione di ripartizione Pr(X<x) è nulla al di fuori di un

insieme compatto• Condizione di Cramer: stabilisce la convergenza dei momenti esponenziali

e la possibilità di scrivere la funzione generatrice dei momenti.Sotto opportune condizioni anche le e-tx possono essere viste come membri di uno spazio Lp.

• Distribuzioni subesponenziali: distribuzioni introdotte da Chistiakov (1964)[75] che soddisfano il seguente limite

Però usate da Taleb nella forma di Teugels(1975): [76] siano Pr(X ≤ x)=F(x) lafunzione di ripartizione e Pr(X1+X2 < x)=F2(x)=∫f2(x)dx ottenuta comeconvoluzione13 , si considerino distribuzioni sub-esponenziali quelle che

13

n71

soddisfano il seguente limite:

Questo limite è equivalente al seguente:

Il quale stabilisce che il ruolo degli eventi estremi nella somma di tutti glieventi non è trascurabile.

• La norma più conosciuta è la deviazione media assoluta:

Questa è sempre finita purché la media μ lo sia anche. Per una distribuzionecon coda paretiana, se α≤1 la media diverge.

Per quanto riguarda l'utilizzo di grafici, l'uso congiunto di un boxplot (box-and-whisker) e un istogramma, 2 rappresentazioni classiche, è un buon compromesso trariuscire a rappresentare insieme il centro e la coda della distribuzione: l'istogrammarappresenta le frequenze empiriche nella loro proporzione naturale, quindi damaggior evidenza al centro della distribuzione, mentre essendo il boxplot costruitousando i quartili della distribuzione, da un lato i bordi (baffi) risentono poco dellapresenza di eventi estremi e dall' altro lato tutti gli eventi che cadono al di fuori deibaffi sono sufficientemente rari da essere considerati eventi estremi14. Gli istogrammipresentano una linea verticale sulla destra e/o sulla sinistra. Queste linee indicano ipercentili ricavabili dal teorema di Chebyshev, il quale stabilisce che dato un numeroγ ≥ 1 e un insieme di m valori, una percentuale pari ad almeno (1-1/ γ2) di questivalori si colloca all'interno dell'intervallo (μ - γ σ, μ + γ σ). Per effetto dicomparazione con altri lavori usiamo γ=6 (six sigma), quindi 1-1/ γ2 = 97,2% incondizioni di media e varianza finita. Questo consente di identificare rapidamentevalori particolarmente estremi e allo stesso tempo verificare se la media campionariae la deviazione campionaria sono simili a quelle della popolazione, altrimenti se 1-1/γ2 < 97,2% siamo motivati a concludere che ciò è falso.

14 Il termine outlier è poco preciso, una volta che la sua definizione matematica tende a variare a seconda del contesto, in questo caso verrebbero considerati come outlier tutti gli eventi che cadono al di fuori dell'intervallo [qL-3Dq, qU+3Dq], dove qL è il quartile inferiore, qU il quartile superiore e Dq= qU - qL

n72

Avendo quindi chiarificato lo schema concettuale di analisi che stiamo usando,illustriamo il metodo di analisi per la variabile numero di eventi per paziente deldatabase degli eventi.

####### Attributo V1 per codice identificativo univoco del paziente #######Numero osservazioni = 371766.00Media = 1.89Media Troncata (0.05) (coda destra) = 1.58Mediana = 1.00Range = 70.00Deviazione Media Assoluta = 1.32

n73

Deviazione Standard = 1.78Coefficiente Variazione = 0.94Asimmetria = 4.53

Notiamo innanzitutto che anche se nell'istogramma gli eventi estremi non si riesconoa vedere, il boxplot ci mostra che sono anche abbastanza numerosi. La grande partedegli eventi si colloca dentro all'upper bound di Chebyshev (μ + 6 σ =12.57), ealmeno 50% dei pazienti ha un solo evento. Il restante degli eventi è collocato entroun intervallo abbastanza ampio (70-12.57=57.43). La distribuzione è chiaramenteasimmetrica.Iniziamo l'analisi delle code col verificare la rapidità di convergenza della mediaÊ[X] e della varianza Ê[(X-μ)2]. Ora siccome in un campione non abbiamo mai a chefare con una serie infinita di numeri, è necessario per scopi pratici definire un erroredi stima ammesso:

1. Se usiamo un limite inferiore L per il valore stimato, allora commettiamo unerrore di stima al massimo pari a L-Ê[g(X)]

2. Se invece ammettiamo che ci possa essere divergenza, allora E(g(X))=∞

n74

In questo caso i grafici non forniscono evidenza di divergenza per la media e per lavarianza. Un modo per verificare scalabilità nelle code è studiare la mediacondizionata[69]:

Se la coda è paretiana abbiamo che:

Mentre per distribuzioni a coda sottile, come la gaussiana:

Indichiamo α*, il valore recuperato dal campione, per differenziarlo da α, ilparametro teorico della distribuzione. Sul grafico che segue la linea rossa corrispondeal limite pari ad 1 e la linea blu al limite α/(α-1)=1.5 per α=3. Ricordiamo dal tableaudi Taleb, che per α*>3, come in questo caso, abbiamo una convergenza non troppolenta per i limiti classici della statistica.

n75

Un'altra tecnica tiene conto del risultato seguente: se Pr(X>x)= Cx-α, è noto chepassando ai logaritmi si ha log(Pr(X>x)=logC-αlogx, il che fornisce un agile metodoper verificare la presenza di una legge di potenza: in un grafico log-log la presenza diuna retta è sinonimo di legge di potenza. Il primo grafico è per la densità mentre ilsecondo è per Pr(X>x)=1-Pr(X<x)

α*=4.27 α*=4.27

Abbiamo quindi che α≤α*=4.27. Il modo più diretto per verificare la subesponezialitàè il calcolo del limite:

n76

Come prima accennato F2(x) è la convoluzione di una variaibile aleatoria X con sestessa. La convoluzione è agevole da implementare su R, l'unico problema è che secalcolata sul campione per le frequenze sulla coda si arriva a un punto dove c'èinsufficienza di sampling e il limite si comporta come quello di una distribuzione concoda sottile. Vediamo nel primo grafico che il limite tende a convergere verso 2 primadel valore 40 sull'ascissa, suggerendo un'evidenza significativa per il caso sub-esponenziale. Confrontiamo con il grafico sulla destra che fa riferimento alla codadestra della variabile età del paziente (V3) : vediamo chiaramente che il limite disubesponenzialità cresce monotonamente e diverge, questo è tipico delle distribuzioninon subesponenziali (ossia superesponenziali).

Limite di subesponenzialità per la variabile numero di eventi per paziente (V1)

n77

Limite di subesponenzialità per la variabile numero età del paziente (V3)

In base a questi risultati la variabile eventi per paziente (V1) verrebe collocatasezione subesponenziale del tableau di Taleb.

5.1.2 Risultati dell'analisi esplorativa

Procediamo adesso alla descrizione sistematica del data-set. Il data-set ha un totale di701701 eventi distribuiti tra 371766 pazienti. L'estrazione ha riguardatoprincipalmente i pazienti residenti in Lombardia, in quanto la percentuale di quanticambiano residenza durante la durata dello studio è trascurabile.Vediamo, dal prossimo grafico, che il numero di pazienti scompensati che il sistemasanitario lombardo deve assorbire continua a crescere, questo è dovuto in parte alfatto che i pazienti che entrano nello studio dal 2000 e vi rimangono fino al 2012 ameno che non muoiano nel frattempo. La riduzione del tasso di crescita degli eventinegli ultimi anni dello studio suggerisce l'avvicinarsi del raggiungimentodell'equilibrio rispetto agli ingressi e alle uscite dei pazienti nella popolazione distudio.

n78

La media dell'età è 76 anni, la quale è poco distorta dagli eventi estremi e comeprima accennato la distribuzione può essere considerata super-esponenziale. Vediamodal prossimo grafico, come l'età media tende ad aumentare negli anni:

n79

Unità di misura: anni

####### Attributo V3 età #######Numero osservazioni = 701700.00Media = 76.09Media Troncata (0.05) = 76.85Mediana = 78.00Range = 110.00Deviazione Media Assoluta = 10.51Deviazione Standard = 11.91Coefficiente Variazione = 0.16Assimetria = -1.35

L'età media dei pazienti tende a crescere come vediamo dal prossimo grafico, infattipassa da 74.13 a 77.69 in 12 anni.

Unità di misura: anni

Il prossimo grafico mostra una piccola predominanza dei casi di scompenso maschilirispetto a quelli femminili, la quale oscilla tra il 52 e 53%.

n80

Il numero di eventi con un ricovero è il 99.99%, infatti qui l'analisi delle code è dipoco valore, però possiamo comunque notare che la media è poco affetta dagli eventiestremi.

Unita di misura: numero di ricoveri

####### Attributo V16 numero di ricoveri contenuti nell'evento #######Numero osservazioni = 701701.00

n81

Media = 1.16Media Troncata (0.05) (coda destra) = 1.07Mediana = 1.00Range = 21.00Deviazione Media Assoluta = 0.23Deviazione Standard = 0.50Coefficiente Variazione = 0.43Assimetria = 4.34

La prossima variabile è la durata di ospedalizzazione in giorni, la cui media è pari13.74 giorni, sulla quale l'impatto degli eventi estremi è di 2.48 giorni. Notiamo cheche ci sono eventi particolarmente estremi, arrivando a eventi di 1000 giorni didurata.

n82

####### Attributo V41 durata in giorni dell'evento #######Numero osservazioni = 701700.00Media = 13.74Media Troncata (0.05) (coda destra) = 11.26Mediana = 10.00Range = 1597.00Deviazione Media Assoluta = 9.99Deviazione Standard = 14.86Coefficiente Variazione = 1.08Assimetria = 6.30

Sia la media e la varianza convergono (la varianza più lentamente), e il test di sub-esponenzialità indica convergenza:

n83

n84

La media si comporta in maniera apparentemente non stabile, oscillando con il valorepiù basso nel 2003 (13 giorni) e quello più alto nel 2011 (14.5 giorni), periodo entro ilquale cresce ininterrottamente. Questo non è però necessariamente un cattivo segnale,visto che se usata come covariata in un modello, quello che conta è tutta ladistribuzione e non solo il primo momento, ovvero la media. Infatti mostreremo piùavanti come il modello di regressione sviluppato riesce a tenere conto di questeoscillazioni.

n85

Per quanto riguarda il problema del censoring, notiamo che la grande parte deipazienti sono deceduti entro la fine dello studio (64%) , ciononostante una parte nontrascurabile (34%) è stata troncata.

La mortalità ha un andamento temporale interessante, infatti il tasso di mortalità èdecrescente per i pazienti di questo studio, il tasso di mortalità intraospedalieromostra un trend di crescita fino al 2008, il quale viene successivamente invertito. Èda notare che il decesso intraospedaliero rappresenta solo il 9,4% degli eventi.

n86

Unità di misura: percentuale sul totale annuale

n87

La percentuale di eventi con riabilitazione a fine ricovero è maggiore rispetto allariabilitazione a inizio ricovero, anche se prevalgono i casi senza nessun ricovero.

n88

Vediamo che la percentuale di eventi con riabilitazione continua ad aumentare neglianni, tranne per l'ultimo anno:

n89

Da quattro grafici che seguono risulta chiaro che il numero di eventi con terapiaintensiva è più alta di quello per eventi con DRG chirurgico, che a sua volta è più altodegli eventi con cardiochirurgia, che è maggiore degli eventi in cui è stato necessarioimpiantare un defibrillatore.

n90

n91

Guardiamo adesso l'andamento di queste variabili nei diversi anni, abbiamo ingenerale una situazione di declino, tranne per l'uso dei defibrillatori che presentanoun chiaro trend di crescita.

n92

n93

I reparti medici più visitati sono, in ordine decrescente, la medicina generale,cardiologia, l'unità coronarica, l'allergologia e l'unità “recupero e riabilitazionefunzionale”. Sotto sono riportati i reparti visitati sia in ingresso e in uscita:

Disciplina del reparto di ingressoRepparto Percentuale di eventiMedicina generale 46.54%Cardiologia (a) 24.08%Unità coronarica (l) 8.64%Allergologia 4.18%Recupero e riabilitazione funzionale (g) 3.17%Geriatria 2.57%Astanteria 2.00%Terapia intensiva (i) 1.97%Nefrologia 1.28%Cardiochirurgia 1.08%Chirurgia generale 0.78%Neurologia (d) 0.71%Urologia 0.35%Chirurgia vascolare 0.35%Oncologia 0.32%Nefrologia (abilitazione trapianto rene) 0.26%Ortopedia e traumatologia 0.25%Lungodegenti 0.25%Neurochirurgia 0.14%Pensionanti 0.12%Malattie infettive e tropicali 0.12%Gastroenterologia 0.12%Otorinolaringoiatria 0.09%Pediatria (e) 0.08%

n94

Chirurgia toracica 0.07%Malattie endocrine, del ricambio e della nutrizione 0.06%Ostetricia e ginecologia 0.06%Medicina del lavoro 0.05%Oculistica 0.05%Neuro-riabilitazione 0.05%Angiologia 0.03%Ematologia 0.03%Reumatologia 0.03%Anatomia ed istologia patologica 0.02%Terapia intensiva neonatale 0.02%Dermatologia 0.02%Chirurgia maxillo facciale 0.01%Cure palliative/hospice 0.01%Cardiochirurgia pediatrica 0.01%Neonatologia 0.01%Riabilitazione di mantenimento 0.01%Chirurgia plastica 0.01%Radioterapia 0.00%Nido 0.00%Chirurgia pediatrica 0.00%Unità spinale 0.00%Neuropsichiatria infantile 0.00%Grandi ustioni 0.00%Oncoematologia 0.00%Emodialisi 0.00%Medicina nucleare 0.00%Nefrologia pediatrica 0.00%Oncoematologia pediatrica 0.00%Odontoiatria e stomatologia 0.00%

Disciplina del reparto di uscita

Repparto Percentuale di eventiMedicina generale 44.33%Cardiologia (a) 30.48%Recupero e riabilitazione funzionale (g) 9.12%Allergologia 4.27%Geriatria 2.90%Nefrologia 1.57%Unità coronarica (l) 1.51%Terapia intensiva (i) 1.00%Astanteria 0.78%Lungodegenti 0.62%

n95

Cardiochirurgia 0.59%Neurologia (d) 0.37%Nefrologia (abilitazione trapianto rene) 0.33%Chirurgia generale 0.32%Oncologia 0.31%Chirurgia vascolare 0.25%Urologia 0.23%Pensionanti 0.13%Neuro-riabilitazione 0.13%Malattie infettive e tropicali 0.12%Neurochirurgia 0.08%Gastroenterologia 0.08%Pediatria (e) 0.06%Malattie endocrine, del ricambio e della nutrizione 0.06%Chirurgia toracica 0.05%Cure palliative/hospice 0.05%Medicina del lavoro 0.04%Angiologia 0.03%Riabilitazione di mantenimento 0.03%Ortopedia e traumatologia 0.03%Reumatologia 0.02%Ematologia 0.02%Neonatologia 0.02%Cardiochirurgia pediatrica 0.02%Otorinolaringoiatria 0.02%Ostetricia e ginecologia 0.01%Anatomia ed istologia patologica 0.01%Unità spinale 0.00%Terapia intensiva neonatale 0.00%Dermatologia 0.00%Radioterapia 0.00%Chirurgia plastica 0.00%Chirurgia pediatrica 0.00%Neuropsichiatria infantile 0.00%Oculistica 0.00%Chirurgia maxillo facciale 0.00%Nido 0.00%80 0.00%82 0.00%Emodialisi 0.00%Nefrologia pediatrica 0.00%Day surgery (b) 0.00%Grandi ustioni 0.00%Neurochirurgia pediatrica 0.00%Odontoiatria e stomatologia 0.00%Oncoematologia 0.00%

Le comorbidità più presenti sono in ordine decrescente: lo scompenso cardiacocongestizio, l'ipertensione, l'aritmia cardiaca e le malattie polmonari croniche.

n96

comorbidità Percentuale di eventiscompenso cardiaco congestizio 31.52%aritmia cardiaca 15.69%ipertensione 13.96%malattie polmonare croniche 10.66%insufficienza renale 7.23%arteriopatia obliterante periferica 4.06%disturbi della circolazione polmonar 3.07%complicazioni del diabete 2.77%deficienze legate all'anemia 2.67%tumore 2.26%malattia del fegato 1.64%demenza 1.24%disturbi legati ai fluidi e elettroliti 1.20%emiplegia 0.82%cancro metastatico 0.48%coagulopatia 0.25%psicosi 0.25%perdita di peso 0.12%abuso di alcool 0.10%

Notiamo che non c'è un netta dominanza di una ASL nella distribuzione degli eventi:i codici più numerosi corrispondo rispettivamente alla città di Milano, Brescia eVarese. Le eventi di cura sono sostenuti principalmente in Aziende Ospedaliere enelle case di cura private accreditate.

Struttura del primo ricovero dell'evento_V65 Percentual di eventiA.O. 68.73%ASL 2.23%CASA DI CURA PRIVATA ACCREDITATA 14.60%CASA DI CURA PRIVATA NON ACCREDITATA 0.07%IRCCS PRIVATO 6.89%IRCCS PUBBLICO 3.89%OSPEDALE CLASSIFICATO EQUIPARATO 3.59%

n97

Descrizione ASL della prima struttura dell'evento_V67 Percentuale di eventiA.S.L. DELLA CITTA' DI MILANO 11.52%A.S.L. DELLA PROVINCIA DI BERGAMO 9.11%A.S.L. DELLA PROVINCIA DI BRESCIA 14.05%A.S.L. DELLA PROVINCIA DI COMO 4.55%A.S.L. DELLA PROVINCIA DI CREMONA 5.06%A.S.L. DELLA PROVINCIA DI LECCO 3.45%A.S.L. DELLA PROVINCIA DI LODI 2.15%A.S.L. DELLA PROVINCIA DI MANTOVA 5.21%A.S.L. DELLA PROVINCIA DI MILANO 1 6.84%A.S.L. DELLA PROVINCIA DI MILANO 2 4.74%A.S.L. DELLA PROVINCIA DI MILANO 3 5.28%A.S.L. DELLA PROVINCIA DI MONZA E BRIANZA 1.96%A.S.L. DELLA PROVINCIA DI PAVIA 7.25%A.S.L. DELLA PROVINCIA DI SONDRIO 1.59%A.S.L. DELLA PROVINCIA DI VARESE 9.74%A.S.L. DI MILANO 6.44%A.S.L. DI VALLECAMONICA - SEBINO 0.98%A.S.L. DI VALLECAMONICA-SEBINO 0.08%

Sotto abbiamo le percentuali di eventi per i gruppi definiti nella sezione 3.2. Abbiamoche il 70% appartiene al gruppo 1, il 17% al gruppo 2, il 13% al gruppo 4 e il restoappartiene al gruppo 4.

Gli eventi MDC5 sono quasi l'80% degli eventi, gli eventi per sindrome di shock,angioplastica coronarica e bypass aortocoronarico sono una percentuale molto piccoladel data-set (<3%), ciononostante è interessante osservare come queste ultime sicomportano nel tempo, abbiamo una situazione di diminuzione, tranne perl'angioplastica coronarica:

n98

n99

La variabile ricoveri valorizzati è la variabile più importante del data-set:

####### Attributo V29 somma delle valorizzazioni in euro dei ricoveri appartenenti all'evento#######Numero osservazioni available= 699252.00Numero di osservazioni NA = 2449.00Percentuale di osservazioni NA = 0.00Media = 6030.94Media Troncata (0.05) (coda destra) = 4609.21Mediana = 3380.00Range = 498581.99

n100

Deviazione Media Assoluta = 4271.29Deviazione Standard = 8104.00Coefficiente Variazione = 1.34Assimetria = 6.47

I primi 2 anni del dataset corrispondo a un diverso regime di DRG, come possiamovedere dal prossimo grafico. É da notare che anche dopo il taglio del valore dei DRGavvenuto nel 2002, il costo medio continua comunque ad aumentare negli anni.

È necessario quindi correggere questa variabile, rimuovendo gli anni 2000 e 2001.Rimuovendo i primi 2 anni, infatti il range si riduce della metà, rimanendo comunquealto a 250.000 €. Per la variabile corretta la media dei ricoveri è di 5688 €, con unadifferenza sulla media di 1312 euro quando rimuoviamo il 5% degli eventi estremi.Del resto abbiamo per questa variabile una situazione simile a quella della durata diospedalizzazione.

n101

####### Attributo V29 somma delle valorizzazioni in euro dei ricoveri appartenentiall'evento_corretto #######Numero osservazioni = 701701.00Media = 5688.66Media Troncata (0.05) (coda destra) = 4376.80Mediana = 3281.00Range = 256248.99Deviazione Media Assoluta = 3891.33Deviazione Standard = 7327.80Coefficiente Variazione = 1.29Assimetria = 5.23

n102

Sia la media e la varianza convergono (la varianza più lentamente):

Siccome l'analisi delle code per questa variabile è particolarmente importante, questaverrà ripresa nel punto 5.2.5, qui ci limitiamo a notare che α*>3 e che non possiamoescludere il caso sub-esponenziale, sulla base dei seguenti grafici:

n103

Oltre al costo del ricovero, sono presenti nel data-set anche i costi delle protesi,variabile per la quale rimuovere i primi 2 anni è indifferente. Qui notiamo che ladistribuzione è chiaramente bimodale (che in questo caso spiega una mediana più altarispetto alla situazione precedente), diversamente da quanto avveniva per la variabileprecedente. Procedendo con il confronto notiamo che l'impatto sulla media del 5%degli eventi estremi è pari a 334 euro, quindi con una differenza di 978 € rispetto aprima.

n104

####### Attributo V30 somma valuta dei ricoveri appartenenti all'evento_corretto #######Numero osservazioni = 701701.00Media = 2836.27Media Troncata (0.05) (coda destra) = 2502.36Mediana = 3586.00Range = 39453.00Deviazione Media Assoluta = 2383.62Deviazione Standard = 2764.38Coefficiente Variazione = 0.97Assimetria = 3.10

La media e la varianza convergono:

n105

Notiamo che c'è un aumento costante della media dal 2003 fino al 2011, per unadifferenza di 1575 €:

n106

Dal prossimo che vediamo che il limite di subesponenzialità è soddisfatto:

5.2 Costruzione del modello

Come abbiamo visto prima i modelli di regressione dipendono fortementedall'assunzione che si fa rispetto alla distribuzione del responso e più in particolaredei residui. Quindi è buona norma scegliere una distribuzione adatta a descrivere idati di costo.Per realizzare un fit di distribuzione esistono diversi criteri in competizione nella

n107

letteratura statistica (statistica test χ2 di Pearson, statistica di Kolmogorov, criterio diCramer-Von Mises, etc.), tipicamente questi possono essere espressi attraverso ilmetodo comune del calcolo del p-value. Il p-value è un valore probabilistico, unnumero tra zero e uno, e per quanto riguarda il test di ipotesi è la probabilità dirifiutare l'ipotesi nulla H0 quando l'ipotesi nulla H0 è vera (probabilità dell'errore delprimo tipo), quindi un valore del p-value prossimo a zero ci motiva fortemente arifiutare H0 . Questa analisi sarà basata sull'ormai classico test di Kolmogorov-Smirnov[79][2] cheusa come statistica test la differenza massima tra la funzione di ripartizione teoricaproposta e la funzione di ripartizione empirica calcolata sul campione:

In termini probabilistici se il valore stimato della statistica test d*n >λγ/√n , con 1-

γ=ampiezza dell'errore di tipo I, λγ il γ-quantile della distribuzione di Kolmogorov e nla dimensione del campione. Il nostro procedimento sarà quello di ottenere il fit di una/o più distribuzioni scelte usando il metodo classico della massima verosimiglianza e poi usare il test di Kolmogorov-Smirnov come criterio per valutare il “goodness-of-fit”. Facciamo notare insieme a Taleb[82][69] che per le decisioni economiche è di interesse oltre all'errore probabilistico, anche l'errore monetario di previsione, che qui è calcolato con il seguente metodo:

• Si realizza un fit per i dati di costo dal 2002 al 2007• Il costo atteso per ogni evento è E[X]=μ• Per il 2008 possiamo stimare il costo atteso con neventi μ, siccome in un primo

momento non è l'interesse principale valutare la stocasticità del numero di eventi annuali, neventi lo otteniamo come il numero effettivo di eventi verificatosinel 2008

• Il costo effettivamente verificatosi per il 2008 è dato dalla somma degli eventi di costo Σxi

• La differenza Δ=Σxi- neventi μ è l'errore monetario annuale assoluto di previsione• Il rapporto Δ/ Σxi ci da lo stesso errore però in termini percentuali

Il contenimento di questo errore è una delle principali priorità del ciclo di pianificazione e controllo tradizionale.[83][84] Qui in realtà abbiamo a che fare con un ciclo di controllo non aziendale ma regionale, però con le dovute precauzioni alcuni risultati della letteratura possono essere estesi alla nostro problema.

La parola outlier non ha una chiara definizione matematica e il significato varia

n108

spesso a seconda del contesto. In esperimenti scientifici tradizionali viene spesso fatto riferimento alla distribuzione gaussiana e si parla di outlier per osservazioni che sono distanti dalla media un certo numero σ di volte, oppure nel boxplot si usano i quartili per definire gli estremi del diagramma a baffi. Finché abbiamo a che fare con distribuzioni a coda sottile (aleatorietà lieve/ gaussiana per approssimazione) questo modo di procedere non presenta molti problemi, però per distribuzioni con code più spesse si finisce per scartare sezioni anche significative della coda.Il metodo di analisi qui proposto è il seguente:

• Usando la rappresentazione di funzione a tratti, il campione viene separato in 2zone A e B, la zona B contiene gli outlier rispetto alla zona A

• Il fit di distribuzione viene ottenuto per la zona A e il goodness-of-fit viene giudicato con la statistica di Kolmogorov-Smirnov

• 2 fit di distribuzione diversa comparati con la statistica di Kolmogorov-Smirnov stabiliscono in genere bound diversi per A e B

• Avendo ottenuto separatamente Pr(X<x|A) e Pr(X<x|B), siccome per costruzione Pr(A)+Pr(B)=1, allora possiamo ricostruire la funzione di ripartizione Pr(X<x) usando il classico teorema delle probabilità totali

Questo metodo di analisi rivela la sua vera potenza implementato in un linguaggio diprogrammazione, una volta che permette un'analisi di tipo iterativo che permette distudiare come l'assenza o la presenza di certi valori del campione impattano sullaperformance di un modello. Siccome ci concentreremo prima sullo studio degli eventidel dataset, avremo la situazione semplificata perché i valori mancanti sonotrascurabili (0,34%) e non ci sono valori nulli di costo. Vediamo alcuni fit didistribuzione ottenuti per la variabile V29 (valorizzazione in euro dei ricoveri), in blula funzione di ripartizione proposta:

n109

n110

n111

Quindi per esempio per la distribuzione gamma avremo che:

Nel caso di questo metodo ∫Agamma(x;r,λ)dx=1 però più in generale si può realizzare un fit il cui integrale non sia pari ad 1, questa normalizzazione provvede a rettificare il problema.Alcune osservazioni importanti:

• Notiamo che nessun fit è in grado di generare il valore standard del p-value >0.5, il motivo è la presenza dei valori ripetuti (ties nel linguaggio statistico) che fanno in modo che la funzione di ripartizione empirica sia discontinua. Il problema non è facilmente rimediabile con una correzione di continuità, infatti se prendiamo un valore tie particolarmente elevato e i 2 valori di campione più vicini a e b sulla sinistra e sulla destra e sostituiamo i valori presenti nell'intervallo [a,b] con una distribuzione omogenea di valori, la frequenza di questi valori rispetto al resto del campione è cosi elevata che anche usando questo metodo per tutta la funzione di ripartizione la statistica di Kolmogorov-Smirnov rimane praticamente inalterata.

• Questo esercizio illustra anche un'altra problematica dell'uso dei p-value/statistica test: per uno stesso campione possiamo ottenere fit di distribuzione diverse equivalenti in termine di errore. Vediamo per esempio che il fit con la distribuzione gamma e quello con la distribuzione lognormale sono equivalenti, senza che ci sia un motivo a priori per scegliere tra una o l'altra. Il fit con la distribuzione normale inversa gaussiana (NIG) è preferibile rispetto a quello con la distribuzione gamma, una volta che riuscendo a descrivere più dati con lo stesso errore di fit, risulta un fit più potente. La distribuzione alfa-stabile di Levy genera un errore di fit comparabile agli altri, con il problema però che se calibrata su questi valori la media risulta non definita.

• In realtà in questo caso specifico abbiamo un elemento ulteriore con il quale giudicare la bontà del modello, ovvero confrontare gli errori previsionali monetari come prima annunciato:

n112

Tabella con errore monetarioannuale di previsione

Δ (errore percentuale) Δ(Y A) (errore percentuale)ϵ

gamma 159.147.634 € (49%) 7.067.650 € (6%)

lognormale 159.132.022 € (49%) 7.057.495 € (6%)

NIG 127.861.614 € (39%) 13.203.926 € (7%)

alfa-stabile Non definito Non definito

Zoom sul fit con la distribuzione gamma, la differenza massima si ha intorno al valore 2650.

In base a queste considerazioni è stata scelta la distribuzione normale inversa gaussiana per l'intervallo [0,12000€]. La NIG è un caso particolare della distribuzione iperbolica generalizzata introdotta da Barndorff-Nielsen (1977)[80], perle proprietà sia della NIG e dell'iperbolica generalizzata si veda Barndorff-Nielsen (2013)[81]. La fdp della NIG è data da:

n113

Mentre la fdp dell'iperbolica generalizzata è data:

α è un parametro per la pesantezza delle code, β una parametro di asimmetria, δ un parametro di scala, μ un parametro di centralità, γ=√α2-β2

. λ è un parametro specifico della distribuzione iperbolica generalizzata, per λ=-1/2 abbiamo la NIG. Kλ è una funzione di Bessel del secondo tipo modificata15, che per scopi computazionali può essere rappresentata come[85][86]:

Ricapitolando abbiamo come densità per l'intervallo A:

In questo caso ancora una volta ∫ANIG(y;α,β,μ,δ,γ)dx=1. Useremo Y come variabile d'ora in avanti in ottica della costruzione del modello di regressione. La calibrazione della NIG ottenuta su R è , μ=2897, σ=1866, γ=57816.

L'algoritmo per realizzare il fit di distribuzione e l'analisi di regressione corrispondono a circa 2334 linee di codice in R con un'attività di debugging ancora più intensa rispetto all'analisi esplorativa.Nella costruzione di un modello di regressione è necessario decidere quali variabili

15 La funzione modificata di Bessel del terzo tipo è la stessa funzione però con un nome diverso. É un problema di carattere puramente storico e non matematico.16 Nella libreria ghyp di R la parametrizzazione standard per la NIG è diversa da quella tradizionale

n114

usare come covariate:• Le candidate a covariate sono state scelte in base alla revisione della letteratura• Queste variabili candidate sono o numeriche o binarie• Sotto è riportata la matrice di correlazione per le variabili numeriche con

evidenziate in ciano le variabili con correlazione rispetto al responso ritenuta non trascurabile, facendo però notare che la variabile V30 (protesi valorizzate) è una variabile di costo e quindi più adatta a servire come responso insieme alla variabile V29, tra le variabili V16 (numero di ricoveri per evento) e V41(durata dell'evento) sussiste una correlazione non trascurabile (come è ragionevole aspettarsi) evidenziata in giallo e quindi come prima spiegato è sconsigliabile l'uso simultaneo di entrambe come covariate.

Matrice di correlazione

V29 V30 V3 V16 V41 V64

V29 1 0.3650 -0.09 0.3785 0.577 0.094

V30 0.365 1 -0.036 -0.005 0.041 0.108

V3 -0.09 -0.03 1 0.0647 0.087 0.166

V16 0.378 -0.0051 0.065 1 0.682 0.085

V41 0.577 0.041 0.087 0.6825 1 0.157

V64 0.094 0.1076 0.166 0.0852 0.157 1

La filosofia di modellazione qui adoperata è inspirata dall'approccio di Mandelbrot:

• Cercare di costruire il modello più semplice possibile che emuli il fenomeno sotto studio con la maggior efficacia possibile, questo quindi diventa il benchmark per i miglioramenti successivi

• Inserire i miglioramenti uno alla volta in modo da valutare quali sono le modifiche più significative

In questo modo si ottiene un buon compromesso tra l'economicità del modello da un lato e l'efficacia rappresentativa dall'altro. Per efficacia rappresentativa qui si intende la capacità di contenere gli errori di modellazione, che sono comunque inevitabili visto che un modello è sempre e comunque una rappresentazione semplificata della realtà. Il tema degli errori di modellazione è un tema vasto che richiederebbe un livello di approfondimento che esula dagli scopi di questo lavoro. Alcuni dei problemi connessi a questo tema verranno illustrati in concomitanza al nostro esercizio di modellazione.

n115

5.2.1 Regressione con responso modellizzato tramite distribuzione gamma

Illustriamo la problematica dell'analisi di regressione prima con la distribuzione gamma e poi con la distribuzione NIG. Come già spiegato prima se un responso non gaussiano appartiene alla famiglia esponenziale, il metodo preferenziale è una regressione con modello lineare generalizzato. Su R il fit è ottenuto con la funzione glm(), la quale ottiene il fit usando il metodo dei minimi quadrati con ripesamento iterativo.Il primo fit illustrativo sarà ottenuto con la funzione di link più semplice immaginabile, ovvero una funzione lineare tra responso (V29=ricoveri valorizzati) e covariata (V41=durata di ospedalizzazione). La durata di ospedalizzazione è stata scelta come variabile prioritaria per la maggior correlazione con il responso e per il suo valore storico come proxy per l'assorbimento di risorse nella letteratura dei DRG.La funzione di regressione è per definizione[2]:

Quindi nel nostro caso un glm con funzione di link lineare e responso distribuito secondo una gamma:

La gamma è una distribuzione con 2 parametri r e λ, effettuando la regressione la media della distribuzione diventa funzione della(e) covariata(e), quindi se E(Y)=μ=r/ λ, allora r= λμ. Questo possiamo vederlo in una forma più famigliare:

Dove ε è il rumore17 previsionale, si confronti questa rappresentazione con quella che abbiamo chiamato regressione lineare multivariata nella sezione 2.4.2: li avevamo un rumore ε~N(0,σ2) simmetrico additivo, qui abbiamo un rumore sempre dipendente daun solo parametro, però asimmetrico e moltiplicativo. La media del rumore è 1 e la varianza è 1/λ2.Avendo quindi chiarificato la struttura del modello vediamo di valutarne la bontà o efficacia rappresentativa. La funzione glm() di R fornisce di default già alcune indicazioni: la significatività dei coefficienti e il criterio di informazione di Akaike (AIC).

17 È diventato usuale chiamare ε errore (terminologia ereditata dalla fisica), però qui facciamo notare che ε è una variabile casuale la cui distribuzione è un output dell'analisi di regressione e non un errore nel senso di errore di modellazione prima accennato

n116

Call:glm(formula = var ~ V41_x, family = Gamma(link = "identity"))

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2825.2069 1.6082 1756.8 <2e-16 ***V41_x 9.8190 0.1409 69.7 <2e-16 ***---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Gamma family taken to be 0.01346338)

Null deviance: 2496.3 on 176441 degrees of freedomResidual deviance: 2426.7 on 176440 degrees of freedomAIC: 2558630

Number of Fisher Scoring iterations: 4

La significatività Pr(>|t|) dei coefficienti è il p-value ottenuto dalla statistica test t-value. Osserviamo che l'interpretazione corretta di questo numero dipende da quale è l'ipotesi nulla H0 usata nel test: Piazza(2007) e Vercellis (2006) [2] usano H0 : βJ=0, un valore alto del t-value corrisponde quindi a un valore basso di Pr(>|t|), il che implica un coefficiente probabilmente diverso da zero. Quindi diversamente dal test di Kolmogorov-Smirnov un p-value basso in questo caso è buon segnale. Ricordiamo, come è consueto nella letteratura statistica, che un coefficiente significativo dal punto di vista probabilistico non necessariamente contribuisce con grandi variazioni nel responso.L' AIC è stato proposto da Akaike (1974)[87]:

Questo criterio è un compromesso tra il goodness-of-fit (funzione di massimaverosimiglianza) e l'overffiting, in questo caso rappresentato da un termine di penalità(2*numero di parametri calibrati del modello).Questi sono criteri probabilistici, l'errore monetario di previsione per il 2008 èΔ(2000≤Y≤4000)=7.114.710 € (6,1%), dell'errore Δ non condizionato ne parleremopiù avanti nell'analisi di rischio. Nei grafici che seguono possiamo avere un'idea dicome opera il modello: c'è una retta di regressione intorno a cui si dispongono leosservazioni, più grande è la differenza tra l'osservazione e la retta di regressione (iresidui) e più il colore tende dal rosso verso il verde: i grafici nella parte superiorefanno riferimento al rumore relativo ε prima esposto, i grafici nella parte inferiorerappresentano il rumore invece in termini assoluti.

n117

Facciamo notare che i grafici sulla sinistra non sono dei grafici di test, ma sono solo una rappresentazione alternativa al fit di regressione, mentre i grafici sulla destra costituiscono dei veri test una volta che il modello viene usato fuori dall'insieme di training. Dal punto di vista del controllo di gestione alla fine del 2007 abbiamo 5 annidi storia sui quali calibrare il modello, e dobbiamo fare una previsione per il 2008: una previsione evidentemente per il costo aggregato tra i valori compresi tra 2000 e

n118

4000 euro (questa è ottenibile senza il modello di regressione), poi abbiamo una seconda affermazione logica, più forte della prima, che consiste nell'affermare che nel2008 la media dei ricoveri valorizzati tenderà ad aumentare con l'aumentare della durata di ospedalizzazione cosi come è successo negli anni precedenti seguendo una curva di regressione stimata.Facendo delle prove su R con altre funzioni di link si vede che una soluzione quasi-ottima per il criterio AIC è data una funzione quadratica.

5.2.2 Regressione con responso modellizzato tramite distribuzione NormaleInversa Gaussiana (NIG)

Vediamo adesso di migliorare il modello usando la distribuzione NIG al posto delladistribuzione gamma. Qui ottenere la funzione di regressione risulta menoimmediato; usare un modello lineare generalizzato presuppone infatti la possibilità discrivere il responso come membro della famiglia esponenziale:

Supponendo anche che ciò sia possibile, resta poi il problema di avere un algoritmo su R capace di ottenere il fit con un responso distribuito come la NIG. Dopo una vasta ricerca su Internet si può concludere che probabilmente non è mai stata scritta una funzione per fare questo. Pertanto la funzione di regressione è stata ottenuta ricorrendo alla definizione:

Vediamo sotto il grafico di un fit di iperbolica generalizzata bivariata per le 2 variabili: questa è la funzione di densità congiunta fXY(x,y), la funzione condizionata fX|Y(y|x)=fXY(x,y)/fX(x) è più difficile da immaginare, però essa si dispone intorno alla funzione di regressione, per più dettagli su questo punto si veda Piazza (2007)[2]. La funzione di regressione è ottenuta per via numerica con un ciclo for, una rappresentazione analitica molto buona può essere ottenuta con una b-spline (curva inrosso)

n119

Curiosamente l'errore monetario di previsione Δ(0≤Y≤12000)=10.930.185 € (5,7%) èaddirittura minore rispetto al modello precedente. Seguono gli stessi graficiillustrativi di prima (notiamo che il dominio della covariata è stato ampliato insieme aquello del responso)

n120

Ora facciamo un'ulteriore test particolarmente istruttivo, siccome il modello affermache la media dei ricoveri valorizzati dipende positivamente dalla durata diospedalizzazione, il modello dovrebbe essere in grado di riprodurre il seguentefenomeno che emerge dall'analisi esplorativa del dataset.

Media della durata di ospedalizzazione neidiversi anni

Media dei ricoveri valorizzati nei diversi anni(2000 e 2001 sono esclusi da questa analisi)

Adesso estendiamo quindi le previsioni del nostro modello fino al 2012, ricordandoci di rimanere dentro il dominio [0,100]X[0,12000] entro cui è stato ricavato il modello:

n121

Vediamo che il modello riesce a riprodurre l'aumento del costo medio negli anni, peròsi nota chiaramente la presenza di un bias: il modello prevede più accuratamentedentro l'insieme di training e l'errore previsionale tende ad aumentare col passare deltempo.

5.2.3 Regressione (NIG) con insieme di training mobile

Esistono diversi modi per rimediare al problema del bias, qui proponiamo un modellodinamico ispirato dal ciclo controllo classico, detto anche ciclo cibernetico dicontrollo[83][84]:

n122

In generale possiamo possiamo dire che all'oggetto del controllo è associata una variabile di output O, mentre la fase di programmazione/previsione elabora un modello O=f(I,A) usando delle variabili input I e delle variabili ambientali A. La misura dei risultati a consuntivo prevede la misurazione di O', I' e A', in generale diversi dai valori previsti. Questo ciclo cibernetico in realtà comprende alcune dimensioni che non sono presenti in questa analisi:

• non stiamo valutando la possibilità di intervenire con delle azioni correttive sulpiano operativo concreto ma solo a livello modellistico, quindi più tecnicamente non stiamo modelizzando il double loop ma solo un single loop

• qui non stiamo prevedendo il valore delle variabili I e A, in realtà abbiamo una sola variabile, la durata di ospedalizzazione, la quale e conosciuta a posteriori, quindi I=I' e A=A': segue necessariamente che se O≠O' allora f≠f'

• I meccanismi di apprendimento immaginabili sono i più diversi, qui useremo quello più semplice: un insieme di training mobile con durata annuale fissa, che cercherà di evitare un overfitting del modello rispetto al passato.

Vediamo alcune prove:

n123

5 anni di training

n124

1 anno di training

n125

3 anni di training

n126

5.2.4 Calcolo dei p-value

La metodologia di selezione delle variabili binarie che useremo qui può essereclassificata come un metodo di tipo filter[3]: le variabili poco significative vengonoeliminate prima della costruzione del modello di regressione multivariato, l'approccioè simile al metodo di Hosmer e Lemeshow[31]. Consideriamo un generico modello diregressione:

Il calcolo del p-value è stato reso popolare dalla metodologia di Ronald Fisher[99], laquale consiste nei seguenti passi:

1) Identificare una ipotesi nulla H0 : nel caso della regressione di solito si pone bi=0e si cerca di rifiutare questa ipotesi2) Scegliere uno stimatore del parametro incognito bi

3) Calcolare la probabilità Pr(rifiutare H0|H0 è vera) = , cioè

il p-value (Rc=Regione Critica sono i valori di che ci portano a rifiutare l'ipotesinulla)

Lo stimatore secondo il metodo dei minimi quadrati è , che può essere

semplificato per le variabili dummy in:

Dove A è l'insieme di indici per le dummy=1 e B è l'insieme di indici per ledummy=0. Segue necessariamente che:

Cioè lo stimatore è la differenza delle medie condizionali. Per il caso multivariatoabbiamo invece che .Quindi per il problema della regressione abbiamo che se Y è il responso,

è il vettore di parametri incogniti da stimare:

n127

Dopodiché per ricavare il p-value è necessario calcolare:

Se f è la distribuzione normale, allora g è anche esso normalmente distribuito el'espressione analitica di g non è particolarmente complicata, se invece abbiamo a chefare con altre distribuzioni, come la NIG, allora è possibile ricavare g tramite ilmetodo di bootstrapping, proposto da Bradley Effron[100]. Il bootstrap è un metododi resampling a partire dal campione originale, quindi nel nostro caso se sonole nostre osservazioni possiamo ottenere un resampling del campione originale

con kj=k1,....,kn . Dal punto di vista computazione è necessario quindieffettuare B campionamenti casuali con ripetizione degli indici (k1,....,kn ) delleosservazioni, in questo modo possiamo ottenere una forma approssimata per g:

L'approssimazione è tanto migliore quanto più grande sono n e B. A questo puntoperò non abbiamo ancora la distribuzione per bi=0, ma questa può essere ottenutaattraverso una operazione di traslazione:

Quindi per calcolare il p-value:

Dove è il valore assoluto della statistica calcolato sul campione originale .Sotto riportiamo il calcolo del bi (dimensione dell'effetto) e del p-value calcolatiindividualmente per tutte le variabili usando l'insieme di training:

n128

Quindi, per esempio, siccome per la variabile V21 =6493 cade fuori dall'intervallo(-2000,1500) abbiamo che , quindi è altamente improbabileosservare la statistica test sotto l'ipotesi nulla.

n129

Da questa analisi preliminare possiamo già scartare alcune variabili con pvalue>0.05, seripetiamo l'analisi con le variabili rimanenti con i dati separati in diversi anni ecalcoliamo i p-value pvalue(2002), pvalue(2003), …, pvalue(2011), pvalue(2012) abbiamo untest più stringente. Per esempio per la variabile V19 abbiamo una situazione di questotipo:

n130

Siccome però stiamo facendo comparazioni di multiple di ipotesi H t simultaneamenteaumentiamo il rischio di osservare un risultato significativo per puro caso. Unapossibilità abbastanza rigorosa è quella di ricorrere alla disuguaglianza di Bonferroni,la quale garantisce che:

Il lato sinistro della disuguaglianza fornisce quindi un lower bound per la probabilitàdi ∩Ht

n131

Avendo quindi fatto il calcolo con un criterio più stringente vediamo che dobbiamoeliminare la variabile V21 dall'analisi. Se inseriamo le variabili che hanno passato ilfiltro in un modello di regressione multi-variato abbiamo che:

n132

Due variabili risultano ancora poco significative, se vengono eliminate, il modellorimanente per le variabili binarie è:

Variabili Effetti Valori_pb_0 Intercetta 3363.7554225571 0b_1 V9 -2481.4105890059 0b_2 V18 4342.3938642392 0b_3 V19 974.2301072898 0b_4 V20 190.320726345 0.01b_5 V22 4230.9906532779 0b_6 V33 -92.7825495301 0b_7 V55 109.816275999 0b_8 V61 111.8483361829 0b_9 V63 31.3164292649 0b_10 V86 -534.5174610578 0

Uniamo questi risultati a quelli ottenuti nell'analisi della variabile durata diospedalizzazione in un unico modello:

dove f() è la forma funzionale ottenuta nella sezione 5.2.2.

n133

Variabili Effetti Valori_pb_0 Intercetta 652.6653911951 0b_1 V9 -1593.0537422637 0b_2 V18 2689.7161877192 0b_3 V19 856.295957588 0b_4 V20 -89.5975343975 0.78b_5 V22 3824.7113938856 0b_6 V33 -29.9925279557 0b_7 f(V41) 0.8152413392 0b_8 V55 5.5776711406 0.06b_9 V61 -21.7924518722 1b_10 V63 -6.6270982894 0.34b_11 V86 -332.8774311965 0

Siccome la variabile f(V41) è particolarmente significativa e importante, abbiamo chealcune variabili binarie si rivelano superflue nel modello finale. È bene tenerepresente che non tutte le variabili rimosse sono inutili dal punto di vista del controllodella spesa, ma siccome nella logica previsionale quello che conta sono solo levariabili che generano variazioni significative abbiamo come modello di regressionefinale:

Variabili Descrizione Effetti Valori_pb_0 Intercetta Termine costante 662.27315704 0b_1 V9 Inizio evento con riabilitazione -1540.575717 0b_2 V18 Fine evento con riabilitazione 2662.9551861 0b_3 V19 Paziente è stato in terapia intensiva nell'evento 837.88766037 0b_4 V22 Paziente ha un DRG chirugico nell'evento 3839.3251122 0b_5 V33 Sesso del paziente -21.93291006 0b_6 f(V41) f(Durata di ospedalizzazione) 0.807408389 0b_7 V86 Evento MDC5 -332.9717588 0

Un ultimo check riguarda l'errore monetario previsionale:

Δ(0≤Y≤12000) (errore percentuale) = -2.297.287 € (-1.20%)

Essendo questo valore più basso di quello ottenuto nella sezione 5.2.2 possiamoconcludere che l'aggiunta delle variabili binarie della tabella precedentecontribuiscono a ridurre l'errore previsionale.

n134

5.2.5 Analisi di rischio

Avvedo lasciato in sospeso l'analisi delle code per la variabile ricoveri valorizzati (V29), riprendiamo questo tema, sia per completezza del modello, sia per l'idoneità dell'analisi che stiamo effettuando. Vedremo che il tema ha una importanza sia dal punto di vista matematico, sia dal punto di vista manageriale. I sistemi di controllo di gestione prevedono anche un interfaccia con la funzione del risk management, in particolare secondo Azzone (2007)[84] i sistemi di budgeting possono essere classificati in:

• Sistemi che escludono il rischio• Sistemi che considerano il rischio a valle della definizione del target (approccio

sequenziale)• Sistemi che definiscono in modo integrato target e livello di rischio

La scrittura matematica della funzione fY(y|Y B), nel framework che stiamo usando, ϵcorrisponde a un'analisi di rischio di tipo sequenziale. Ricordando che qui non stiamolavorando con la possibilità di intervenire sul fenomeno sotto studio, è interessante notare che la scrittura di fY(y) o Pr(Y<y) usando il teorema delle probabilità totali, è concettualmente analogo a integrare un sistema di controllo di gestione con un sistema di risk management. In base al metodo proposto prima l'insieme B=[12000,+∞) contiene gli eventi estremi, i quali oltre ad essere estremi sono anche rari e quindi hanno una probabilità non stimabile con precisione: la probabilità nelle code è piccola rispetto al centro e quindi l'errore esclusivamente probabilistico rimane piccolo in termini assoluti, ma alto in termini relativi, in più siccome gli eventi estremi sono rari si avrà tipicamente insufficienza di sampling e la probabilità degli eventi estremi viene necessariamente sottovalutata in base alle frequenze empiriche osservabili.Il termine rischio indica però una funzione di probabilità*magnitudo [84][43], ovvero della pesantezza della coda della distribuzione, segue quindi che per eventi estremi piccoli errori di stima nella loro probabilità possono risultare in errori giganteschi nella stima del rischio, questo può essere agevolmente verificato con un'analisi di sensitività dei parametri della distribuzione per la coda.Resta però il problema di quale forma funzionale scegliere per rappresentare la coda. Un procedimento conservatore è supporre una distribuzione di Pareto, intuitivamentesiccome non sappiamo quanto pesante sia la coda supponiamo di avere la coda più pesante (in realtà esistono code più pesanti di quelle paretiane però con α≤1 abbiamo già il caso della media infinita che è completamente distruttivo senza il bisogno di andare oltre). D'altra parte con α→∞ si possono ottenere code molto sottili, la distribuzione di Pareto si mostra quindi molto versatile per le analisi di sensitività di rischio. Abbiamo quindi la seguente rappresentazione:

n135

Con ym=12000 euro e α il parametro dell'analisi di sensitività.Quindi per il nostro problema abbiamo che:

n136

Quindi dal primo grafico segue che 3.93<α*+1<5.76, dove α*+1 ricavabile dal grafico, dal secondo grafico abbiamo che α*=5.17, avevamo già che α*>3, quindi sicuramente α<5.76 (qui α inteso come il parametro “vero” della distribuzione). È importante non escludere a priori che α, possa essere più basso di 5.76, una volta che è stato verificato le calibrazioni del parametro α* hanno una tendenza a ricavare un α*>α[116].È possibile che ci sia un limite massimo al rischio, però questo limite non può essere identificato con tecniche di tipo campionario, ma necessitano di un giudizio a priori sulle proprietà della variabile aleatoria[88].Alla luce di queste considerazioni avremo quindi che le funzioni di ripartizione condizionate sono:

n137

Vediamo adesso di ricostruire Pr(Y<y) e fY(y) usando il teorema delle probabilità totali e imponendo:

Indichiamo per brevità Pr(A)=p e Pr(B)=1-p. Dal teorema della probabilità totali segue che:

Per il nostro caso p≈0.8. La funzione di densità più semplice immaginabile per questafunzione di ripartizione è:

Ricaviamo anche la media per questa fdp:

La funzione di densità così costruita è in generale discontinua, che non èincompatibile con una variabile aleatoria continua o assolutamente continua, vistoche la caratteristica fondamentale di quest'ultima è la continuità della funzione diripartizione che non viene violata. Se si volesse comunque costruire una fdp continuasi impone:

Dove O=12000 euro nel nostro caso. Così facendo il parametro della Pareto diventalegato a un parametro della NIG.Il rischio (per Regione Lombardia) nel nostro caso però non è patito a livello diPr(Y<y) ma a livello di costo totale ∑Y. Illustriamo la problematica con una versione

n138

modificata del tradizionale grafico ABC, in questa formulazione nelle ascisseabbiamo valori assoluti e non percentuali mentre nelle ordinate abbiamo la solitafunzione di costo cumulato totale in termini percentuali18, questa rappresentazionerisulta più intuitiva per il nostro problema:

Entrambi i grafici fanno riferimento agli anni 2002-2007, con la differenza che ilprimo è per l'insieme A=[0,12000] mentre il secondo è per tutto il dominio di analisi.La linea rossa individua il quantile 80% e la linea gialla individua il quantile 95%.Qui possiamo avere una idea intuitiva del ruolo che giocano le code:

• La funzione di costo totale cumulata è convessa-concava, mentre nell'analisi

18 Tradizionalmente si usa l'analisi ABC per il fatturato, però qui vedremo che risulta molto illutrativa anche per le analisi di costo

n139

ABC tradizionale abbiamo funzioni solitamente concave• Possiamo vedere che diversamente da quanto avviene nella distribuzione

gaussiana o nella distribuzione gamma, gli eventi sulla coda offrono uncontribuito non trascurabile al totale, inoltre nel nostro caso questo contributo èdifficile da stimare a priori, vediamo come l'aggiunta dell'insieme B sposta lapercentuale di costo cumulato totale. Per più dettagli su questo punto si vedaTaleb e Doaudy (2014) [89]

Qui ipotizziamo un controllo di gestione che funzioni con cadenza annuale, il costoannuo sarà quindi dato per definizione da:

L'espressione analitica esatta per la distribuzione di probabilità del costo annuo non èagevole da ricavare, infatti anche supponendo di conoscere con certezza NTOT ilteorema del limite centrale non ci può aiutare per i motivi prima esposti nellapresentazione del tableau di Taleb. Se ammettiamo di usare una coda paretiana comeprima esposto, allora vale il teorema generalizzato del limite centrale, il qualegarantisce la convergenza debole della somma ΣY alla distribuzione stabile.[91]D'altra parte per quanto riguarda l'ottica del controllo di gestione si valutano sempreprevisioni puntuali e la differenza di queste dai valori che si sono effettivamenteverificati e non di valori che potrebbero aver avuto luogo ma non si sono verificati,quest'ultima funzione è tipicamente riservata al risk management. Alla luce di ciò,quello che proponiamo qui è un buon compromesso tra l'interpretazionefrequentistica (Richard von Mises) e soggettiva (Bruno de Finetti)19 del concetto diprobabilità, che consiste nel comparare il valore atteso del costo annuo con il suovalore effettivo misurato. Avremo quindi che il costo annuo sarà dato da:

La funzione di interesse qui è la seguente:

19 Nell'interpretazione soggettiva l'uso del baricentro E(X) come valore atteso è standard (da qui infatti proviene la denominazione “valore atteso”), per quanto riguarda invece l'interpretazione frequentistica abbiamo che NTOT che è un numero abbastanza alto se prendiamo un intervallo annuale e quindi NTOTE(Y) approssima con accuratezza il valore misurato di costo totale.

n140

Si vede rapidamente che costo annuo=F(b)-F(a). Per questa analisi saraà sufficiente considerare a=0 e b=max(Y). La F(ybound) può essere stimata nel nostro caso usando una b-spline (curva in rosso):

Riassumiamo quanto detto finora:

• Usiamo come ipotesi di modellazione per la variabile V29 (valore economicodegli eventi di ricovero) una distribuzione NIG+Pareto

• Abbiamo stimato una curva di regressione di V29 su V41 (durata diospedalizzazione)

• Abbiamo stimato una curva di costo totale cumulato che ci consente dianalizzare l'impatto di eventi estremi

• Questa rappresentazione matematica unifica sia alcuni aspetti dell'approcciocontrollo di gestione, sia dell'approccio risk management

Cosi come per il controllo di gestione si cerca di prendere decisioni in grado di modificare le variabili di input I, nel risk management propriamente detto si può agire[84]:

• evitando i rischi

n141

• riducendo il rischio • trasferendo il rischio (tipicamente con forme di tipo assicurativo)

Qui valuteremo semplicemente lo status quo, senza considerare altre alternative. Perscrivere matematicamente il rischio associato alla curva di costo totale cumulato,consideriamo quest'ultima come una funzione di perdita, cioè -∑Y:

Evidentemente c'è anche una funzione di incasso per la Regione, la quale però non èl'oggetto di studio di questo lavoro. L'approccio che presenteremo qui è abbastanzagenerico e quindi indicheremo per semplicità la funzione di perdita come Z=φ(Y)senza fare riferimento esplicito alle variabili che abbiamo considerato finora. Inquesta ottica possiamo definire il rischio shortfall ζ per la variabile Z come:

E il rischio per la variabile Y come:

n142

Dove e sono la semi

deviazione assoluta per la coda sinistra (per la variabile Y e Z rispettivamente) e Ωuna costante di centralità della distribuzione. Ora il punto principale: quanto più spesse sono le code tanto più l'errore di stima di ζcresce, e per insufficienza di sampling nelle code non possiamo stimare stimaredirettatemente ζ. La situazione può sembrare paradossale in un primo momento, perònon è senza soluzione. Esiste un paper di Taleb e Douady (2013)[90] il quale stabilisce che, definendo

come la fragilità20 di Z “ereditata” da Y, il

teorema di aggravamento della fragilità garantisce che:

1. Se φ è concava nell'intervallo (-∞,k'] con K<k'<Ω e lineare in [k',Ω]2. Allora Z è più fragile di Y in L=φ(K) che Y in K

E viceversa:

1.Se φ è convessa nell'intervallo (-∞,k'] con K<k'<Ω e lineare in [k',Ω]2.Allora Z è più robusta di Y in L=φ(K) che Y in K

Riassumiamo questo risultato in un grafico:

20 A volte gli autori usano anche il termine sensitività semi-vega

n143

Nel nostro caso la parte convessa dominasulla coda sinistra, il che implica unariduzione del rischio di F(ybound) rispetto aquello della variabile V29, è da notareperò che il punto di inversione dellaconcavità della funzione F(ybound)corrisponde al parametro di centralitàμ, che deve essere tenuto opportunamentesotto controllo per evitare un aumento deicosti totali in generale.

Evidentemente esiste un livello di costo totale insopportabile per la Regione, però inuna situazione robusta questo livello è più difficile da raggiungere rispetto asituazione fragile.

5.3 Analisi delle sub-popolazioni21

Avendo quindi completato questo esercizio di modellazione, passiamo alladescrizione delle sub-popolazioni costruite usando le variabili binarie individuatecome significative dall'analisi di regressione nella sezione 5.2.4, cioè studiando solo isottoinsiemi del data-set dove queste variabili assumono il valore di flag=1. La primasub-popolazione che analizziamo è quella identificata dagli eventi con presenza diriabilitazione, il termine Δ(1) indica la differenza ottenuta nella variabile dentro lasub-popolazione rispetto all'analisi esplorativa condotta nella sezione 5.1.2:

21 Nel marketing si usa il termine analisi di segmentazione, però qui il termine risulta poco preciso perchè abbiamo a che fare con sub-popolazioni non necessariamente mutuamente esclusive

n144

Notiamo che gli eventi che iniziano con riabilitazione finiscono quasi tutti conriabilitazione e sono quasi tutti eventi dove è stata realizzato un intervento di bypassaortocoronarico.

n145

n146

Per la sub-popolazione identificata dagli eventi che finiscono con una riabilitazione,possiamo osservare che come nella sub-popolazione precedente c'è una grandepresenza di pazienti che hanno avuto un ricovero pediatrico con meno di 2 anni dietà. Ricordiamo che il codice di reparto 56 corrisponde al reparto recupero eriabilitazione funzionale, atto a preparare il paziente per la dimissione dall'ospedale.La prossima sub-popolazione è quella identificata dalla terapia intensiva:

n147

n148

Qui notiamo che i codici di reparto 50, 26 e 08 corrispondono rispettivamenteall'unità coronarica, alla medicina generale e alla cardiologia.

n149

Dall'analisi di quest'ultima sub-popolazione, possiamo osservare che il DRGchirurgico è collegato principalmente all'angioplastica coronarica e al bypass. Per la prossima sub-popolazione il termine Δ(M-F) indica la differenza a favore dellasub-popolazione maschile in relazione a quella femminile:

n150

Ovvero i maschi hanno meno DRG chirurgici e ricoveri pediatrici rispetto allefemmine. Per la sub-popolazione individuata dagli eventi MDC5 non è statoosservato nessuna differenza particolarmente significativa per le variabili binarie ecategoriche.L'utilità dell'analisi condotta in questa sezione è mostrare come i problemi nella curadello scompenso cardiaco in Regione Lombardia non siano scollegati tra di loro, peresempio se si volesse ridurre la percentuale di pazienti soggetti a riabilitazione èmolto probabile che sia anche necessario cercare di diminuire la percentuale dipazienti soggetti a intervento con bypass aortocoronarico. Se si volesse ridurre lapercentuale di DRG chirurgici da un lato sembra opportuno ridurre la percentuale dimaschi ricoverati e dall'altro la presenza di DRG chirurgici sembra essere collegatacon una maggior presenza di interventi di angioplastica e bypass aortocoronarico.Questa sezione fornisce evidenza che i casi di recupero con le terapie e interventicitati, oltre ad essere più statisticamente costosi, tendono a coinvolgere più di unatipologia di intervento e terapia, il che suggerisce da un lato che per ridurre il costo dicura sarà necessario adottare una visione di tipo sistemico agendo su più variabiliallo stesso tempo, dall'altro lato la presenza di ricoveri pediatrici segnala l'importanzadi problemi non esclusivamente legati alla età del paziente e la necessità di usaremetodi di medicina preventiva, incluso il home care del paziente, ove possibile.

n151

Conclusione e sviluppi futuri

Lo scompenso cardiaco è una sindrome complessa che porta a una serie didebilitazioni e ed un fattore di rischio particolarmente importante: riuscire amigliorare le modalità di cura (efficcacia) senza un consumo eccessivo di risorse(effficienza) è quindi fondamentale per i policy maker sanitari.Siccome la popolazione arruolata nei trial clinici non è sufficiente per capire ladinamica di una malattia dentro una popolazione di una regione geografica, si fanecessario l'utilizzo di database amministrativi per studiare problemi di naturaepidemiologia e per fare analisi di tipo economico a livello di sistema sanitario.Il progetto HF DATA, finanziato dal Ministero della Salute e Regione Lombardia èstato creato proprio per supportare i policy maker regionale e nazionale nello studiodi questo problema. Sulla base dei dati disponibili è stato sviluppato un modello di regressione che sia ingrado di cogliere i fenomeni principali della cura dello scompenso cardiaco. Tra ipregi del modello sviluppato possiamo citare l'accuratezza della calibrazioneeffettuata, con l'obbiettivo principale di evitare che gli errori previsionaliaumentassero in modo incontrollabile, mostrando anche come ridurre il biasprevisionale attraverso l'utilizzo di un insieme di training mobile per il modello.Abbiamo anche mostrato come condurre un'analisi delle code compatibile con lostato attuale dell'arte e abbiamo mostrato l'importanza di questo tema per problemi dinatura economica e gestionale. Il modello sviluppato può quindi servire di supporto alcontrollo di gestione regionale, purché le previsioni del modello rimangano entrol'orizzonte annuale.Tra i possibili sviluppi futuri, a partire da questo tema, possiamo individuare:

• Uno studio più approfondito delle funzioni di probabilità cumulate a gradino,le quali renderebbero possibile l'ottenimento di un fit di distribuzione diprobabilità per i costi medici praticamente perfetto. Questo è un problemainteressante dal punto di vista della teoria della probabilità, perché ci obbliga aconsiderare variabili aleatorie di tipo misto, quindi né discrete, né continue.

• La distribuzione Normale Inversa Gaussiana (NIG), e più in generale ladistribuzione iperbolica generalizzata, sono distribuzioni che si sono presentatepiù volte nella letteratura statistica come fit di dati reali: fenomeni legati aiprocessi eolici, distribuzioni di rocce e cristalli nel sottosuolo, serie finanziariee adesso anche dati di costo in ambito medico. Se ci sia qualcosa di comune atutti questi fenomeni è difficile dire, però sarebbe interessante ottenere unaformula esplicita per la distribuzione della somma di variabili condistribuzione generale iperbolica, particolarmente in ambito economico, doveè necessario sommare una serie di realizzazioni di una variabile di costo.

n152

• L'analisi qui effettuata può essere estesa nella costruzione di un modello joint,che seppure più complicato, permetterebbe di tenere conto del fenomeno dellamortalità intraospedaliera. Per un modello joint è necessario stimare lafunzione di Kaplan-Meier, il che è particolarmente oneroso22, probabilmentesarà necessario fornire all'analista la funzione già calcolata insieme al data-set.Non è da scontare l'utilizzo delle potenzialità del cloud.

• Per effettuare un'analisi aggregando i dati nella logica del paziente èopportuno, oltre alla stima della funzione di Kaplan-Meier, correggere i valoriottenuti con la regressione in base alla durata temporale dell'insieme di traininge di validazione (vedere appendice).

22 È stato usato un notebook Intel core i7 per fare l'analisi

n153

Appendice

Passaggio matematico dalla logica della prevalenza alla logica dell'incidenza

Per passare dalla logica degli eventi alla logica del paziente ci sono 2 strade possibili.Sia Y la variabile di costo per evento, N il numero di eventi per paziente e Z lavariabile di costo per paziente, la prima possibilità è costruire le funzioni di costo perogni singolo paziente e usare un modello di tipo joint:

Oppure aggregare i dati di costo per paziente e usare il valore totale di costo dentro ilmodello:

É noto dalla letteratura la seguente proprietà della distribuzione NIG[81]: seY1~NIG(α,β,δ1,μ1) e Y2~NIG(α,β,δ2,μ2) e Y1 Y2 sono stocastica-mente indipendenti,allora Y1+Y2~NIG(α,β,δ1+δ2,μ1+μ2), segue quindi che per N Yi idd abbiamonecessariamente che Σyi~NIG(α,β,Nδ,Nμ). Questa proprietà semplifica notevolmentel'analisi, infatti la funzione di ripartizione è data da :

E il valore atteso può essere ricavato attraverso:

Quindi, per esempio se disponiamo di 8 anni di dati e vogliamo ottenere unaprevisione per i costi dei pazienti per 3 anni usando 5 anni come insieme di training,basta dividere l'output della regressione per il numero medio di eventi in 5 anni e poimoltiplicare il quoziente per il numero medio di eventi in 3 anni.

n154

Bibliografia

[1]Michael Y. Henein, Heart Failure in Clinical Practice, Springer, 2010[2]Piazza E., Probabilità e Statistica, Progetto Leonardo 2007[3]Vercellis C., Business Intelligence, McGrawHill 2006[4]Dobson, A.J. , Introduction to Generalized Linear Models, Chapman andHall/CRC 2002[5]Sito Ministero della salute http://www.salute.gov.it/portale/temi/p2_6.jsp?lingua=italiano&id=1232&area=ricoveriOspedalieri&menu=vuoto[6]Diagnosis Related Groups and the Medicare Program: Implications for MedicalTechnology, Washington D.C. U.S. Congress, Office of Technology Assessment 1983[7]Fetter Robert B., Shin Y, Freeman JL, Averill RF, Thompson JD , Case MixDefinition by Diagnosis-Related Groups, Medical Care 1980[8]Mills RE, Fetter RB, Riedel DC, Averill RF, AUTOGRP: an interactive computersystem for the analysis of health care data, Medicare 1976[9]Jing Fang, George A. Mensah, Janet B. Croft, Nora L. Keenan, Heart Failure-Related Hospitalization in the U.S., 1979 to 2004, Journal of the American College ofCardiology, 2008[10]Simon Stewart, Andrew Jenkins, Scot Buchan, Alistair McGuire, SimonCapewell, John J.J.V. McMurray, The current cost of heart failure to the NationalHealth Service in the UK, The european journal of heart failure 2002[11]Dunlay S.M., Shah N.D., Shi Q., Morlan B., VanHouten H., Long K.H., RogerV.L., Lifetime costs of medical care after heart failure diagnosis, Circulation:Cardiovascular Quality and Outcomes, 2011[12]Lu Tian, Jie Huang, A two-part model for censored medical cost data, Statisticsin Medicine 2007[13]James M. Robins, Andrea Rotnitzky and Lue Ping Zhao, Estimation ofRegression Coefficients When Some Regressors Are Not Always Observed, Journal ofthe American Statistical Association 1994[14]Manning WG, Mullahy J. Estimating log models: to transform or not to transform? Journal of Health Economics 2001[15]Esposito D., Bagchi A.D., Verdier J.M., Bencio D.S., Kim M.S., Medicaidbeneficiaries with congestive heart failure: Association of medication adherence withhealthcare use and costs, American Journal of Managed Care 2009[16]Melinda Beeuwkes Buntin, Alan M. Zaslavsky, Too much ado about two-partmodels and transformation? Comparing methods of modeling Medicare expenditures,Journal of Health Economics 2004[17]Titler M.G., Jensen G.A., Dochterman J.M., Xie X.-J., Kanak M., Reed D.,Shever L.L.,Cost of hospital care for older adults with heart failure: Medical,pharmaceutical, and nursing costs, Health Services Research 2008[18]Kung-Yee Liang, Scott L. Zeger, Longitudinal Data Analysis Using Generalized

n155

Linear Models, Biometrika 1986[19]Setoguchi S., Glynn R.J., Stedman M., Flavell C.M., Levin R., Stevenson L.W.,Hospice, opiates, and acute care service use among the elderly before death fromheart failure or cancer, American Heart Journal 2010[20]Guangyong Zou, A Modified Poisson Regression Approach to ProspectiveStudies with Binary Data, American Journal of Epidemiology 2004[21]Jayadevappa R., Chhatre S., Weiner M., Raziano D.B., Health resourceutilization and medical care cost of acute care elderly unit patients, Value in Health2006[22]Student (pseudonimo di W. S. Gosset), The probable error of a mean Biometrika1908 [23]Somnath Datta, Bootstrap, Encyclopedia of Statistical Sciences Wiley 2006[24]Frank Wilcoxon, Individual Comparison by Ranking Methods, BiometricsBulletin 1945[25]Naihua Duan, Smearing Estimate: A Nonparametric Retransformation Method,Journal of the American Statistical Association 1983[26]Kaul P., McAlister F.A., Ezekowitz J.A., Bakal J.A., Curtis L.H., Quan H.,Knudtson M.L., Armstrong P.W., Resource use in the last 6 months of life amongpatients with heart failure in Canada, Archives of Internal Medicine 2011[27]William G. Cochran, Some Methods for Strengthening the Common χ2 Tests,Biometrics 1954[28]Mantel N., Haenszel W., Statistical aspects of the analysis of data fromretrospective studies of disease, J Natl Cancer Inst. 1959[29]Mantel N., Chi-Square Tests with One Degree of Freedom, Extensions of theMantel- Haenszel Procedure, Journal of the American Statistical Association 1963[30]Howell S., Coory M., Martin J., Duckett S., Using routine inpatient data toidentify patients at risk of hospital readmission, BMC Health Services Research 2009[31]Hosmer D, Lemeshow S, Applied Logistic Regression Second edition. New York, John Wiley & Sons 2000[32]Whitfield M.D., Gillett M., Holmes M., Ogden E., Predicting the impact ofpopulation level risk reduction in cardio-vascular disease and stroke on acutehospital admission rates over a 5 year period-a pilot study, Public Health 2006[33]Anderson K.M., Odell P.M., Wilson P.W., Kannel W.B., Cardiovascular diseaserisk profiles, Am Heart J 1991[34]Richard J. Stevens, Viti Kothari, Amanda I. Adler, Irene M. Stratton and Rury R.Holman, The UKPDS risk engine: a model for the risk of coronary heart disease inType II diabetes (UKPDS 56), Clinical Science 2001[35]Graham Kalton, Systematic Sampling, Encyclopedia of Statistical Sciences Wiley2006[36]Hung-Ir Li, Longitudinal Data, Encyclopedia of Statistical Sciences Wiley 2006[37]Raymond S. Greenberg, Retrospective Studies, Encyclopedia of Statistical

n156

Sciences Wiley 2006[38]Wijeysundera H.C., Wang X., Tomlinson G., Ko D.T., Krahn M.D., Techniquesfor estimating health care costs with censored data: An overview for the healthservices researcher, ClinicoEconomics and Outcomes Research 2012[39]E. L. Kaplan and Paul Meier, Nonparametric Estimation from IncompleteObservations, Journal of the American Statistical Association 1958[40]Lin DY, Feuer EJ, Etzioni R, Wax Y., Estimating medical costs from incompletefollow-up data, Biometrics 1997[41]Bang H, Tsiatis AA., Estimating medical costs with censored data, Biometrika2000[42]Raikou M, McGuire A., Estimating medical care costs under conditions ofcensoring., J Health Econ. 2004[43]Antonio Calabrese, Gestione degli impianti industriali - Volume I, CUSL 2004[44]Roderick J.A. Little, Donald B.Rubin, Statistical Analysis with Missing Data,Wiley-Interscience 2002[45]Rubin D., Multiple Imputation for Nonresponse in Surveys, John Wiley and Sons1987[46]Gregory D., DeNofrio D., Konstam M.A., The economic effect of a tertiaryhospital-based heart failure program, Journal of the American College of Cardiology2005[47]Meeker WQ, Escobar LA, Statistical Methods for Reliability Data, John Wileyand Sons Inc. 1998[48]Friedman B, Pauly MV, Cost functions for a service firm with variable quality and stochastic demand, Rev Econ Stat 1981[49]Sun L., Song X., Zhou J., Liu L., Joint analysis of longitudinal data withinformative observation times and a dependent terminal event, Journal of theAmerican Statistical Association 2012[50]http://hfdata.cefriel.it/[51]Testo aggiornato del decreto legislativo 30 dicembre 1992, n. 502 recante:“Riordino della disciplina in materia sanitaria, a norma dell’art.1 della legge 23ottobre 1992, n. 421” [52]Regione Lombardia, Progetto CRS-SISS: secondo aggiornamento del FlussoUnico di Rendicontazione delle Farmacie (FUR 03), protocollo H1.2008.0029278[53]https://support.3mhis.com/app/answers/detail/a_id/9907/~/definitions-manuals[54]R Core Team (2014). R: A language and environment for statistical computing. RFoundation for Statistical Computing, Vienna, Austria. URL http://www.R-project.org/[55]ALL PATIENT REFINED DIAGNOSIS RELATED GROUPS (APR-DRGs),Version 20.0, Methodology Overview, 3M Health Information Systemshttp://www.hcup-us.ahrq.gov/db/nation/nis/APR-DRGsV20MethodologyOverviewandBibliography.pdf

n157

[56]History of the development of the ICD,http://www.who.int/classifications/icd/en/HistoryOfICD.pdf[57]Heart Failure: Codice 428 nella codifica ICD9-CMhttp://www.icd9data.com/2014/Volume1/390-459/420-429/428/default.htm[58]Deliberazione della Giunta Regionale 11637[59]Deliberazione della Giunta Regionale 10804, vedi allegato 12[60]http://www.sanita.regione.lombardia.it/cs/Satellite?c=Page&childpagename=DG_Sanita%2FDGLayout&cid=1213480043116&p=1213480043116&pagename=DG_SANWrapper Ricerca testo nell'oggetto=Grouper[61]Richard F. Averill, M.S., John H. Muldoon, M.H.A., James C. Vertrees, Ph.D.,Norbert I. Goldfield, M.D, Robert L. Mullin, M.D., Elizabeth C. Fineran, M.S., MonaZ. Zhang, B.S., Barbara Steinbeck, ART, Thelma Grant, RRA, The Evolution ofCasemix Measurement Using Diagnosis Related Groups (DRGs), Stud HealthTechnol Inform. 1994[62]Guida tecnica, Descrizione delle principali modificazioni tra la versione 19.0 ela versione 24.0 del sistema di classificazione Medicare DRG, 3M SistemiInformativi per la Salute 2008[63]Marino Nonis , Enrico Rosati, Guida alla classificazione degli interventichirurgici, ISTITUTO POLIGRAFICO E ZECCA DELLO STATO 2013[64]http://www.cdc.gov/nchs/icd/icd9cm.htm[65]Dall T.M., Askarinam Wagner R.C., Zhang Y., Yang W., Arday D.R., Gantt C.J.,Outcomes and lessons learned from evaluating TRICARE's disease managementprograms, American Journal of Managed Care 2010[66]Adrian G Barnett, Jolieke C van der Pols, Annette J Dobson, Regression to themean: what it is and how to deal with it, International Journal of Epidemiology 2005[67]Garrett N., Martini E.M., The boomers are coming: A total cost of care model ofthe impact of population aging on the cost of chronic conditions in the United States,Disease Management 2007[68]U.S. Census Bureau. Interim Projections of the U.S. population by age, sex, raceand Hispanic Origin. www.census.gov[69]Taleb, Nassim N., Silent Risk: Lectures on Fat Tails, (Anti)Fragility, Precaution,and Asymmetric Exposures (August 5, 2014). Available at SSRN:http://dx.doi.org/10.2139/ssrn.2392310[70]Li Z., Amstrong E.J., Parker J.P., Danielsen B., Romano P.S., Hospital variationin readmission after coronary artery bypass surgery in California, Circulation:Cardiovascular Quality and Outcomes 2012[71]Liu L., Huang X., O'Quigley J., Analysis of longitudinal data in the presence ofinformative observational times and a dependent terminal event, with application tomedical cost data, Biometrics 2008[72]Liu L., Joint modeling longitudinal semi-continuous data and survival, with

n158

application to longitudinal medical cost data, Statistics in Medicine 2009[73]Paul Embrechts, Modelling extremal events: for insurance and finance, Springer1997[74]Mandelbrot, Benoit B., Fractals and Scaling in Finance: Discontinuity,Concentration, Risk, Springer, 1997[75]V. P. Chistyakov, A Theorem on Sums of Independent Positive Random Variablesand Its Applications to Branching Random Processes, Theory of Probability and ItsApplications 1964[76]J.L. Teugels, The Class of Subexponential Distributions, Annals of Probability1975[77]Morris H. Degroot, Record Linkage and Matching Systems, Encyclopedia ofStatistical Sciences Wiley 2006[78]Goldstein H, Browne W, Rasbash J., Partitioning variation in multilevel models, Understand Statistics 2002[79]Smirnov N., Table for estimating the goodness of fit of empirical distributions,Annals of Mathematical Statistics 1948[80]Barndorff-Nielsen O., Exponentially Decreasing Distributions for the Logarithmof Particle Size, Proceedings of the Royal Society of London 1977[81]Ole E Barndorff-Nielsen, Thomas Mikosch and Sidney I. Resnick, LévyProcesses: Theory and Applications, Birkhäuser 2013[82]Taleb, Nassim Nicholas and Tetlock, Philip E., On the Difference between BinaryPrediction and True Exposure with Implications for Forecasting Tournaments andDecision Making Research (November 27, 2013). Available at SSRN:http://ssrn.com/abstract=2284964[83]Azzone G., Bertelè U., L'impresa. Sistemi di governo, valutazione e controllo,ETAS 2007[84]Giovanni Azzone, Sistemi di controllo di gestione, ETAS 2006[85]Wolfgang Breymann, David Lu thi, ghyp: A package on generalized hyperbolic distributions, 2013 http://cran.at.r-project.org/web/packages/ghyp/index.html[86]Abramowitz M., Stegun I.A., Handbook of Mathematical Functions withFormulas, Graphs, and Mathematical Tables, National Bureau of Standards AppliedMathematics Series 1972[87]Hirotugu Akaikae, A new look at the statistical model identification, IEEETransactions on Automatic Control 1974[88]Douady, Raphael and Taleb, Nassim Nicholas, Statistical Undecidability(October 12, 2010). Available at SSRN: http://ssrn.com/abstract=1691165[89]Taleb, Nassim Nicholas and Douady, Raphael, On the Super-Additivity andEstimation Biases of Quantile Contributions (November 11, 2014). Available atSSRN: http://ssrn.com/abstract=2434363[90]N.N.Taleb, R.Douady, Mathematical definition, mapping, and detection of(anti)fragility, Quantitative Finance 2013

n159

[91]Vladimir V. Uchaikin, Vladimir M. Zolotarev, Chance and Stability: StableDistributions and their Applications, De Gruyter 1999[92]S De Lusignan, An educational intervention, involving feedback of routinelycollected computer data, to improve cardiovascular disease management in UKprimary care, Methods Inf Med 2007[93]C van Walraven, P Austin, Administrative database research has uniquecharacteristics that can risk biased results, Journal of clinical epidemiology 2012[94]Kuwabara K1, Matsuda S, Anan M, Fushimi K, Ishikawa KB, Horiguchi H,Hayashida K, Fujimori K., Difference in resource utilization between patients withacute and chronic heart failure from Japanese administrative database, Int J Cardiol2010[95]Vijaya Sundararajana, Toni Henderson, Catherine Perry, Amanda Muggivan,Hude Quan, William A. Ghali New ICD-10 version of the Charlson ComorbidityIndex predicted in-hospital mortality, J Clin Epidemiol 2004[96]Yu-Chun Chen, Hsiao-Yun Yeh, Jau-Ching Wu, Ingo Haschler, Tzeng-Ji Chen,Thomas Wetter, Taiwan's National Health Insurance Research Database:administrative health care database as study object in bibliometrics, Scientometrics2011[97]Saczynski JS1, Andrade SE, Harrold LR, Tjia J, Cutrona SL, Dodd KS, GoldbergRJ, Gurwitz JH., A systematic review of validated methods for identifying heartfailure using administrative data, Pharmacoepidemiol Drug Saf. 2012[98]Quach S1, Blais C, Quan H., Administrative data have high variation in validityfor recording heart failure, Can J Cardiol 2010[99]Ronald Fisher, Statistical Methods for Research Workers 1925[100]Efron B., R. J. Tibshirani, An introduction to the bootstrap, Chapman &Hall/CRC 1998[101]Sebastian Schneeweiss, Jerry Avorn, A review of uses of health care utilizationdatabases for epidemiologic research on therapeutics, Journal of ClinicalEpidemiology 2005[102]Joshua J. Gagne, Robert J. Glynn, Jerry Avorn, Raisa Levin, SebastianSchneeweiss, A combined comorbidity score predicted mortality in elderly patientsbetter than existing scores, Journal of Clinical Epidemiology 2011[103]Riunione progetto HFData-AO Niguarda Ca’ Granda-27/03/2014[104]Alan S. Go, Dariush Mozaffarian, Véronique L. Roger, Emelia J. Benjamin,Jarett D. Berry, William B. Borden, Dawn M. Bravata, Shifan Dai, Earl S. Ford,Caroline S. Fox, Sheila Franco, Heather J. Fullerton, Cathleen Gillespie, Susan M.Hailpern, John A. Heit, Virginia J. Howard, Mark D. Huffman, Brett M. Kissela,Steven J. Kittner, Daniel T. Lackland, Judith H. Lichtman, Lynda D. Lisabeth, DavidMagid, Gregory M. Marcus, Ariane Marelli, David B. Matchar, Darren K. McGuire,Emile R. Mohler, Claudia S. Moy, Michael E. Mussolino, Graham Nichol, Nina P.Paynter, Pamela J. Schreiner, Paul D. Sorlie, Joel Stein, Tanya N. Turan, Salim S.

n160

Virani, Nathan D. Wong, Daniel Woo and Melanie B. Turner, Heart Disease andStroke Statistics −− 2013 Update: A Report From the American Heart Association,Circulation 2012[105]Frieder Braunschweig, Martin R. Cowie, Angelo Auricchio, What are the costsof heart failure?, Eurospace 2011[106]Douglas Gregory, David DeNofrio, Marvin A. Konstam, Douglas Gregory,David DeNofrio, Marvin A. Konstam, The Economic Effect of a Tertiary Hospital-Based Heart Failure Program, Journal of the American College of Cardiology 2005[107]Paul J. Hauptman, Jason Swindle, Thomas E. Burroughs, Mark A. Schnitzler,Resource utilization in patients hospitalized with heart failure: Insights from acontemporary national hospital database, Am Heart J. 2008[108]Karen Hogg, Karl Swedberg, John McMurray, Heart Failure With PreservedLeft Ventricular Systolic Function Epidemiology, Clinical Characteristics, andPrognosis, Journal of the American College of Cardiology 2004[109]Shelby D. Reed, Padma Kaul, Yanhong Li, Zubin J. Eapen, Linda Davidson-Ray, Kevin A. Schulman, Barry M. Massie, Paul W. Armstrong, Randall C. Starling,Christopher M. O’Connor, Adrian F. Hernandez, Robert M. Califf, Medical ResourceUse, Costs, and Quality of Life in Patients With Acute Decompensated Heart Failure:Findings From SCEND-HF, Journal of Cardiac Failure 2013[110]Barbora Rihova, Jiri Parenica, Jiri Jarkovsky, Roman Miklik, AlexandraSulcova, Simona Littnerova, Marin Felsoci, Petr Kala, Jindrich Spinar, Cardiologydepartment hospitalization costs in patients with acute heart failure vary accordingto the etiology of the acute heart failure: Data from the AHEAD Core registry 2005–2009, Cor et Vasa 2013[111]Rodriguez F, Wang Y, Johnson CE, Foody JM, National patterns of heart failurehospitalizations and mortality by sex and age, J Card Failure 2013[112]Rossano JW, Kim JJ, Decker JA, Price JF, Zafar F, Graves DE, Morales DL,Heinle JS, Bozkurt B, Towbin JA, Denfield SW, Dreyer WJ, Jefferies JL, Prevalence,morbidity, and mortality of heart failure-related hospitalizations in children in theUnited States: a population-based study, J Card Fail. 2012[113]Wijeysundera HC, Wang X, Tomlinson G, Ko DT, Krahn MD., Techniques forestimating health care costs with censored data: an overview for the health servicesresearcher, Clinicoecon Outcomes Res. 2012[114]http://www.heart.org/HEARTORG/General/History-of-the-American-Heart-Association_UCM_308120_Article.jsp[115]http://data.worldbank.org/indicator/NY.GDP.MKTP.CD[116]Manuela Casula, Elena Tragni, Alberico L. Catapano, Utilizzo dei databaseamministrativi nella ricerca sanitaria, CARE 2, 2011[117]Pollastri, A., Analysis of the kurtosis in a skew distribution, Metron 1999