Parte III, a cura di: Stefano Parodi, Servizio di Epidemiologia e Biostatistica, Direzione Scientifica, Istituto G. Gaslini, Largo G. Gaslini, 5 – 16147 Genova Tel: 010 5636301, e-mail: stefanoparodi@ospedale-gaslini.ge.it

Esercitazione pratica nell’ambito del Corso di Modelli Statistici

Anno Accademico 2005-2006

Applicazione del modello di regressione logistica nell’ambito di uno studio caso-

controllo

Corso di Laurea in

Statistica Matematica e trattamento Informatico dei Datiati

Università di Genova

Metodi per la selezione di un modello di regressione multivariata

Avendo a disposizione numerosi predittori si possono ottenere moltissimi modelli.

Tuttavia il ricercatore desidera pervenire ad un solo oppure a pochissimi modelli in grado di descrivere in modo sintetico ed esaustivo l’effetto delle variabili in studio.

La selezione del modello può essere effettuata attraverso diversi metodi, ovvero: backward, forward, stepwise, “best model”.

Il metodo backward, che verrà utilizzato nell’esercitazione, inserisce in blocco tutte le variabili nel modello e quindi le elimina una alla volta in base alla significatività dei rispettivi coefficienti.

In genere il test utilizzato è il LR test.

Tale procedura viene ripetuta ricorsivamente finché solamente le variabili con coefficiente statisticamente significativo rimangono all’interno del modello.

Tuttavia molti Autori suggeriscono di lasciare all’interno del modello anche predittori non significativi, che però siano noti confondenti.

Il metodo forward consiste nell’inserire le variabili nel modello una alla volta sula base della loro significatività.

Tale metodo è più conservativo, e particolarmente adatto quando il numero dei predittori è relativamente elevato rispetto al numero dellle osservazioni.

Il metodo bakward tende invece a sovrastimare il fitting del modello, ma risulta più adatto per l’analisi di data set contenenti numerose osservazioni.

Presenta inoltre il vantaggio di poter controllare per il confondimento simultaneo da due o più variabili (joint confounding).

Il metodo stepwise, consiste nell’inserire e rimuovere ricorsivamente le variabili dal modello sulla base della loro significatività statistica.

Può essere quindi applicato sia in al metodo backward che a quello forward, costituendo di fatto una variante di tali metodi.

Il metodo “best model” consiste nel fittare tutti i modelli possibili e nello scegliere quello con la migliore bontà di adattamento.

È caduto ormai quasi completamente in disuso, poiché tende a sovrastimare largamente il fitting del modello, generando associazioni spurie.

Modelli con variabili nested

Talvolta il ricercatore deve analizzare variabili che sono presenti solamente all’interno di alcuni livelli di altri predittori.

Ad esempio, tali variabili possono rappresentare una misura dell’intensità dell’esposizione, che esiste ovviamente solamente nel sottogruppo degli esposti.

In tale caso è possibile evitare il modellamento di tipo nested attribuendo ai non esposti un’intensità di esposizione pari a 0.

Tale approccio non è sempre consigliabile, in quanto spesso gli esposti differiscono dai non esposti per una serie di variabili il cui effetto non è controllabile dal ricercatore (in quanto molte di esse non sono esattamente note).

Ad esempio, i fumatori differiscono dai non fumatori per una serie di abitudini di vita e dietetiche che non sono tutte facilmente registrabili in fase di raccolta dati.

In sintesi, si suppone che l’effetto del fumo possa racchiudere anche effetti non misurati di altre variabili e, di conseguenza, che non si possa assumere come zero il valore (teorico) del rischio di patologia in assenza dell’esposizione.

Inoltre esistono altre variabili, in genere anch’esse legate a caratteristiche dell’esposizione, che non possono essere introdotte nel modello se non come variabili nested (cioè annidate in altre variabili)

Un esempio è dato dal tempo di cessazione dall’abitudine al fumo, che si suppone debba essere associato in modo inverso al rischio di malattia.

Infatti tale variabile è “annidata” entro la variabile “abitudine al fumo”, essendo misurata esclusivamente per i soggetti ex- fumatori.

Per i fumatori correnti sarebbe teoricamente possibile calcolare tale variabile, pur con i limiti sopra accennati, attribuendole valore zero.

Tuttavia per i non fumatori tale variabile non è calcolabile.

In tal caso un modello “nested” permette di valutarne l’effetto senza dover ricorrere ad una restrizione (cioè senza togliere dal modello i non fumatori).

Una procedura adeguata per fittare modelli con variabili nested su scala continua, consiste, come primo passaggio, nel centering della variabile stessa.

Ovvero alla variabile viene sottratto il suo valor medio entro gruppo.

