Regression Binomiale e Logistica PDF

LEZIONE 12: MODELLI CHE RIENTRANO IN GLM ** REG LIN NORMALE (file 10) REGRESSIONE LINEARE BINOMIALE Analizzare l’influenza di una o più vv.ii. su una v.d. dummy significa chiedersi in che misura le prime incidono...

LEZIONE 12: MODELLI CHE RIENTRANO IN GLM ** REG LIN NORMALE (file 10) REGRESSIONE LINEARE BINOMIALE Analizzare l’influenza di una o più vv.ii. su una v.d. dummy significa chiedersi in che misura le prime incidono sulla frequenza relativa di valori ‘1’ nella v.d. Vediamo un esempio su dati Istat 2009, giovani 20-39enni. Y indica chi «vive fuori dalla famiglia di origine», 1=sì , 0=no. La numerosità campionaria è pari a 145'952; di questi, 74'532 hanno risposto sì. La media, E(Y), è la quota (frequenza relativa) di intervistati fuori dalla famiglia di origine: 0.511 (o 51.1%). Pertanto: Se il campione è casuale, la frequenza relativa E(Y) si interpreta come probabilità. Nel nostro esempio, la probabilità di estrarre dalla popolazione di 20-39enni nel 2009 un giovane che abita fuori dalla famiglia di origine (FFO) è: P(Y=1) = pigrecoi = 0.511 L’equazione di regressione della slide precedente può perciò essere riscritta come segue: Il modello di Regressione Lineare Binomiale è più conosciuto come Linear Probability Model. ES 1: Consideriamo dapprima solo i primi 20 casi in matrice. Stimiamo ora la dipendenza di v1 (v.d.) da v2 (v.i.) con la regressione lineare. Se vogliamo stimare l’impatto netto della condizione occupazionale sull’indipendenza abitativa dobbiamo tenere sotto controllo altre variabili antecedenti e concomitanti rilevanti. Passiamo quindi a un LPM multiplo. Calcoliamo ora, sulla base dell’output della slide precedente, le probabilità di uscita da casa per alcuni profili. Conseguenze di una specificazione lineare quando Y è dicotomica - alcuni valori predetti ( i) possono essere ‘out of range’ -la distribuzione degli errori non è normale; gli errori di previsione commessi dal modello (Y-Ŷ) possono assumere solo due valori: ad es., in un modello con un solo regressore, per xi=1 gli errori sono 1- - (se Y=1) oppure - - (se Y=0); le stime intervallari e i test parametrici sono meno affidabili (specie in campioni piccoli); - c’è eteroschedasticità: Var(e| xi)=s2 i=(beta0+beta1+xi)[1-(beta0+beta1+ xi)]; pertanto, le stime prodotte con OLS sono ancora corrette ma non efficienti; - l’incremento di pii non è costante in tutto il campo di variazione delle vv.ii.: l’impatto di queste è sottostimato per i valori centrali e sovrastimato per quelli estremi. Le due principali incongruenze del Linear Probability Model REGRESSIONE LOGISTICA BINOMIALE La regressione logistica binomiale, che si usa quando la v.d. y può assumere solo due valori Per ovviare ai limiti del LPM, colleghiamo i e il predittore lineare con una diversa link function, non lineare. Quella più usata è la RLB. L’andamento nella prima parte è quasi esponenziale poi, dopo il punto di flesso, rallenta e la curva diventa asintotica al valore 1. Nella prima equazione osserviamo che la distribuzione teorica di riferimento della variabile aleatoria Yi è quella bernoulliana (o binomiale con n=1)3 e non più, come nel caso della regressione lineare, quella normale. Quando la v.d. y è dicotomica, infatti, anche la corrispondente variabile aleatoria Yi potrà assumere solo una coppia di valori discreti, convenzionalmente indicati con 0 e 1. Se stabiliamo che πi è la probabilità che Yi assuma il valore 1, la probabilità che assuma valore 0 è uguale a 1-πi. La distribuzione bernoulliana ha quindi solo un parametro indipendente, πi. Forma moltiplicativa Nella seconda equazione, osserviamo che il collegamento tra la probabilità πi e il predittore lineare η! non è più, come era nella regressione lineare normale, la funzione identità ma una funzione logistica4, indicata con la lettera greca Λ. Questa funzione vede il predittore lineare come esponente della costante di Nepero e (2.718…); essa restituisce, a seconda dei valori assunti dal predittore lineare η! compresi nel dominio [-∞ , +∞] una probabilità πi compresa nell’intervallo tra 0 e 1 (0≤πi≤1). L’equazione in forma moltiplicativa, può essere riscritta in forma adittiva ed è la link function La terza equazione descrive ηi come l’esito di una somma lineare definita dai coefficienti β da stimare Forma additiva: interpretazione semplice del membro a destra ma non di quello a sinistra La quantità a sinistra di quest’ultima equazione è detta logit (logarithmic odds unit) ossia un’unità di misura espressa dal log naturale dell’odds, termine che indica il rapporto tra la probabilità di accadimento di un evento e il suo contrario. La regressione logistica binomiale non consente di calcolare l’impatto di una v.i. X su Pr(Y=1) in modo univoco; l’effetto marginale, ossia l’effetto su Pr(Y=1) di una variazione unitaria di X non è costante, ma dipende: a) dal livello della X (salvo che sia dicotomica); b) dalle combinazioni dei valori assunti dalle altre vv.ii. del modello. Metodo stima parametri β: Non essendo lineare la relazione tra πi e x, il predittore lineare, la stima dei parametri non può essere effettuata con il metodo dei minimi quadrati. Il metodo usato è quello della massima verosimiglianza (maximum likelihood estimation, ML). Questo metodo permette di ottenere i coefficienti di regressione per i quali la probabilità di ottenere i dati che abbiamo osservato nella v.d. è massima. A differenza di OLS non fornisce una formula analitica per calcolare le stime ma è implementato da un algoritmo iterativo di ottimizzazione numerica. Lo stimatore massimizza la seguente funzione di log-verosimiglianza (LL): ln[flogit ( 0 … k; Y1…Y n|X1n…Xkn)] Semplificando: tra tutti i valori che le incognite possono assumere l’algoritmo individua quelli che più probabilmente hanno prodotto i valori osservati nel campione. Con campioni sufficientemente ampi (100 casi più 10 per ogni parametro oltre i primi due) il metodo ML produce stime dei parametri corrette, normalmente distribuite ed efficienti. Esempio : Riprendiamo il solito esempio delle ore di studio (X) e dell’esame di matematica; questa volta la variabile dipendente è dicotomica (1/0) «esame superato/non superato» e il modello è una regressione Nel diagramma sono rappresentati i casi con le seguenti coordinate: sulle ascisse le ore aggiuntive di studio per la preparazione all’esame; sulle ordinate gli esiti (successo o insuccesso). Allo scatterplot è stata sovrapposta la curva sigmoidale corrispondente alla funzione logistica A(ηi ) La v.d. y definisce l’appartenenza del singolo studente a uno dei due gruppi (promosso o respinto). Obiettivo: stimare la probabilità di uno studente di appartenere a uno dei due gruppi, a partire dalle ore dedicate alla preparazione dell’esame. Per semplificare, consideriamo binaria anche la variabile indipendente, riaccorpando le ore aggiuntive in due modalità in base al valore mediano della distribuzione (60 ore) e osserviamone la distribuzione congiunta con l’esito all’esame (tab. 8).  Per prima cosa, osserviamo la percentuale di successo all’esame indipendentemente dalle ore dedicate allo studio % Promossi=108 200 ∗ 100 = 54.0% Ipotizzando che il campione sia casuale, possiamo usare la frequenza relativa come stima della probabilità π di essere promossi indipendentemente dalle ore di studio: π =0.54 e, poiché la probabilità varia in un intervallo (0,1), calcoliamo la probabilità di essere respinti come complemento a 1: 1 −π=0.46 Le medesime informazioni possono essere espresse in modo alternativo attraverso il rapporto di probabilità, l’odds. Quando la proporzione dei successi (Y = 1) è uguale a quella degli insuccessi (Y = 0), l’odds è evidentemente pari a 1. Quando il numero dei successi è superiore a quello degli insuccessi, l’odds assume maggiori di 1 e nel caso opposto, compresi tra 0 e 1. Nel nostro esempio: Il valore dell’odds indica che ogni 10 studenti respinti, 12 superano l’esame. Con la stessa formula possiamo calcolare l’odds di superare l’esame tra gli studenti che hanno dedicato più di 60 ore a prepararlo: il risultato è pari a 4.1, ossia per ogni studente che fallisce 4 hanno successo. Il medesimo rapporto calcolato tra chi ha dedicato meno ore a studiare è pari a 0.7, il che indica che ogni 7 studenti promossi 10 vengono respinti. Utilizzando il rapporto tra gli odds, è possibile infine calcolare un indice di associazione tra le due variabili considerate, che si chiama odds ratio (OR). Questo indice, come vedremo, è utile per interpretare i coefficienti stimati. L’odds ratio è un rapporto di probabilità condizionate che si interpreta come segue: le chance di superare l’esame anziché non superarlo tra chi fa più di 60 ore aggiuntive di studio supera in media di oltre 5 volte le chance di superarlo anziché non superarlo tra chi fa meno di 60 ore di studio aggiuntive. I logit, presentati nell’ultima riga della tabella, sono infine un altro modo di esprimere la medesima informazione contenuta negli odds. Il logit è il logaritmo naturale (ln) dell’odds. Il campo di variazione del logit è [-∞,+∞], invece quello dell’odds è [0, +∞]. Si noti che, quando l’odds è compreso tra 0 e 1, il logit ha segno negativo, mentre quando è maggiore di 1 ha segno positivo. Nel caso di odds pari a 1 (equiprobabilità), il logit assume valore 0. Es 2 + grado interesse: Riprendiamo l’esempio precedente e aggiungiamo al modello una variabile indipendente, il grado di interesse mostrato dagli studenti per la materia: basso, medio e elevato. Possiamo inserire nel modello un ulteriore regressore anche se la numerosità campionaria (200 casi) non è particolarmente elevata: per la regressione logistica multipla servono almeno 100 casi di base più 10 per ciascun regressore aggiuntivo Metodo 1 per interpretare i risultati regressione logistica →Sostituiamo ora i valori dei parametri stimati all’interno dell’equazione usando ogni combinazione delle due vv.ii. per calcolare le probabilità predette di essere promosso Y «esito esame» (1: ‘positivo’, 0: ‘negativo’) X «ore studio aggiuntivo» (1: ‘>60h.’, 0: ‘ 60h.’) C «interesse per la materia» (var. ordinale ricodificata in tre regressori indicatori per ‘basso’, ‘medio’, ‘elevato’). Per calcolare l’effetto prodotto da X su Pr(Y=1) dobbiamo: a) sceglierne due livelli da confrontare; b) fissare le variabili di controllo a valori di interesse, C = c. Successivamente calcoliamo: - Pr(Y|X=x; C = c ) - Pr(Y|X=x+ x; C = c ) - infine la differenza tra le due probabilità  Probabilita % di essere promosso all’esame per le 6 combinazioni delle vv.ii: il modello utilizza una funzione logistica per trasformare le variabili indipendenti (interesse/ore studio) in una probabilità compresa tra 0 e 1 Es Interesse basso/fino a 60ore →calcolo del valore della funzione logistica: z = -1.06 + exp(-1.06) Trasformo il valore in una probabilità: p = 1 / (1 + exp(-z)) Il calcolo della probabilità nella tabella si basa su un modello di reg log stimato, i coefficienti del modello sono (-1.06, exp(-1.06)) e rappresentano l’effetto di ogni vv indipendente sulla probabilità di essere promossi. L’impatto netto della v.i. X «ore studio aggiuntive» su Pr(Y=1, cioè esito esame positivo) cambia al variare della variabile di controllo (Interesse). La differenza di probabilità di successo all’esame dovuta alle ore di studio aggiuntivo è elevata tra coloro che hanno un interesse basso (+41.1 p.p.) e per chi ha un interesse medio (+39.9 p.p.) mentre è molto inferiore per chi ha un interesse elevato (+7.5 p.p.). Quindi nei tre livelli di interesse per la materia - tra quanti si sono impegnati di più e chi ha dedicato meno tempo alla preparazione dell’esame, si osserva che quando la materia interessa poco, l’impegno profuso nello studio può davvero fare la differenza (41 punti di differenza nella probabilità di superare l’esame anziché essere respinti). Avere un interesse medio, l’impegno anche fa la differenza 39.9 di superare l’esame invece di essere respinti. Per chi ha un interesse elevato, solo il 7.5 punti di differenza di superare l’esame se ci si impegna di più. Metodo 2 interpretazione regressione logistica →Lettura parametri: Nella tab 19 sono riassunti i risultati della stima della probabilità di superare l’esame dato il numero di ore aggiuntive dedicate allo studio e l’interesse per la materia. 1. vedere la significatività statistica dei parametri stimati: tutti i parametri sono stat significativi eccetto, quello associato a interesse medio (0.06) 2. Interpretiamo i parametri facendo rif alla colonna in cui i logit sono trasformati in odds ratio: i coefficienti di regressione stimati vanno considerati al netto dell’influenza esercitata dagli altri regressori inclusi nel modello. Diremo quindi che, a parità di interesse mostrato per la materia, la probabilità di superare l’esame anziché essere respinti tra chi ha studiato più di 60 ore è superiore di quasi sei volte rispetto a chi ha studiato meno ore. Ma anche l’interesse per la materia influenza la probabilità di superare o meno l’esame. A parità di tempo dedicato allo studio, infatti, un elevato interesse per i contenuti studiati (rispetto a uno basso) aumenta di 27 volte la probabilità di superare l’esame. 1 Nella prima colonna (“b stimati”) sono riportate le stime di massima verisimiglianza dei parametri, in quella successiva le stime degli errori standard. Nelle quattro colonne successive si trovano le informazioni per valutare la significatività statistica dei parametri stimati: prima attraverso la stima intervallare e poi con la statistica-test W e l’associato p-value. Abbiamo per completezza incluso entrambi i modi per fare inferenza sui parametri stimati anche se nella pratica l’analista solitamente opta per uno dei due. Nell’ultima colonna, infine, sono riportati gli odds ratio. Va precisato che la trasformazione dei b in odds ratio mediante la loro esponentazione non aggiunge alcuna informazione, è solo un altro modo di esprimere il rapporto tra v.d. e v.i. Gli odds ratio non sono, in altri termini, analoghi ai b standardizzati. Analogamente a quanto visto per la regressione lineare multipla, anche per la logistica si possono standardizzare i parametri. Riportiamo i b stimati del modello con cui abbiamo calcolato i dati della tabella precedente. L’impatto delle ore di studio aggiuntive sul logit(πi), al netto dell’interesse per la materia è b1=1.76. L’interpretazione di questo b non è diversa rispetto a un b di regressione parziale del modello lineare. Tuttavia, in questo caso la v.d. è logit(πi), ossia una quantità non intuitiva. Invece data la non linearità di Ʌ(ηi) l’impatto b di X1 su πi non è un valore costante. Calcolo valori predetti dal modello ESERCIZIO 3: Fissiamo le variabili di controllo: uomo, res. NO, età20_24, scuola obbligo Consideriamo un OCCUPATO: η1= b0 (intercetta)+ b1(coeff di occupato)=-2.647+.446= -2.201 Fare es per : - Donna, res. Centro, eta25_29, post diploma OCCUPATA η2= b0 (intercetta)+ b2(coeff di occupato)=-2.647+.911= -1.736 tra le donne residenti nel Centro Italia, 25_29enni, con titolo postdiploma, la propensione all’indipendenza abitativa è [maggiore?/minore?] tra le occupate rispetto alle non occupate di …[quanti?]. punti percentuali - NON OCCUPATO: tra i maschi residenti nel N-O, 20_24enni, con scuola dell’obbligo, la propensione all’indipendenza abitativa è [maggiore?/minore?] per gli occupati rispetto ai non occupati di …[quanti?].. punti percentuali (Vedere soluzione) Gli exp(b) sono Odds Ratio: quanto varia l’odds(probabilità che si verifichi l’evento) di Y=1 quando X cresce di una unità: Odds(x+1)/Odds(x). Variano tra 0 e +infinito, il valore 1 indica totale indipendenza da X della P(Y). I valori negativi del logit, una volta trasformati in odds ratio, variano nell’intervallo [0,1] mentre quelli positivi variano nell’intervallo [0, +∞]. Al netto delle variabili di controllo, la propensione relativa degli occupati a uscire dalla FO, anziché rimanervi, è circa una volta e mezza di quella dei non occupati (o anche: la propensione relativa a uscire degli occupati è il 56% in più dei non occupati). Per facilitare il confronto di O.R. compresi tra 0 e 1 (relazione inversa) con O.R. >1 (relazione diretta) si calcolano i reciproci dei primi e si aggiunge il segno meno. BONTÀ ADATTAMENTO AL MODELLO La valutazione della bontà di adattamento del modello di regressione logistica ai dati segue una logica solo in parte simile a quella in uso nel caso della regressione lineare normale. Ciò è dovuto a varie ragioni: ad esempio al fatto che la variabile predetta e quella osservata non solo presentano, ovviamente, valori diversi, ma hanno una scala diversa (una è dicotomica e l’altra è continua); la varianza di una variabile dipendente dicotomica è logicamente diversa da quella di una variabile dipendente continua ecc.: pertanto, misure di varianza spiegata esattamente corrispondenti a R2 non sono realizzabili, se ne possono al più pensare delle approssimazioni. Indici di bontà 1. Rapporto di verosimiglianza G2 (talvolta indicato come 𝜒2) in quanto segue la distribuzione teorica del chi-quadrato; basato sul confronto di due equazioni: la prima equazione descrive il modello nullo (M0) in cui si stima solo il parametro b0, la seconda descrive il modello M1 che include anche i regressori. 2. Poi ci sono gli Indici indicati come Pseudo-R2 e possono essere interpretati, analogamente a quanto si fa con il coefficiente di determinazione nella regressione lineare, cioè come misura della variabilità della v.d. spiegata dai regressori. ◼ R2 di McFadden, noto anche come likelihood ratio index, è una misura di adattamento del modello che usa la funzione di log verosimiglianza (log Likelihood), mettendo a confronto il modello pieno, con tutti i regressori (M1) con il modello intercetta (M 0). La formula per calcolarlo è la seguente: Per costruzione, l’indice ha valore minimo 0(caso in cui M1=M2) ma il suo massimo non raggiunge mai il valore di 1. Il valore aumenta con l’inserimento di nuovi regressori nel modello, senza signidicare un miglioramento stat significativo nella bontà di adattamento. Per questo Ben-Akiva e Lerman hanno proposto una variante adjusted, che penalizza i modelli non parsimoniosi ◼ Un altro Pseudo-R2 molto utilizzato è stato proposto da Cox e Snell (1989) e si basa sempre sul confronto tra M1 e M0 3. Hosmer-Lemeshow Goodness-Of-Fit test CAPACITÀ PREDITTIVA DI UN MODELLO LOGBIN Per capacità predittiva si intende in che misura il modello specificato è in grado di prevedere l’appartenenza dei casi osservati alle modalità della vv dipendente. - A tal fine, in primo luogo, vanno individuate delle vv.ii. di cui si ipotizza una buona capacità predittiva della v.d.: in altri termini, che permettano di classificare correttamente il numero più elevato possibile di casi sulla v.d - Dopo aver stimato i parametri del modello, costruiamo una matrice contenente i casi che il modello prevede correttamente o no. - Preliminarmente, occorre rendere confrontabili i valori osservati della variabile dummy Y con i valori predetti i che sono continui (tra 0 e 1). - Fissato un valore-soglia c, discretizziamo i in una variabile dummy Ŷ: se i  c allora Ŷ=1 , in caso contrario Ŷ=0. La soglia c è detta anche cut off point. Il valore di default dei software statistici è c=0.5 ; in tal caso: Ŷ=0 | i < 0.5 e Ŷ=1 | i> 0.5. Chiamiamo Ŷ=0 «caso negativo», Ŷ=1 «caso positivo», lo stesso facciamo per gli Y. Incrociando i valori predetti e osservati di y avremo 4 situazioni: due per i casi correttamente classificati, i «veri positivi» e i «veri negativi» due per quelli classificati erroneamente, i «falsi positivi» e i «falsi negativi» Matrice di confusione: è una tab delle classificazioni, per costruirla è necessario tradurre in dicotomie le probabilità predette dal modello. Torniamo al cut off point: scegliere un punto di taglio dei valori di probabilità c=0.5 significa dare uguale importanza ai due tipi di errori (FN e FP). Talvolta uno dei due è ritenuto più ‘grave’ o ‘costoso’ dell’altro. Ad es. scegliere c=0.2 equivale a considerare gli errori di I tipo (FN) quattro volte più costosi di quelli di II tipo (FP). In generale, dato k, il rapporto tra i costi di FN e FP, c=1/(1+k) 77+74/77+74+15+34 *100= 151/200=0.755*100= 75.5% Ciò sarebbe però improprio, perché includeremmo una parte di casi concordanti per puro effetto del caso. La percentuale di studenti promossi è infatti, nel nostro esempio, pari al 54%: quindi, oltre la metà degli studenti verrebbero classificati correttamente banalmente assegnando tutti gli studenti alla modalità ‘Promossi’. Alcuni indici di potere predittivo di un modello di regressione LogBin La buona ragione per non rinunciare al LPM è che i suoi parametri sono più facilmente interpretabili di quelli ottenuti con la regressione logistica. Peraltro, il numero di casi con probabilità predette come ‘insensate’ è spesso assai contenuto. Quando l’uso della regressione è a fini esplorativi, si possono per semplicità commentare i parametri del LPM a condizione che: i segni, le significatività e l’ordine di importanza dei b del LPM corrispondano a quelli ottenuti con la regressione logbin; la % di valori ŷi “out of range” nel LPM sia < 5%. Cenni alla Regressione Logistica Multinomiale È un’estensione della regressione LogBin. In quest’ultima abbiamo una equazione per un solo logit che pone a confronto la probabilità del verificarsi di Y con quella del suo non verificarsi: Nella multinomiale, con una Y a j modalità, abbiamo j-1 logit ed equazioni, tenendo una modalità come riferimento (pivot). Ad es. con 3 modalità si fanno 2 confronti con la pivot. Esempio: Quanto influiscono alcune variabili sulla probabilità di aderire a uno di tre orientamenti religiosi: j=1 «cattolico senza riserve» (cat. rif.), 2 «cattolico con riserve», 3 «credente senza appartenenza». I regressori sono quattro: X1 «età» (continua), X 2 «femmina» (dummy), X 3 «Italia centrale» (dummy), X 4 «Italia meridionale» (dummy). Possiamo poi mettere a confronto le probabilità predette per differenti combinazioni di valori delle variabili esplicative. Es.: come cambia la probabilità di aderire a tre diversi orientamenti religiosi al variare dell’età (X) per le donne residenti in Italia meridionale. Si nota che tra le donne del Sud al crescere dell’età aumenta la probabilità di essere CSR. L’età invece influisce (inversamente) poco sulla probabilità di essere CSA e pochissimo su quella di essere CCR.