Autor opracowania: Marek Walesiak
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ć z
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
1
Autor opracowania: Marek Walesiak
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.
y
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)
maddala_data_pwn.csv
Plik maddala_data_pwn
W modelu logitowym prawdopodobieństwo przyjmowania przez zmienną yi jednej z dwóch
możliwych wartości jest funkcją xTb (i  numer obserwacji):
i
1 exp(xTb) exp(b0 +� b1T +� b2Y +� b3LF +� b4NW )
i
Pi =� F(xTb) =� =� =� ,
i
1+� exp(-�xTb) 1+� exp(xTb) 1+� exp(b0 +� b1T +� b2Y +� b3LF +� b4NW )
i i
gdzie: F  dystrybuanta rozkładu logistycznego.
Wartości funkcji odwrotnej do F dla tego modelu nazywa się logitami:
-�1
F (Pi ) =� ln(Pi 1-� P1) =� xTb =� b0 +� b1T +� b2Y +� b3LF +� b4NW
i
2
Autor opracowania: Marek Walesiak
2. Wykorzystując w programie R procedurę logit.r (dla modelu logitowego):
a) oszacować parametry modeli logitowego ln(Pi 1-� Pi ) =� b0 +� b1T +� b2Y +� b3LF +� b4NW . 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
ln(Pi 1-� Pi ) =� b0 +� b1T +� b2Y +� b3LF +� b4NW (interpretacja w kategoriach prawdopodobień-
stwa oraz interpretacja ilorazu szans Pi 1-� Pi ),
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 Pi dla obserwacji znajdujących się w próbie,
 prognozę prawdopodobieństwa Pi dla obserwacji znajdujących się poza próbą
 prognozę wartości yi (1 lub 0), tj. prognozę zmiennej objaśnianej dla i-tej obserwacji 
próba zbilansowana,
 prognozę wartości yi (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 ln(Pi 1-� Pi ) =� b0 +� b1T +� b2Y +� b3LF +� b4NW
[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 ln(Pi 1-� Pi ) =� b0 +� b1T +� b2Y +� b3LF +� b4NW z oszacowa-
nymi parametrami podając w nawiasach pod ocenami estymatorów parametrów ich błędy
ln(Pi 1-� Pi ) =�16,567+� 0,017T +� 9,132Y -� 0,715LF +� 85,362 NW
(19,639) (0,010) (5,053) (0,479) (35,854)
3
Autor opracowania: Marek Walesiak
b) podać interpretację parametrów modelu logitowego
ln(Pi 1-� Pi ) =� b0 +� b1T +� b2Y +� b3LF +� b4NW (interpretacja w kategoriach prawdopodobień-
stwa oraz interpretacja ilorazu szans Pi 1-� Pi ),
Interpretacja w kategoriach prawdopodobieństwa
Znak oszacowanego parametru przy zmiennej X (kolumna 3) określa kierunek wpływu zmien-
j
nej X na prawdopodobieństwo Pi :
j
 dla dodatniego bj wzrost (spadek) X wiąże się ze wzrostem (spadkiem) szans na to, że
j
yi =�1;
 dla ujemnego bj wzrost (spadek) X wiąże się ze spadkiem (wzrostem) szans na to, że
j
yi =�1.
Znak pa- Kolej-
Zmienna Parametr rametru Iloraz szans ność
T
b1=0,016516 + exp(b1) =�1,017 3
Y
b2 =9,131548 + exp(b2) =� 9242,321 2
LF b3 = 0,715389  4
exp(b3) =� 0,489
NW
b4 =85,3616 + exp(b4) =�1,181e +� 37 (przecinek 37 miejsc w prawo) 1
Interpretacja ilorazu szans Pi 1-� Pi
[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 xij (ceteris paribus) exp(bj ) pokazuje
krotność o jaką się zmienia iloraz szans (stosunek szansy (prawdopodobieństwa) na to, że yi =�1 do
szansy na to, że yi =� 0 ).
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
ln(Pi 1-� Pi ) =� b0 +� b1T +� b2Y +� b3LF +� b4NW
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 b0 a� =� 0,10 <� 0,3989 nie ma podstaw do odrzucenia hipotezy zerowej.
Oznacza to, że parametr b0 nieistotnie różni się od zera.
Z uwagi na to, że dla b1 a� =� 0,10 >� 0,0857 hipotezę zerową odrzucamy. Oznacza to, że parametr
b1 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 b2 a� =� 0,10 >� 0,0707 hipotezę zerową odrzucamy. Oznacza to, że parametr
b2 istotnie różni się od zera. Zmienna objaśniająca Y ma istotny wpływ na zmienną objaśnianą.
4
Autor opracowania: Marek Walesiak
Z uwagi na to, że dla b3 a� =� 0,10 <� 0,1355 nie ma podstaw do odrzucenia hipotezy zerowej. Za-
tem parametr b3 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 b4 a� =� 0,10 >� 0,0173 hipotezę zerową odrzucamy. Oznacza to, że parametr
b4 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
H0 : b1 =�K�=� bm =� 0 ,
H1 : przynajmniej jeden bj ą� 0 ( j =�1,K�,m ).
LUR
Statystyka testu ma postać: LR =� 2�� ln ,
LR
gdzie: LUR  maksimum funkcji wiarygodności, przy maksymalizacji względem wszystkich para-
metrów (dla pełnego modelu),
LR  maksimum funkcji wiarygodności przy maksymalizacji z warunkiem "�bj =� 0 (dla mo-
j
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 a� =� 0,10 >� 2,666e -� 05 (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 a� =� 0,10 )
[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ł [�-�15,737; 48,871]� pokryje nieznaną wartość parametru
b0 z modelu ln(Pi 1-� Pi ) =� b0 +� b1T +� b2Y +� b3LF +� b4NW .
Z prawdopodobieństwem 0,90 przedział [�0,0007; 0,032]� pokryje nieznaną wartość parametru b1
z modelu ln(Pi 1-� Pi ) =� b0 +� b1T +� b2Y +� b3LF +� b4NW .
Z prawdopodobieństwem 0,90 przedział [�0,820; 17,443]� pokryje nieznaną wartość parametru b2
z modelu ln(Pi 1-� Pi ) =� b0 +� b1T +� b2Y +� b3LF +� b4NW .
5
Autor opracowania: Marek Walesiak
Z prawdopodobieństwem 0,90 przedział [�-�1,504; 0,073]� pokryje nieznaną wartość parametru b3
z modelu ln(Pi 1-� Pi ) =� b0 +� b1T +� b2Y +� b3LF +� b4NW .
Z prawdopodobieństwem 0,90 przedział [�2,639; 144,337]� pokryje nieznaną wartość parametru
b4 z modelu ln(Pi 1-� Pi ) =� b0 +� b1T +� b2Y +� b3LF +� b4NW .
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
Wartość miary
Miara dopasowania R2
Kwadrat współczynnika korelacji między wartościami empirycznymi yi i
0,5982244
teoretycznymi w
i .
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 Pi 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
6
Autor opracowania: Marek Walesiak
f) wyznaczyć prognozy na podstawie modelu logitowego:
 prognozę prawdopodobieństwa Pi 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 ( Pi =� 0,2653).
f) wyznaczyć prognozy na podstawie modelu logitowego:
 prognozę wartości yi (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 yi (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
7