Transcript Lab_R_6

LABORATORIO DI STATISTICA
AZIENDALE
6. Modello di Regressione di Poisson con R
Enrico Properzi - [email protected]
A.A. 2010/2011
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
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 ...
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
-4.6716 -2.2039 -0.5921
3Q
0.8632
Max
4.0113
Coefficients:
Estimate Std. Error z value
(Intercept)
1.7918
0.1543 11.612
typeB
1.7957
0.1666 10.777
typeC
-1.2528
0.3273 -3.827
typeD
-0.9045
0.2875 -3.146
typeE
-0.1178
0.2346 -0.502
--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01
Pr(>|z|)
< 2e-16
< 2e-16
0.000130
0.001653
0.615698
***
***
***
**
‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 614.54
Residual deviance: 170.71
AIC: 278.58
on 33
on 29
degrees of freedom
degrees of freedom
Number of Fisher Scoring iterations: 6
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)
Deviance Residuals:
Min
1Q
Median
-4.77934 -3.99640 -2.47040
3Q
0.09608
Max
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
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 614.54
Residual deviance: 610.79
AIC: 712.65
on 33
on 32
degrees of freedom
degrees of freedom
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)
Deviance Residuals:
Min
1Q
Median
-5.15752 -3.84193 -2.54327
3Q
0.05806
Max
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
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 614.54
Residual deviance: 592.51
AIC: 698.38
on 33
on 30
degrees of freedom
degrees of freedom
Number of Fisher Scoring iterations: 6
***
***
**
‘ ’ 1
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)
Deviance Residuals:
Min
1Q Median
-5.694 -3.251 -1.305
3Q
1.346
Max
6.745
Coefficients:
Estimate Std. Error z
(Intercept) 1.781e+00 7.207e-02
months
5.959e-05 2.921e-06
--Signif. codes: 0 ‘***’ 0.001 ‘**’
value Pr(>|z|)
24.71
<2e-16 ***
20.40
<2e-16 ***
0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 614.54
Residual deviance: 309.15
AIC: 411.02
on 33
on 32
degrees of freedom
degrees of freedom
Number of Fisher Scoring iterations: 5
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
-3.8174 -1.4399 -0.6102
3Q
0.7901
Max
3.8942
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
Residual deviance: 101.66
AIC: 217.53
on 33
on 25
degrees of freedom
degrees of freedom
Number of Fisher Scoring iterations: 6
N.B: il coefficiente della variabile dummy typeE 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
> 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
-2.5484 -1.3867 -0.4307
3Q
0.5222
Max
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 ‘
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 614.539
Residual deviance: 70.498
AIC: 188.36
on 33
on 24
degrees of freedom
degrees of freedom
Number of Fisher Scoring iterations: 5
**
***
**
***
***
**
***
***
’ 1
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
Df Deviance
<none>
70.498
- operation
1 101.663
- type
4 122.926
- construction 3 135.854
- months
1 139.085
AIC
188.364
217.529
232.793
247.721
254.952
Call: glm(formula = damage ~ type + construction + operation + months,
poisson(link = "log"), data = barche)
Coefficients:
(Intercept)
1.786e-01
typeD
-8.294e-01
construction1970-74
1.500e+00
months
6.697e-05
typeB
6.701e-01
typeE
-1.493e-01
construction1975-79
8.545e-01
typeC
-1.192e+00
construction1965-69
1.087e+00
operation1975-79
7.284e-01
Degrees of Freedom: 33 Total (i.e. Null); 24 Residual
Null Deviance:
614.5
Residual Deviance: 70.5
AIC: 188.4
family =
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?
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
-2.5484 -1.3867 -0.4307
3Q
0.5222
Max
3.1152
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
Residual deviance: 70.498
AIC: NA
on 33
on 24
degrees of freedom
degrees of freedom
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.