Slide 1

5. Modello di Regressione logistica con R Enrico Properzi - enrico.properzi3@unibo.itenrico.properzi3@unibo.it A.A. 2010/2011

Slide 2

I modelli in cui la variabile dipendente dicotomica rientrano come caso particolare dei modelli di regressione generalizzata. In R la stima di questa tipologia di modelli viene realizzata per mezzo del comando glm (formula, family = gaussian, data, weights, subset, na.action, ) Dove il parametro family pu assumere le seguenti specifiche: binomial(link= logit) gaussian(link=identity) Gamma(link=inverse) inverse.gaussian(link=1/mu^2) poisson(link=log) quasi (link=identity, variance=constant) quasibinomial(link=logit) Quasipoisson(link=log) Nel caso particolare di variabile dipendente dicotomica si utilizza il parametro: Family= binomial(link=logit)

Slide 3 table(rank) rank 1 2 3 4 61 151 121 67">

Caso di studio: Un ricercatore interessato a come alcune variabili, tra cui il punteggio di laurea (GRE), la media degli esami (GPA) e il prestigio delluniversit (rank) influiscano sullammissione alla scuola di specializzazione post-laurea. La variabile risposta (admit) quindi una variabile binaria (0/1) graduate table(rank) rank 1 2 3 4 61 151 121 67

Slide 4

> table(admit) admit 0 1 273 127 > table(rank,admit) admit rank 0 1 1 28 33 2 97 54 3 93 28 4 55 12

