LABORATORIO DI STATISTICA AZIENDALE

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

2. Se la variabile risposta è una variabile di conteggio (e può assumere solo valori interi) e segue la distribuzione di Poisson e si vuole studiare la seguente relazione: E(yi|xi) = λi = exp(xi’ß) ossia la dipendenza di Yi da un insieme di variabili Xi (regressione di Poisson) si possono utilizzare i GLM

3. Caso di studio: Si vuole studiare se il numero di incidenti subiti da un insieme di barche dipende da alcuni regressori: type -> tipo di imbarcazione construction -> anno di costruzione operation -> periodo di operatività months -> mesi di servizio effettuati barche <- read.table("ships.txt", header=T) str(barche) 'data.frame': 34 obs. of 5 variables: $ type : Factor w/ 5 levels "A","B","C","D",..: 1 1 1 1 1 1 1 2 2 2 ... $ construction: Factor w/ 4 levels "1960-64","1965-69",..: 1 1 2 2 3 3 4 1 1 2 ... $ operation : Factor w/ 2 levels "1960-74","1975-79": 1 2 1 2 1 2 2 1 2 1 ... $ months : int 127 63 1095 1095 1512 3353 2244 44882 17176 28609 ... $ damage : int 0 0 3 4 6 18 11 39 29 58 ...

4. Costruiamo un primo modello in cui ipotizziamo che il numero dei danni dipenda soltanto dal tipo di imbarcazione: > mod_type <- glm(damage~type, family=poisson(link="log"), data=barche) > summary(mod_type) Call: glm(formula = damage ~ type, family = poisson(link = "log"), data = barche) Deviance Residuals: Min 1Q Median 3Q Max -4.6716 -2.2039 -0.5921 0.8632 4.0113 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 1.7918 0.1543 11.612 < 2e-16 *** typeB 1.7957 0.1666 10.777 < 2e-16 *** typeC -1.2528 0.3273 -3.827 0.000130 *** typeD -0.9045 0.2875 -3.146 0.001653 ** typeE -0.1178 0.2346 -0.502 0.615698 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 614.54 on 33 degrees of freedom Residual deviance: 170.71 on 29 degrees of freedom AIC: 278.58 Number of Fisher Scoring iterations: 6

5. Ora costruiamo un secondo modello in cui ipotizziamo che il numero di incidenti dipenda soltanto dal periodo di operatività dell’imbarcazione > mod_op <- glm(damage~operation, family=poisson(link="log"), data=barche) > summary(mod_op) Call: glm(formula = damage ~ operation, family = poisson(link = "log"), data = barche) DevianceResiduals: Min 1Q Median 3Q Max -4.77934 -3.99640 -2.47040 0.09608 10.73683 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 2.22642 0.08482 26.249 <2e-16 *** operation1975-79 0.20903 0.10864 1.924 0.0543 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersionparameterforpoisson family takentobe 1) Nulldeviance: 614.54 on 33 degreesoffreedom Residualdeviance: 610.79 on 32 degreesoffreedom AIC: 712.65

6. Costruiamo un altro modello in cui ipotizziamo che il numero di incidenti dipenda soltanto dall’anno di costruzione dell’imbarcazione > mod_con <- glm(damage~construction, family=poisson(link="log"), data=barche) > summary(mod_con) Call: glm(formula = damage ~ construction, family = poisson(link = "log"), data = barche) DevianceResiduals: Min 1Q Median 3Q Max -5.15752 -3.84193 -2.54327 0.05806 9.02390 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 2.0513 0.1195 17.163 < 2e-16 *** construction1965-69 0.5365 0.1477 3.633 0.00028 *** construction1970-74 0.4168 0.1509 2.763 0.00573 ** construction1975-79 -0.1054 0.2070 -0.509 0.61079 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersionparameterforpoisson family takentobe 1) Nulldeviance: 614.54 on 33 degreesoffreedom Residualdeviance: 592.51 on 30 degreesoffreedom AIC: 698.38 Numberof Fisher Scoringiterations: 6

7. Costruiamo infine un modello in cui ipotizziamo che il numero di incidenti dipenda soltanto dal numero di mesi di attività dell’imbarcazione: mod_mon <- glm(damage~months, family=poisson(link="log"), data=barche) > summary(mod_mon) Call: glm(formula = damage ~ months, family = poisson(link = "log"), data = barche) DevianceResiduals: Min 1Q Median 3Q Max -5.694 -3.251 -1.305 1.346 6.745 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 1.781e+00 7.207e-02 24.71 <2e-16 *** months 5.959e-05 2.921e-06 20.40 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersionparameterforpoisson family takentobe 1) Nulldeviance: 614.54 on 33 degreesoffreedom Residualdeviance: 309.15 on 32 degreesoffreedom AIC: 411.02 Numberof Fisher Scoringiterations: 5

8. Da questa prima serie di analisi l’unico regressore che non sembra essere molto significativo è operation. Il modello che finora si adatta meglio ai dati è mod_mon (AIC è il più basso) Proviamo ora a costruire un modello completo che contiene tutti e quattro I regressori e uno che li contiene tutti tranne operation e verifichiamo quale è il migliore. mod_tcm <- glm(damage~type+construction+months, family=poisson(link="log"), data=barche) > summary(mod_tcm) Call: glm(formula = damage ~ type + construction + months, family = poisson(link = "log"), data = barche) Deviance Residuals: Min 1Q Median 3Q Max -3.8174 -1.4399 -0.6102 0.7901 3.8942

9. Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 7.762e-01 2.422e-01 3.204 0.001353 ** typeB 9.492e-01 2.109e-01 4.500 6.81e-06 *** typeC -1.215e+00 3.274e-01 -3.711 0.000206 *** typeD -8.559e-01 2.876e-01 -2.976 0.002921 ** typeE -1.810e-01 2.347e-01 -0.771 0.440639 construction1965-69 1.011e+00 1.750e-01 5.778 7.54e-09 *** construction1970-74 1.354e+00 2.252e-01 6.012 1.83e-09 *** construction1975-79 9.279e-01 2.750e-01 3.374 0.000740 *** months 4.775e-05 7.385e-06 6.465 1.01e-10 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 614.54 on 33 degrees of freedom Residual deviance: 101.66 on 25 degrees of freedom AIC: 217.53 Number of Fisher Scoring iterations: 6 N.B: il coefficiente della variabile dummytypeE risulta non significativo: ciò indica che il numero di danni subiti dalle barche ti tipo E non è significativamente differente da quello delle barche di tipo A

10. > mod <- glm(damage ~type + construction+operation+months, family=poisson(link="log"), data=barche) > summary(mod) Call: glm(formula = damage ~ type + construction + operation + months, family = poisson(link = "log"), data = barche) Deviance Residuals: Min 1Q Median 3Q Max -2.5484 -1.3867 -0.4307 0.5222 3.1152 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 1.786e-01 2.771e-01 0.645 0.519201 typeB 6.701e-01 2.172e-01 3.085 0.002034 ** typeC -1.192e+00 3.275e-01 -3.638 0.000275 *** typeD -8.294e-01 2.877e-01 -2.883 0.003941 ** typeE -1.493e-01 2.349e-01 -0.636 0.524853 construction1965-69 1.087e+00 1.792e-01 6.067 1.30e-09 *** construction1970-74 1.500e+00 2.247e-01 6.673 2.50e-11 *** construction1975-79 8.545e-01 2.759e-01 3.097 0.001955 ** operation1975-79 7.284e-01 1.357e-01 5.369 7.93e-08 *** months 6.697e-05 8.523e-06 7.857 3.92e-15 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 614.539 on 33 degrees of freedom Residual deviance: 70.498 on 24 degrees of freedom AIC: 188.36 Number of Fisher Scoring iterations: 5

11. Considerando l’AIC il modello completo sembra essere il migliore. In quest’ultimo caso oltre al coefficiente di typeE anche il coefficiente dell’intercetta (β0) risulta non essere significativo. La conferma arriva anche dall’analisi stepwise: > step(mod, direction="backward") Start: AIC=188.36 damage ~ type + construction + operation + months DfDeviance AIC <none> 70.498 188.364 - operation 1 101.663 217.529 - type 4 122.926 232.793 - construction 3 135.854 247.721 - months 1 139.085 254.952 Call: glm(formula = damage ~ type + construction + operation + months, family = poisson(link = "log"), data = barche) Coefficients: (Intercept) typeBtypeC 1.786e-01 6.701e-01 -1.192e+00 typeDtypeE construction1965-69 -8.294e-01 -1.493e-01 1.087e+00 construction1970-74 construction1975-79 operation1975-79 1.500e+00 8.545e-01 7.284e-01 months 6.697e-05 DegreesofFreedom: 33 Total (i.e. Null); 24 Residual NullDeviance: 614.5 ResidualDeviance: 70.5 AIC: 188.4

12. Significato dei coefficienti: Covariate quantitative: Il coefficiente β9 relativo alla variabile quantitativa months indica che il logaritmo del valore atteso dei danni subiti da un imbarcazione aumenta di 6.697e-05 per ogni mese di attività in più. Di conseguenza il valore exp(β9) indica... Covariate qualitative (dicotomiche o politomiche): Il coefficiente β1 relativo alla variabile dummy typeB, invece, indica che le imbarcazioni di tipo B sono caratterizzate da un valore atteso condizionale di numero di danni subiti che è exp(0.67) volte quello delle barche di tipo A (prese come riferimento, nel caso in esame). Qual è il modello definitivo che possiamo stimare?

13. Problema della sovradispersione: La distribuzione di Poisson (e il relativo modello di regressione) presenta un difetto importante: implica la condizione di equidispersione. Spesso questa condizione non è supportata dai dati. Risulta infatti frequente il caso di sovradispersione in cui V(y|x) > E(y|x) = exp(x’β) In R la sovradispersione si riconosce nel caso in cui la residual deviance è molto più grande dei gradi di libertà. In questo caso si può ricorrere all’utilizzo della quasi-verosimiglianza: > mod_quasi <- glm(damage~type+construction+operation+months, family=quasipoisson(), data=barche) > summary(mod_quasi) Call: glm(formula = damage ~ type + construction + operation + months, family = quasipoisson(), data = barche) Deviance Residuals: Min 1Q Median 3Q Max -2.5484 -1.3867 -0.4307 0.5222 3.1152

14. Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.786e-01 4.587e-01 0.389 0.700447 typeB 6.701e-01 3.595e-01 1.864 0.074653 . typeC -1.192e+00 5.422e-01 -2.198 0.037863 * typeD -8.294e-01 4.763e-01 -1.741 0.094420 . typeE -1.493e-01 3.888e-01 -0.384 0.704284 construction1965-69 1.087e+00 2.967e-01 3.665 0.001224 ** construction1970-74 1.500e+00 3.721e-01 4.031 0.000487 *** construction1975-79 8.545e-01 4.568e-01 1.871 0.073628 . operation1975-79 7.284e-01 2.246e-01 3.243 0.003461 ** months 6.697e-05 1.411e-05 4.746 7.92e-05 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for quasipoisson family taken to be 2.740679) Null deviance: 614.539 on 33 degrees of freedom Residual deviance: 70.498 on 24 degrees of freedom AIC: NA Number of Fisher Scoring iterations: 5 Osserviamo che alcuni coefficienti di regressione non risultano più significativi e che le stime dei loro standard error sono ora molto più grandi.