LABORATORIO DI STATISTICA AZIENDALE

5. Modello di Regressione logistica con R Enrico Properzi A.A. 2010/2011 [email protected]

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

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

> 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

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

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

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

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

CONFRONTO TRA MODELLI

Considerando due modelli:   Modello completo (C) con k+r variabili esplicative Modello ridotto (R) con k variabili esplicative L kr L k : verosimiglianza relativa al modello stimato C : 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 parametri che aggiunti Pertanto se si vuole nel rifiutiamo verificare modello questa ipotesi è completo che siano i allora coefficienti congiuntamente degli r nulli.

il modello aggiunge qualcosa di significativo al modello più semplice.

completo

> L0 <- logLik(mod1) > L0 'log Lik.' -243.0281 (df=2) > L1 <- logLik(modlogit) > L1 'log Lik.' -229.2587 (df=6) > L01 <- as.vector(-2*(L0-L1)) > L01 [1] 27.53865

> df<- attr(L1,"df")-attr(L0,"df") > df [1] 4 log-verosimiglianza modello semplice log-verosimiglianza modello completo test LR calcolo dei gradi di libertà > pchisq(L01, df, lower.tail=F) [1] 1.546749e-05 p-value Il p-value piccolo porta a rifiutare l’ipotesi nulla -> il modello completo è preferibile!

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

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

Dall’output risulta evidente che la variabile significativo nel prevedere se un individuo sia d’accordo o meno con l’affermazione oggetto dell’indagine.

education gioca un ruolo 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)

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

significativo.

risulta essere altamente 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.