Si supponga che la variabile principale di esposizione sia a tre livelli, ad esempio: fumatori, non fumatori ed ex-fumatori.

Il secondo passaggio dell’analisi consiste nel sostituire il valore zero ai soggetti con valore mancante (nell’esempio i non fumatori).

Si procede quindi a generare due nuove variabili realizzando il prodotto tra le due dummy variables che descrivono l’esposizione e i valori della variabile nested centrata.

Si supponga inoltre che la variabile sia annidata entro due livelli di esposizione (es: intensità dell’esposizione a fumo di tabacco entro i due gruppi di fumatori, ex- e attuali)

Siano D1 e D2 le due dummy per l’esposizione, con D1 che

assume valore 1 se il soggetto è un ex-fumatore e 0 altrimenti e la dummy D2 che vale 1 se il soggetto è un fumatore corrente e

vale 0 altrimenti. .

Siano D1Fc e D2Fc le due variabili per l’effetto nested

dell’esposizione a fumo (dove il pedice c ricorda il fatto che la variabile è stata centrata, per cui presenta media pari a zero).

Un modello di regressione logistica può quindi essere formulato nella maniera seguente:

CFDFDDDYit cc 5,224,11322110ˆˆˆˆˆˆˆlog

dove 5C indica l’effetto di uno o più confondenti.

CFDFDDDYit cc 5,224,11322110ˆˆˆˆˆˆˆlog

Poiché la variabile F è centrata sulla sua media il coefficiente stimato 1 opportunamente esponenziato, rappresenta una stima

di rischio relativo (OR) tra gli ex-fumatori che hanno consumato un numero medio di sigarette e i non-fumatori.

Analogamente esponenziando 2 si ottiene l’OR tra i fumatori che

hanno fumato una quantità media di sigarette e i non-fumatori.

I due coefficienti per la variabile nested (3 e 4) rappresentano

invece una stima dell’effetto dell’esposizione entro i due sottogruppi di esposti (ovvero una misura della forza della relazione dose-risposta).

Essendo Fc una variabile continua, i due coefficienti 3 e 4

rappresentano una stima della variazione media nel valore

dell’OR in ognuno dei due gruppi di esposti al variare di un’unità della variabile di esposizione (nell’esempio, migliaia di sigarette fumate per anno).

CFDFDDDYit cc 5,224,11322110ˆˆˆˆˆˆˆlog

Valutare la bontà di adattamento di un modello di regressione logistica

Una volta selezionato un modello che descriva i dati sulla base delle ipotesi formulate nel disegno dello studio occorre procedere ad una valutazione della bontà del modello stesso.

È quindi necessario ricavare misure quantitative della capacità del modello di interpolare i dati osservati.

L’analisi della bontà di adattamento di un modello (goodness-of-fit, GOF) utilizza essenzialmente due step:

a) verifica della bontà di adattamento globale

b) verifica che tutte le osservazioni contribuiscano al fit, almeno approssimativamente, nella stessa misura.

L’analisi della GOF dipende dai gradi di libertà del modello, che a loro volta dipendono dal numero delle osservazioni e dalla presenza o meno di variabili continue.

Se il modello contiene solamente predittori categorici, la variabile risposta viene analizzata per sottogruppi, definiti dall’insieme dei valori che i predittori possono assumere (pattern di covariate).

Invece se nel database vi sono n osservazioni e almeno un predittore misurato su scala continua, ci si attende che al crescere delle osservazioni i pattern di covariate crescano allo stesso modo.

Nel modello logistico i valori stimati attesi per ogni pattern di covariate j, con mj osservazioni, sono dati dalla seguente espressione:

jjj my ˆ

dove:

jxje

ˆ

1

1ˆ

Nei modelli GLM, come il modello logistico, lo scostamento tra valori osservati e stimati attesi può essere valutato mediante il calcolo della devianza o deviance.

Posto l(0) la log verosimiglianza del modello fittato e l(F) è la log verosimiglianza del modello saturato (Full model), la devianza D viene ottenuta dalla loro differenza moltiplicata per -2.

j

jjjjj ymyl ˆ1lnˆln0

j j

jjj

j

jj m

yym

m

yyFl 1lnln

2 ln lnˆ ˆ1j j j

j j jj j j j j

y m yD y m y

m m

Tale statistica segue asintoticamente una distribuzione 2 con valore atteso pari al numero dei pattern di covariate meno il numero dei parametri nel modello (compresa l’intercetta).

Un test per la GOF può quindi essere formulato utilizzando tale proprietà della deviance, che risulta soddisfatta solo per modelli con predittori categorici.

Inoltre se il numero delle parametri è elevato rispetto a quello dei predittori (ovvero vi sono celle con valori attesi inferiori a 5), l’approsimazione alla distribuzione 2

potrebbe non essere valida.

