Autor opracowania: Marek Walesiak
1
PROJEKT D – MODELOWANIE I PROGNOZOWANIE
ZMIENNYCH JAKOŚCIOWYCH
Nazwisko i imię studenta 1: ..........................................
Kierunek i rok studiów studenta 1: ......
Numer grupy studenta 1: .....
Nazwisko i imię studenta 2: ..........................................
Kierunek i rok studiów studenta 2: ......
Numer grupy studenta 2: .....
Uwagi dla studentów:
1. Program R należy pobrać ze strony: http://cran.r-project.org/
2. Co najmniej jeden projekt (A, B, C, D) należy przesłać na e-mail prowadzącego laboratoria
3. Projekty można wykonywać osobiście lub w zespołach dwuosobowych (liczba zrealizowanych
projektów oraz jakość i estetyka wykonania będzie decydować o ocenie z laboratorium dla
przedmiotu Ekonometria)
4. Liczba obserwacji (dane w postaci szeregów przekrojowych z roku 2009 lub 2010) w projekcie
A, B oraz C musi wynosić co najmniej 12, a w projekcie D co najmniej 30. Dla danych staty-
stycznych należy koniecznie podać źródło
5. Nie wolno w projektach stosować zmiennych użytych w przykładowych projektach prezentowa-
nych na laboratoriach (nie dotyczy projektu C)
6. Wraz z każdym projektem opracowanym w edytorze Word (może też być jego odpowiednik z
pakietu OpenOffice) należy przesłać:
a) plik (pliki) danych w formacie csv
b) odpowiednie procedury w programie R
7. Termin przesłania projektu (projektów): do 03 stycznia 2012 roku
8. Proszę przesyłać projekty z własnych e-maili podając w e-mailu skład zespołu (imię i nazwisko,
rok i forma studiów, numer grupy lub specjalność)
9. Warunkiem przyjęcia projektu (projektów) jest uzyskanie pozytywnej odpowiedzi od prowadzą-
cego laboratoria
10. Odpowiedzi na e-maile informujące o akceptacji projektu lub projektów będą przesyłane w cią-
gu siedmiu dni od ich nadesłania
11. Odrzucane będą projekty, które wykonali inni studenci
Autor opracowania: Marek Walesiak
2
PROJEKT D – MODELOWANIE I PROGNOZOWANIE
ZMIENNYCH JAKOŚCIOWYCH
1. Dane przekrojowe dla stanów USA w 1950 r.
D1 – skłonność danego stanu do utrzymywania w stanowym kodeksie karnym kary śmierci (1 – dla stanów,
w których występuje kara śmierci, 0 – dla stanów bez kary śmierci),
T – mediana czasu spędzonego w więzieniu przez skazanych zabójców zwolnionych w 1951 r. (w miesią-
cach),
Y – mediana dochodu rodzin w 1949 r. (w tys. $),
LF – stopa zatrudnienia w 1950 r. (w %),
NW – udział osób nie-białych w populacji.
Źródło: Maddala [2006], s. 380.
a) wprowadzić dane statystyczne do programu EXCEL w następującym układzie
(plik maddala_data_pwn)
b) zapisać dane w formacie csv na dysku
(podać nazwę pliku maddala_data_pwn.csv)
Plik maddala_data_pwn
maddala_data_pwn.csv
W modelu logitowym prawdopodobieństwo przyjmowania przez zmienną
i
y jednej z dwóch
możliwych wartości jest funkcją
b
x
T
i
(i – numer obserwacji):
)
exp(
1
)
exp(
)
exp(
1
)
exp(
)
exp(
1
1
)
(
4
3
2
1
0
4
3
2
1
0
NW
b
LF
b
Y
b
T
b
b
NW
b
LF
b
Y
b
T
b
b
F
P
T
i
T
i
T
i
T
i
i
b
x
b
x
b
x
b
x
,
gdzie: F – dystrybuanta rozkładu logistycznego.
Wartości funkcji odwrotnej do F dla tego modelu nazywa się logitami:
NW
b
LF
b
Y
b
T
b
b
P
P
P
F
T
i
i
i
4
3
2
1
0
1
1
)
1
ln(
)
(
b
x
Autor opracowania: Marek Walesiak
3
2. Wykorzystując w programie R procedurę logit.r (dla modelu logitowego):
a) oszacować parametry modeli logitowego
NW
b
LF
b
Y
b
T
b
b
P
P
i
i
4
3
2
1
0
)
1
ln(
. Zapi-
sać postać modelu logitowego z oszacowanymi parametrami podając w nawiasach pod oce-
nami estymatorów parametrów ich błędy,
b) podać interpretację parametrów modelu logitowego
NW
b
LF
b
Y
b
T
b
b
P
P
i
i
4
3
2
1
0
)
1
ln(
(interpretacja w kategoriach prawdopodobień-
stwa oraz interpretacja ilorazu szans
i
i
P
P
1
),
c) za pomocą testu z sprawdzić istotność współczynników regresji, a za pomocą testu ilorazu
wiarygodności sprawdzić istotność całego modelu logitowego,
d) wyznaczyć i zinterpretować przedziały ufności dla parametrów strukturalnych modelu logito-
wego,
e) obliczyć i ocenić dopasowanie modelu logitowego,
f) wyznaczyć prognozy na podstawie modelu logitowego:
– prognozę prawdopodobieństwa
i
P dla obserwacji znajdujących się w próbie,
– prognozę prawdopodobieństwa
i
P dla obserwacji znajdujących się poza próbą
– prognozę wartości
i
y (1 lub 0), tj. prognozę zmiennej objaśnianej dla i-tej obserwacji –
próba zbilansowana,
– prognozę wartości
i
y (1 lub 0), tj. prognoza zmiennej objaśnianej dla i-tej obserwacji –
próba niezbilansowana
ODPOWIEDZI Z WYKORZYSTANIEM obliczeń w programie R
a) oszacować parametry modeli logitowego
NW
b
LF
b
Y
b
T
b
b
P
P
i
i
4
3
2
1
0
)
1
ln(
[1] Wyniki estymacji
Call:
glm(formula = D1 ~ T + Y + LF + NW, family = binomial(link = "logit"),
data = d, x = TRUE, y = TRUE)
Deviance Residuals:
Min 1Q Median 3Q Max
-2,072e+00 1,119e-05 2,934e-02 2,479e-01 1,822e+00
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 16,566989 19,639238 0,844 0,3989
T 0,016516 0,009611 1,718 0,0857 .
Y 9,131548 5,053148 1,807 0,0707 .
LF -0,715389 0,479267 -1,493 0,1355
NW 85,361600 35,854424 2,381 0,0173 *
---
Signif. codes: 0 ‘***’ 0,001 ‘**’ 0,01 ‘*’ 0,05 ‘.’ 0,1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 44,584 on 43 degrees of freedom
Residual deviance: 18,215 on 39 degrees of freedom
AIC: 28,215
Number of Fisher Scoring iterations: 9
a) zapisać postać modelu logitowego
NW
b
LF
b
Y
b
T
b
b
P
P
i
i
4
3
2
1
0
)
1
ln(
z oszacowa-
nymi parametrami podając w nawiasach pod ocenami estymatorów parametrów ich błędy
NW
LF
Y
T
P
P
i
i
)
854
,
35
(
)
479
,
0
(
)
053
,
5
(
)
010
,
0
(
)
639
,
19
(
362
,
85
715
,
0
132
,
9
017
,
0
567
,
16
)
1
ln(
Autor opracowania: Marek Walesiak
4
b) podać interpretację parametrów modelu logitowego
NW
b
LF
b
Y
b
T
b
b
P
P
i
i
4
3
2
1
0
)
1
ln(
(interpretacja w kategoriach prawdopodobień-
stwa oraz interpretacja ilorazu szans
i
i
P
P
1
),
Interpretacja w kategoriach prawdopodobieństwa
Znak oszacowanego parametru przy zmiennej
j
X
(kolumna 3) określa kierunek wpływu zmien-
nej
j
X
na prawdopodobieństwo
i
P :
– dla dodatniego
j
b
wzrost (spadek)
j
X
wiąże się ze wzrostem (spadkiem) szans na to, że
1
i
y
;
– dla ujemnego
j
b
wzrost (spadek)
j
X
wiąże się ze spadkiem (wzrostem) szans na to, że
1
i
y
.
Zmienna
Parametr
Znak pa-
rametru
Iloraz szans
Kolej-
ność
T
1
b
=0,016516
+
017
,
1
)
exp(
1
b
3
Y
2
b
=9,131548
+
321
,
9242
)
exp(
2
b
2
LF
3
b
=–0,715389
–
489
,
0
)
exp(
3
b
4
NW
4
b
=85,3616
+
37
181
,
1
)
exp(
4
e
b
(przecinek 37 miejsc w prawo)
1
Interpretacja ilorazu szans
i
i
P
P
1
[1] Interpretacja parametrów modelu logitowego - iloraz szans Pi(1-Pi)
(Intercept) T Y LF NW
1,566577e+07 1,016654e+00 9,242321e+03 4,890017e-01 1,180516e+37
Dla jednostkowego przyrostu zmiennej objaśniającej
ij
x (ceteris paribus)
)
exp(
j
b pokazuje
krotność o jaką się zmienia iloraz szans (stosunek szansy (prawdopodobieństwa) na to, że
1
i
y
do
szansy na to, że
0
i
y
).
Największy iloraz szans występuje dla parametru przy zmiennej NW, najmniejszy zaś przy
zmiennej LF.
c) za pomocą testu z sprawdzić istotność współczynników regresji dla modelu logitowego
NW
b
LF
b
Y
b
T
b
b
P
P
i
i
4
3
2
1
0
)
1
ln(
Test z
Estimate Std. Error z value Pr(>|z|)
(Intercept) 16,566989 19,639238 0,844 0,3989
T 0,016516 0,009611 1,718 0,0857 .
Y 9,131548 5,053148 1,807 0,0707 .
LF -0,715389 0,479267 -1,493 0,1355
NW 85,361600 35,854424 2,381 0,0173 *
Signif. codes: 0 ‘***’ 0,001 ‘**’ 0,01 ‘*’ 0,05 ‘.’ 0,1 ‘ ’ 1
Z uwagi na to, że dla
0
b
3989
,
0
10
,
0
nie ma podstaw do odrzucenia hipotezy zerowej.
Oznacza to, że parametr
0
b nieistotnie różni się od zera.
Z uwagi na to, że dla
1
b
0857
,
0
10
,
0
hipotezę zerową odrzucamy. Oznacza to, że parametr
1
b istotnie różni się od zera. Zmienna objaśniająca T ma istotny wpływ na zmienną objaśnianą.
Z uwagi na to, że dla
2
b
0707
,
0
10
,
0
hipotezę zerową odrzucamy. Oznacza to, że parametr
2
b istotnie różni się od zera. Zmienna objaśniająca Y ma istotny wpływ na zmienną objaśnianą.
Autor opracowania: Marek Walesiak
5
Z uwagi na to, że dla
3
b
1355
,
0
10
,
0
nie ma podstaw do odrzucenia hipotezy zerowej. Za-
tem parametr
3
b nieistotnie różni się od zera i zmienna objaśniająca LF ma nieistotny wpływ na
zmienną objaśnianą.
Z uwagi na to, że dla
4
b
0173
,
0
10
,
0
hipotezę zerową odrzucamy. Oznacza to, że parametr
4
b istotnie różni się od zera. Zmienna objaśniająca NW ma istotny wpływ na zmienną objaśnianą.
c) za pomocą testu ilorazu wiarygodności sprawdzić istotność całego modelu logitowego
:
0
H
0
1
m
b
b
,
:
1
H
przynajmniej jeden
0
j
b
(
m
j
,
,
1
).
Statystyka testu ma postać:
R
UR
L
L
LR
ln
2
,
gdzie:
UR
L
– maksimum funkcji wiarygodności, przy maksymalizacji względem wszystkich para-
metrów (dla pełnego modelu),
R
L – maksimum funkcji wiarygodności przy maksymalizacji z warunkiem
0
j
j
b
(dla mo-
delu tylko z wyrazem wolnym).
Statystyka
LR
ma rozkład chi-kwadrat z liczbą stopni swobody równą liczbie zmiennych obja-
śniających modelu pełnego.
[1] Test ilorazu wiarygodności dla badania istotności całego modelu
(LR=2*Lur/Lr)
Likelihood ratio test
Model 1: D1 ~ T + Y + LF + NW
Model 2: D1 ~ 1
#Df LogLik Df Chisq Pr(>Chisq)
1 5 -9,1076
2 1 -22,2921 -4 26,369 2,666e-05 ***
---
Signif. codes: 0 ‘***’ 0,001 ‘**’ 0,01 ‘*’ 0,05 ‘.’ 0,1 ‘ ’ 1
Z uwagi na to, że dla
e
,
05
666
2
10
,
0
(0,00002666) hipotezę zerową odrzucamy. Ozna-
cza to, że regresja jako całość jest istotna.
d) wyznaczyć i zinterpretować przedziały ufności dla parametrów strukturalnych dla mo-
delu logitowego (dla
10
,
0
)
[1] Przedziały ufności dla parametrów
p_uf_d p_uf_g
(Intercept) -1,573668e+01 48,87066062
T 7,075145e-04 0,03232541
Y 8,198601e-01 17,44323651
LF -1,503714e+00 0,07293508
NW 2,638632e+01 144,33687905
Z prawdopodobieństwem 0,90 przedział
871
,
48
737
,
15
;
pokryje nieznaną wartość parametru
0
b z modelu
NW
b
LF
b
Y
b
T
b
b
P
P
i
i
4
3
2
1
0
)
1
ln(
.
Z prawdopodobieństwem 0,90 przedział
032
,
0
0007
,
0
;
pokryje nieznaną wartość parametru
1
b
z modelu
NW
b
LF
b
Y
b
T
b
b
P
P
i
i
4
3
2
1
0
)
1
ln(
.
Z prawdopodobieństwem 0,90 przedział
443
,
17
;
820
,
0
pokryje nieznaną wartość parametru
2
b
z modelu
NW
b
LF
b
Y
b
T
b
b
P
P
i
i
4
3
2
1
0
)
1
ln(
.
Autor opracowania: Marek Walesiak
6
Z prawdopodobieństwem 0,90 przedział
073
,
0
;
504
,
1
pokryje nieznaną wartość parametru
3
b
z modelu
NW
b
LF
b
Y
b
T
b
b
P
P
i
i
4
3
2
1
0
)
1
ln(
.
Z prawdopodobieństwem 0,90 przedział
337
,
144
;
639
,
2
pokryje nieznaną wartość parametru
4
b z modelu
NW
b
LF
b
Y
b
T
b
b
P
P
i
i
4
3
2
1
0
)
1
ln(
.
Węższe (szersze) przedziały ufności można uzyskać poprzez zmniejszenie (zwiększenie) pozio-
mu ufności.
e) obliczyć i ocenić dopasowanie modelu logitowego
[1] Miary dopasowania
[1] Kwadrat współczynnika korelacji między wartościami empirycznymi i
teoretycznymi
[1] 0,5982244
[1] R kwadrat Efrona
[1] 0,5981534
[1] R kwadrat Nagelkerke
[1] 0,7077272
[1] R kwadrat McFaddena
[1] 0,5914449
Miara dopasowania
2
R
Wartość miary
Kwadrat współczynnika korelacji między wartościami empirycznymi
i
y i
teoretycznymi
i
yˆ .
0,5982244
Miara Efrona
0,5981534
Miara Nagelkerke
0,7077272
Miara McFaddena
0,5914449
Miara Nagelkerke pokazuje na dobre dopasowanie modelu logitowego do danych empirycznych.
Pozostałe miary wskazują, że model jest dość dobrze dopasowany do danych.
f) wyznaczyć prognozy na podstawie modelu logitowego:
– prognozę prawdopodobieństwa
i
P dla obserwacji znajdujących się w próbie
[1] 1a. Prognoza prawdopodobieństwa Pi dla obserwacji znajdujących się w próbie
1 2 3 4 5 6 7
0,99999999 0,99996803 0,99970456 0,98990940 0,62715799 0,83958907 0,99998215
8 9 10 11 12 13 14
0,99999162 0,72803895 0,08571948 0,99973662 1,00000000 0,93407318 0,85886041
15 16 17 18 19 20 21
1,00000000 0,82231552 0,99998469 0,16736354 0,08708977 0,98157700 1,00000000
22 23 24 25 26 27 28
0,37526162 0,99999798 0,51856108 0,88313819 0,99862342 0,97188050 0,96284709
29 30 31 32 33 34 35
0,99940984 0,99839702 0,98885441 0,19028138 0,99710860 0,11119749 1,00000000
36 37 38 39 40 41 42
0,85525045 0,99853555 0,98312092 0,79876223 0,99999693 0,92817998 0,16669814
43 44
0,99836037 0,15447486
Autor opracowania: Marek Walesiak
7
f) wyznaczyć prognozy na podstawie modelu logitowego:
– prognozę prawdopodobieństwa
i
P dla obserwacji znajdujących się poza próbą
Przyjęto założenie, że dla pewnego stanu USA obserwacje na zmiennych objaśniających będą w
okresie prognozowanym następujące (zob. plik maddala_data_pwn_p.csv):
[1] 1b. Prognoza prawdopodobieństwa Pi dla obserwacji znajdujących się poza pró-
bą
45
0,2652799
Zatem dla tego stanu prognozuje się, że będzie niewielka skłonność do umieszczenia kary śmier-
ci w stanowym kodeksie karnym (
2653
,
0
ˆ
i
P
).
f) wyznaczyć prognozy na podstawie modelu logitowego:
– prognozę wartości
i
y (1 lub 0), tj. prognozę zmiennej objaśnianej dla i-tej obserwacji – próba
zbilansowana,
[1] 2a. Prognoza wartości yi (1 lub 0), tj. prognoza zmiennej objaśnianej
dla i-tej obserwacji – próba zbilansowana
[1] Tablica trafności
przewidywane
faktyczne 0 1
0 7 2
1 1 34
[1] Mierniki dokładności prognoz - próba zbilansowana
[1] Procentowa trafność prognozowania
[1] 93,18182
[1] Iloraz szans
[1] 119
f) wyznaczyć prognozy na podstawie modelu logitowego:
– prognozę wartości
i
y (1 lub 0), tj. prognoza zmiennej objaśnianej dla i-tej obserwacji – próba
niezbilansowana (w próbie obejmującej 44 stany dla zmiennej
D1
jest 35 „1” i 9 „0”. Jest to więc
próba niezbilansowana)
[1] 2b. Prognoza wartości yi (1 lub 0), tj. prognoza zmiennej objaśnianej
dla i-tej obserwacji – próba niezbilansowana
[1] Tablica trafności - próba niezbilansowana
przewidywane_n
faktyczne 0 1
0 7 2
1 4 31
[1] Mierniki dokładności prognoz - próba niezbilansowana
[1] Procentowa trafność prognozowania
[1] 86,36364
[1] Iloraz szans
[1] 27,125