Slide 5 summary(mod"> summary(mod1) Call: glm(formula = admit ~ gre, family = binomial(link = "logit")) Deviance Residuals: Min 1Q Median 3Q Max -1.1623 -0.9053 -0.7547 1.3486 1.9879 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.901344 0.606038 -4.787 1.69e-06 *** gre 0.003582 0.000986 3.633 0.00028 *** --- Signif. codes: 0 *** 0.001 ** 0.01 * 0.05. 0.1 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 499.98 on 399 degrees of freedom Residual deviance: 486.06 on 398 degrees of freedom AIC: 490.06 Number of Fisher Scoring iterations: 4"> summary(mod" title="Costruiamo un modello logit con un solo regressore (gre) utilizzando la funzione glm mod1 summary(mod">

Costruiamo un modello logit con un solo regressore (gre) utilizzando la funzione glm mod1 summary(mod1) Call: glm(formula = admit ~ gre, family = binomial(link = "logit")) Deviance Residuals: Min 1Q Median 3Q Max -1.1623 -0.9053 -0.7547 1.3486 1.9879 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.901344 0.606038 -4.787 1.69e-06 *** gre 0.003582 0.000986 3.633 0.00028 *** --- Signif. codes: 0 *** 0.001 ** 0.01 * 0.05. 0.1 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 499.98 on 399 degrees of freedom Residual deviance: 486.06 on 398 degrees of freedom AIC: 490.06 Number of Fisher Scoring iterations: 4

Slide 6

Costruiamo ora un modello di regressione logistica pi completo, utilizzando tutti I regressori a disposizione. Il comando as.factor(rank) indica che la variabile rank viene trattata come un fattore (var. categorica) >modlogit

Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -3.989979 1.139951 -3.500 0.000465 *** gre 0.002264 0.001094 2.070 0.038465 * gpa 0.804038 0.331819 2.423 0.015388 * as.factor(rank)2 -0.675443 0.316490 -2.134 0.032829 * as.factor(rank)3 -1.340204 0.345306 -3.881 0.000104 *** as.factor(rank)4 -1.551464 0.417832 -3.713 0.000205 *** --- Signif. codes: 0 *** 0.001 ** 0.01 * 0.05. 0.1 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 499.98 on 399 degrees of freedom Residual deviance: 458.52 on 394 degrees of freedom AIC: 470.52 Number of Fisher Scoring iterations: 4

Slide 8

Lincremento di ununit della variabile gre determina laumento di 0.002 del log odds della variabile admit Lincremento di ununit della variabile gpa determina laumento di 0.804 del log odds della variabile admit Le variabili dummy associate al rank hanno un significato leggermente diverso. Ad esempio, aver frequentato ununiversit con rank 2 riduce il log odds della variabile admit di 0.675 rispetto allaver frequentato ununiversit con rank 1. Si possono anche esplicitare I coefficienti ed interpretarli come odds- ratio: > exp(modlogit$coef) (Intercept) gre gpa as.factor(rank)2 0.0185001 1.0022670 2.2345448 0.5089310 as.factor(rank)3 as.factor(rank)4 0.2617923 0.2119375 Ora possiamo affermare che laumento di ununit della variabile gpa determina laumento di un fattore pari a 2.23 dellodds di essere ammessi alla scuola di specializzazione.

Slide 9

CONFRONTO TRA MODELLI Considerando due modelli: Modello completo (C) con k+r variabili esplicative Modello ridotto (R) con k variabili esplicative L kr : verosimiglianza relativa al modello stimato C L k : verosimiglianza relativa al modello stimato R Il metodo pi usato per confrontare pi modelli il metodo del rapporto delle verosimiglianze, ovvero il test LR. I confronti tra modelli si fanno costruendo rapporti di verosimiglianze o equivalentemente differenze tra log-verosimiglianze. Nel nostro caso vogliamo confrontare il modello con un solo regressore con quello con tutti i regressori a disposizione. Lipotesi nulla che si vuole verificare che i coefficienti degli r parametri aggiunti nel modello completo siano congiuntamente nulli. Pertanto se rifiutiamo questa ipotesi allora il modello completo aggiunge qualcosa di significativo al modello pi semplice.

Slide 10 df [1] 4 > pchisq(L01, df, lower.tail=F)p-value [1] 1.546749e-05 Il p-value piccolo porta a rifiutare lipotesi nulla -> il modello completo preferibile!">

> L0 L0 'log Lik.' -243.0281 (df=2) > L1 L1 'log Lik.' -229.2587 (df=6) > L01 L01 [1] 27.53865 > df df [1] 4 > pchisq(L01, df, lower.tail=F)p-value [1] 1.546749e-05 Il p-value piccolo porta a rifiutare lipotesi nulla -> il modello completo preferibile!

Slide 11

Caso di studio: In unindagine condotta nel 1974-75 a ciascun intervistato era stato chiesto se era daccordo o in disaccordo con la seguente affermazione: le donne dovrebbero occuparsi di mandare avanti la propria casa lasciando agli uomini il compito di mandare avanti il paese Le risposte sono riassunte nel dataset womenrole.txt. Lobiettivo valutare se le risposte degli uomini e delle donne differiscono e quanto leducazione influisce sulle risposte. donne

Slide 12

Per definire un modello di regressione logistica utilizzando la funzione glm dobbiamo specificare il numero di accordi e disaccordi come una matrice a due colonne che rappresenta la variabile risposta > women1 summary(women1) Call: glm(formula = cbind(agree, disagree) ~ sex + education, family = binomial()) Deviance Residuals: Min 1Q Median 3Q Max -2.72544 -0.86302 -0.06525 0.84340 3.13315 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 2.49793 0.18278 13.666

Dalloutput risulta evidente che la variabile education gioca un ruolo significativo nel prevedere se un individuo sia daccordo o meno con laffermazione oggetto dellindagine. La variabile sex sembra invece non essere importante. Proviamo ora a verificare se c uninterazione tra i due regressori: > women2 summary(women2) Call: glm(formula = cbind(agree, disagree) ~ sex * education, family = binomial()) Deviance Residuals: Min 1Q Median 3Q Max -2.39097 -0.88062 0.01532 0.72783 2.45262 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 3.00294 0.27238 11.025 < 2e-16 *** sexM -0.90474 0.36007 -2.513 0.01198 * education -0.31541 0.02365 -13.338 < 2e-16 *** sexM:education 0.08138 0.03109 2.617 0.00886 ** --- Signif. codes: 0 *** 0.001 ** 0.01 * 0.05. 0.1 1 (Dispersion parameter for binomial family taken to be 1)

Slide 14

Null deviance: 451.722 on 40 degrees of freedom Residual deviance: 57.103 on 37 degrees of freedom AIC: 203.16 Number of Fisher Scoring iterations: 4 Il termine di interazione sex*education risulta essere altamente significativo. Possiamo osservare che nel caso di pochi anni di scolarizzazione le donne manifestano una maggiore probabilit di essere daccordo con laffermazione rispetto agli uomini. Allaumentare degli anni di scolarizzazione superano i 10 la situazione si ribalta.