LABORATORIO DI STATISTICA AZIENDALE

1. 5. Modello di Regressione logistica con R LABORATORIO DI STATISTICA AZIENDALE Enrico Properzi - enrico.properzi3@unibo.it A.A. 2010/2011

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”)

3. 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 dell’università (rank) influiscano sull’ammissione alla scuola di specializzazione post-laurea. La variabile risposta (admit) è quindi una variabile binaria (0/1) graduate <- read.csv("graduate.csv", header=T) str(graduate) 'data.frame': 400 obs. of 4 variables: $ admit: int 0 1 1 1 0 1 1 0 1 0 ... $ gre : int 380 660 800 640 520 760 560 400 540 700 ... $ gpa : num 3.61 3.67 4 3.19 2.93 3 2.98 3.08 3.39 3.92 ... $ rank : int 3 3 1 4 4 2 1 2 3 2 ... attach(graduate) > table(rank) rank 1 2 3 4 61 151 121 67

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

5. Costruiamo un modello logit con un solo regressore (gre) utilizzando la funzione glm mod1 <-glm(admit~gre,family=binomial(link="logit")) > summary(mod1) Call: glm(formula = admit ~ gre, family = binomial(link = "logit")) DevianceResiduals: 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 (Dispersionparameterforbinomial family takentobe 1) Nulldeviance: 499.98 on 399 degreesoffreedom Residualdeviance: 486.06 on 398 degreesoffreedom AIC: 490.06 Numberof Fisher Scoringiterations: 4

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 <-glm(admit~gre+gpa+as.factor(rank), + family=binomial(link="logit"), na.action=na.pass) summary(modlogit) Call: glm(formula = admit ~ gre + gpa + as.factor(rank), family = binomial(link = "logit"), na.action = na.pass) Deviance Residuals: Min 1Q Median 3Q Max -1.6268 -0.8662 -0.6388 1.1490 2.0790

7. 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

8. L’incremento di un’unità della variabile gre determina l’aumento di 0.002 del log odds della variabile admit L’incremento di un’unità della variabile gpa determina l’aumento di 0.804 del log odds della variabile admit Le variabili dummy associate al rank hanno un significato leggermente diverso. Ad esempio, aver frequentato un’università con rank 2 riduce il log odds della variabile admit di 0.675 rispetto all’aver frequentato un’università 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 l’aumento di un’unità della variabile gpa determina l’aumento di un fattore pari a 2.23 dell’odds di essere ammessi alla scuola di specializzazione.

9. CONFRONTO TRA MODELLI Considerando due modelli: Modello completo (C) con k+r variabili esplicative Modello ridotto (R) con k variabili esplicative Lkr : verosimiglianza relativa al modello stimato C Lk : 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. L’ipotesi 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.

10. > L0 <- logLik(mod1) log-verosimiglianza modello semplice > L0 'log Lik.' -243.0281 (df=2) > L1 <- logLik(modlogit) log-verosimiglianza modello completo > L1 'log Lik.' -229.2587 (df=6) > L01 <- as.vector(-2*(L0-L1)) test LR > L01 [1] 27.53865 > df<- attr(L1,"df")-attr(L0,"df") calcolo dei gradi di libertà > df [1] 4 > pchisq(L01, df, lower.tail=F) p-value [1] 1.546749e-05 Il p-value piccolo porta a rifiutare l’ipotesi nulla -> il modello completo è preferibile!

11. Caso di studio: In un’indagine condotta nel 1974-75 a ciascun intervistato era stato chiesto se era d’accordo 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. L’obiettivo è valutare se le risposte degli uomini e delle donne differiscono e quanto l’educazione influisce sulle risposte. donne<- read.table("womenrole.txt", header=T) str(donne) 'data.frame': 42 obs. of 4 variables: $ education: int 0 1 2 3 4 5 6 7 8 9 ... $ sex : Factor w/ 2 levels "F","M": 2 2 2 2 2 2 2 2 2 2 ... $ agree : int 4 2 4 6 5 13 25 27 75 29 ... $ disagree : int 2 0 0 3 5 7 9 15 49 29 ... Abbiamo quindi due variabili esplicative: sex e education

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<- glm(cbind(agree, disagree) ~ sex + education, family=binomial()) > 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 <2e-16 *** sexM 0.01145 0.08415 0.136 0.892 education -0.27062 0.01541 -17.560 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 451.722 on 40 degrees of freedom Residual deviance: 64.007 on 38 degrees of freedom AIC: 208.07 Number of Fisher Scoring iterations: 4

13. Dall’output risulta evidente che la variabile education gioca un ruolo significativo nel prevedere se un individuo sia d’accordo o meno con l’affermazione oggetto dell’indagine. La variabile sex sembra invece non essere importante. Proviamo ora a verificare se c’è un’interazione tra i due regressori: > women2 <- glm(cbind(agree, disagree) ~ sex*education, family=binomial()) > 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)

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 d’accordo con l’affermazione rispetto agli uomini. All’aumentare degli anni di scolarizzazione superano i 10 la situazione si ribalta.