Un altro test per la GOF analogo alla deviance è basato sui residui di Pearson:

jjj

jjjj

m

myr

ˆ1ˆ

ˆ

La sommatoria su j di tali residui al quadrato segue anch’essa asintoticamente una distribuzione 2 e può essere impiegata per testare la GOF del modello.

Anche per tale test valgono le stesse limitazioni del test sulla deviance.

In presenza di variabili continue nel modello oppure di un elevato numero di predittori, un test alternativo per la GOF è stato proposto da Hosmer e Lemeshow (1980).

Tale procedura consiste nel raggruppare le probabilità stimate per ogni soggetto in classi definite sulla base della distribuzione dei percentili.

In genere si utilizzano i decili della distribuzione, detti “decili di rischio”, oppure, meno frequentemente, valori di probabilità stimata fissati a priori (ad es: a step di 0.1).

Il test consiste nel calcolare un 2 di Pearson per la tabella di contingenza ottenuta dall’aggregazione dei dati in colonne corrispondenti, ad esempio, ai decili di rischio e nelle 2 righe, corrispondenti ai due outcome (0 o 1).

g

k KKk

Kkk

n

nyC

1

2

1ˆ

kc

j k

jjk n

m

1

Con:

E mj è il numero di osservazioni entro ogni pattern di covariate j

Infine, un’altra statistica ampiamente utilizzata per la valutazione della GOF è l’area sotto la curva ROC (Area Under the Curve, AUC).

Tale curva viene ottenuta utilizzando i valori stimati attesi del modello logistico e classificando i soggetti come appartenenti alla categoria di variabile risposta y=1 (corrispondente ai Casi) per valori superiori a valori selezionati.

In genere si utilizzano tutti i valori stimati attesi di probabilità (oppure in modo del tutto equivalente, del logit).

La curva viene ottenuta plottando i valori di sensibilità e 1- specificità ottenuti da ogni valore soglia.

In pratica, se una curva ROC presenta un’area elevata (tipicamente al di sopra di 0.8) vi è evidenza di un buon fitting del modello.

Valori di AUC intorno a 0.5 indicano una pessima GOF.

Infatti si dimostra (in modo non banale!) che l’AUC stima la probabilità che un soggetto appartenente al gruppo dei casi presenti valori di probabilità stimata attesa dal modello superiori ad un controllo.

Quindi se il modello presenta un’ottima GOF la grande maggioranza dei casi presenterà valori di probabilità stimata attesa elevati e AUC tenderà a 1.

Se il modello invece presenta una pessima GOF, i casi e i controlli presenteranno mediamente lo stesso valore di probabilità stimata attesa e AUC tenderà a 0.5.

Ricerca dei punti influenti (outliers)

Le statistiche di bontà di un modello dovrebbero essere sempre affiancate da diagnostiche di regressione, che mirano ad identificare le unità statistiche che possono aver influenzato l’esito dell’analisi (“punti influenti” o outliers).

Tra queste misure una delle più utilizzate è il leverage, che stima il peso che ogni osservazione esercita sul suo valore atteso.

Nel modello logistico Pregibon nel 1981 ha dimostrato che il leverage hj presenta la seguente relazione con i valori attesi .

jjjjjjjj bxWXXxxxmh 1ˆ1ˆ

jjj xWXXxb 1

Con:

Si noti che il termine rappresenta la stima della varianza degli yi

Andamento del leverage in funzione dei valori attesi da un modello di regressione logistica (dati simulati):

'1500 cases, 1500 controls'

leve

rag

e

Pr(caco).02254 .981144

.000431

.000977

jjjjjjj bWXXxxxmh 1ˆ1ˆ

Risulta quindi conveniente dividere il leverage per la quantità j, che ne provoca un abbassamento in corrispondenza degli estremi dell’intervallo di probabilità attese stimate.

Altre diagnostiche di regressione utilizzano la delezione delle singole osservazioni o pattern di covariate e ne stimano l’effetto sulle statistiche di GOF.

Una statistica molto utilizzata è la variazione in deviance (D) in corrispondenza dell’eliminazione di un determinato pattern di covariate:

Un plot di D verso i corrispondenti valori attesi permette di identificare visivamente punti con valore più elevato degli altri, che potrebbero aver influenzato il fit del modello.

Hosmer e Lemeshow suggeriscono inoltre di identificare tutti i punti con eventuale valore di D superiore a 4.

In modo del tutto analogo anche la delezione dei residui di Pearson può essere utilizzata per ricercare eventuali punti influenti.

Infine un’altra statistica molto utilizzata è la distanza di Cook che misura la differenza per ogni coefficiente tra il valore stimato nel modello con tutte le osservazioni e quello che si otterrebbe eliminando un determinato pattern di covariate.