07-wazenie.Rmd

# Estymacja w badaniach internetowych -- ważenie danych

## Wstęp

W przypadku badań interentowych możemy wyszczególnić następujące podejścia:

+ oparte na podejściu (pseudo)-randomizacyjnym -- przeważanie danych, w tym modelujemy skłonność do odpowiedzi (ang. propensity to respond),
+ oparte na podejściu modelowym -- budujemy model, którym chcemy wyjaśnić naszą cechę szeregiem cech $\x$.
+ mieszane -- budujemy model, a następnie przeważamy wyniki do znanych wartości globalnych


Niech $U=\{1, 2, ..., N\}$ oznacza populację generalną, a $Y$ oznacza cechę, którą badamy (np. liczba osób pracujących w pełnym wymiarze czasu). Niech $y_i$ oznacza wartości cechy obserwowane dla każdej jednostki populacji. Założmy, że celem jest oszacowanie wartości średniej w populacji

$$
\bar{Y} = \sum_{i=1}^{N}y_k.
$$

Jednak najczęściej nie obserwujemy całej populacji, a pewien jej fragment (próbę). Jeżeli

## Post-stratyfikacja 

Celem post-stratyfikacji jest:

+ sprowadzenie rozkładów próby do znanych rozkładów populacji,
+ zmniejszenie obciażenia estymatorów na podstawie próby przez skorygowanie rozkładów próby do znanych rozkładów z populacji.

Aby dokonać post-stratyfikacji potrzebujemy informacji o wartościach globalnych ($\bX$) zmiennych pomocniczych $\bx$, którymi mogą być płeć, wiek, województwo itp.  

Założmy, że zmienna $\bx$ ma $L$ kategorii (warstw; np. wiek ma 5 kategorii, płeć ma dwie kategorie), które dzielą populację na podpopulacje $U_1, U_2,..., U_L$, a ich wielkości możemy oznaczyć jako $N_h$, gdzie $h=1,2,...,L$, a $N = N_1 + N_2 + ... + N_L$. Zakłdamy, że wielkości warstw są znane (pytanie: **skąd mogą być znane?**). Założmy, że obserwujemy próbę o liczebności $n$, którą również możemy podzielić na podpróby $n_1, n_2, ..., n_h$, które dają $n=n_1 + n_2 + ... + n_h$.

Aby próba była reprezentatywna ze względu na cechę $\bx$ proporcja elementów w warstwie $h$ musi być równa $N_h/N$ dla $h=1,...,L$. W przypadku próby udział poszczególnych kategorii wyznaczamy jako $n_h / n$. W związku z tym aby dokonać korekty musimy dla każdej jednostki $i$ w warstwie $h$ dokonać następującej korekty

$$
c_i = \frac{N_h/N}{n_h/n} = \frac{N_h}{N} \frac{n}{n_h}.
$$

Jeżeli wagi wynikające z losowania są równe $d_i = n / N$, to estymator post-stratyfikowany dla średniej ma następującą postać:

$$
\bar{y}_{ps} = \frac{1}{N}\sum_{h=1}^{L} N_h\bar{y}^{(h)},
(\#eq:y-ps)
$$

gdzie $\bar{y}^{(h)}$ jest średnią obliczoną na podstawie elementów obserwowanych w warstwie $h$. W tym przypadku post-stratyfikowany estymator średniej jest równy ważonej sumie średnich w poszczególnych warstwach.


Teraz założmy, że nie wszystkie jednostki z próby udzieliły odpowiedzi (występuje problem non-response). W takim przypadku estymator  \@ref(eq:y-ps) na następującą postać:

$$
\bar{y}_{R,ps} = \frac{1}{N}\sum_{h=1}^{L} N_h\bar{y}^{(h)_R},
(\#eq:y-ps-r)
$$

gdzie $\bar{y}^{(h)_R}$ to średnia dla jednostek w warstwie $h$, które odpowiedziały. Obciążenie estymatora \@ref(eq:y-ps-r) dane jest wzorem:

$$
B(\bar{y}_{R,ps}) = \frac{1}{N} \sum_{h=1}^{L} N_h \frac{R_{\rho{y}}^{(h)} S_{\rho}^{(h)} S_{y}^{(h)} }{ \bar{\rho}^{(h)}},
$$

gdzie $R_{\rho{y}}^{(h)}$ jest współczynnikiem korelacji miedzy $y$, a $\rho$ w warstwie $h$, $S_{\rho}^{(h)},S_{y}^{(h)}$ jest odchyleniem standardowym $\rho$ oraz $y$ w warstwie $h$, a $\bar{\rho}^{(h)}$ jest średnią $\rho$ w warstwie $h$.

Obciązenie będzie małe w przypadku gdy:


+ Nie ma zależności między badaną cechą, a skłonnością do odpowiedzi w danej warstwie -- korelacja $R_{\rho{y}}^{(h)}$ jest niska,
+ Skłonności do odpowiedzi w warstwie są mniej więcej bliskie -- $S_{\rho}^{(h)}$ jest małe
+ Cecha $y$ jest słabo zróżnicowana w warstwach -- $S_{y}^{(h)}$ jest niskie.


Następnie musimy rozważyć przypadek gdy występuje błąd pokrycia (undercoverage). Wtedy obciązenie post-stratyfikowanego estymatora jest równe

$$
B(\bar{y}_{I,ps}) = \sum_{h=1}^{L} W_h \frac{N_{NI,h}}{N_h} (\bar{y}_{I}^{(h)} - \bar{y}_{NI}^{(h)}),
$$

gdzie $W_h = N_h / N$, $N_{NI,h}$ to liczba osób bez internetu w warstwie $h$, $\bar{y}_{I}^{(h)}$ to średnia dla osób, które mają dostęp do internetu, a $\bar{y}_{NI}^{(h)}$ to średnia dla osób, które nie mają dostępu do Internetu.

Obciażenie estymatora \@ref(eq:y-ps) w przypadku gdy wystąpiła autoselekcja jednostek do próby jest następujące: 

$$
B(\bar{y}_{S,ps}) = \frac{1}{N} \sum_{h=1}^{L} N_h \frac{R_{\rho{y}}^{(h)} S_{\rho}^{(h)} S_{y}^{(h)} }{ \bar{\rho}^{(h)}}.
$$


### Post-stratyfikacja w R

Załóżmy, że naszym celem jest oszacowanie liczby osób pracujących w 2015 roku (zmienna `czy_praca`) będących w wieku produkcyjnym. Aby oszacować te wartości wykorzystamy dwie warstwy płeć (`plec`) oraz województwo (`woj`). W pierwszej kolejności oszacujemy liczbę pracujących na podstawie próby bez korekty.

```{r}
proba_pracujacy <- diagnoza2015 %>%
  filter(wiek <= 59) %>%
  mutate(praca = as.numeric(czy_praca == 1)) %>%
  summarize(m = mean(praca))
proba_pracujacy
```

Liczba pracujących wyniosła `r round(proba_pracujacy[[1]]*100,2)`% w przypadku gdy oszacowanie przedstawimy na podstawie próby bez żadnej korekty. Następnie dokonamy korekty z wykorzystaniem informacji o wielkości populacji w poszczególnych warstwach. W pierwszej kolejności oszacujemy średnie według płci oraz województw.

```{r}
## 1 – mężczyzna, 2 – kobieta
diagnoza_pracujacy <- diagnoza2015 %>%
  filter(wiek <= 59) %>%
  mutate(praca = as.numeric(czy_praca == 1)) %>%
  group_by(plec, woj) %>%
  summarize(praca = mean(praca))

diagnoza_pracujacy
```

W kolejnym kroku potrzebujemy informacji o liczbe osób w wieku produkcyjnym w tych przekrojach. Wykorzystamy w tym celu dane z rejestru PESEL skorygowanego przez informacje ze spisu powszechnego NSP 2011 (dane pobrane z BDL).

```{r}
populacja_N <- populacja %>%
  filter(Płeć != 'ogółem') %>%
  mutate(plec = ifelse(Płeć == 'mężczyźni',1,2),
         woj = gsub('^\\d','',Kod),
         woj = substr(woj,1,2),
         woj = as.numeric(woj)) %>%
  count(woj, plec, wt = Wartosc)

populacja_N
```

Następnie dołączymy te informacje do próby i dokonamy szacunku

```{r}
proba_pracujacy_korekta <- diagnoza_pracujacy %>%
  left_join(populacja_N) %>%
  ungroup() %>% 
  summarize(prac = weighted.mean(x = praca, w = n))
```

Liczba pracujących po korekcie wyniosła `r round(proba_pracujacy_korekta[[1]]*100,2)`%, a przed `r round(proba_pracujacy[[1]]*100,2)`% podczas gdy wartość w populacji wynosiła 68.7% (+/ 0.4 pp). Spójrzmy co będzie jeżeli użyjemy płci i wieku.

```{r}
## w latach według grup wieku
populacja_N <- populacja %>%
  filter(Płeć != 'ogółem') %>%
  mutate(plec = ifelse(Płeć == 'mężczyźni',1,2)) %>%
  count(Wiek, plec, wt = Wartosc) %>%
  ungroup() %>%
  mutate(Wiek = factor(Wiek),
         wiek = as.numeric(Wiek)) %>%
  select(wiek,plec,n)

## wiek
diagnoza_pracujacy <- diagnoza2015 %>%
  filter(wiek <= 59) %>% 
  mutate(praca = as.numeric(czy_praca == 1),
         wiek_cat = cut(x = wiek, breaks = c(15,20,25,30,35,40,45,50,55,60),include.lowest = T, right = F)) %>%
  group_by(wiek_cat, plec) %>%
  summarize(praca = mean(praca)) %>%
  ungroup() %>%
  mutate(wiek = as.numeric(wiek_cat)) %>%
  select(wiek, plec, praca)

## łączymy i szacujemy
prac_wiek <- diagnoza_pracujacy %>%
  left_join(populacja_N) %>%
  ungroup() %>% 
  summarize(prac = weighted.mean(x = praca, w = n))
```

Podsumowując nasze oszacowania (wartość prawdziwa wynosiła 68.7%):

+ bez post-stratyfikacji -- `r round(proba_pracujacy[[1]]*100,2)` % (obciązenie= `r round(proba_pracujacy[[1]]*100,2) - 68.7`) 
+ post-stratyfikacja (płeć x województwo) -- `r round(proba_pracujacy_korekta[[1]]*100,2)`% (obciązenie= `r round(proba_pracujacy_korekta[[1]]*100,2) - 68.7`) 
+ post-stratyfikacja (płeć x grupy wieku) -- `r round(prac_wiek[[1]]*100,2)` % (obciązenie= `r round(prac_wiek[[1]]*100,2) - 68.7`) 

## Kalibracja


## Propensity score weighting

### Poszukiwanie zmiennych diagnostycznych

W związku z tym, że naszym celem jest oszacowanie modelu skłonności do odpowiedzi (ang. *response propensity*) gdzie zmienna $R = \{ 0, 1 \}$


$$
P(R_i = 1 | \boldsymbol{X}_i, I_i = 1),
$$

gdzie $I_i$ to zmienna określająca czy dana osoba korzysta z Internetu, a $\boldsymbol{X}_i$ to zestaw cech opisujących daną jednostkę. Kolejnym etapem jest sprawdzenie, które zmienne związane są ze skłonnością do odpowiedzi. Możemy to zweryfikować biorąc pod uwagę dwie grupy

1. Zmienne objaśniane mają charater jakościowy (np. płeć, grupy wieku, województwo)
    + test t-studenta (test parametryczny, porównanie dwóch średnich), funkcja `t.test()`
    + test Manna-Whitneya (test nieparametryczny, porównanie dwóch rozkładów), funkcja `wilcox.test()`


2. Zmienne objaśniane mają charakter ilościowy (np. wiek, dochody, czas odpowiedzi)
    + Test proporcji (ma sens w tabeli $2 \times 2$), funkcja `prop.test()` -- zastosowanie w przypadku dwóch zmiennych
    + (Dokładny) test Fishera (zwykle dla tabel $2 \times 2$ i małych prób), funkcja `fisher.test()` 
    + Ilorazy szans, funkcja `vcd::oddsratio()` oraz  test Woolf'a dla tabel $2 \times 2 \times K$, funkcja `vcd::woolf_test()` 
    + Testy oparte na statystyce $\chi^2$ -- więcej o tej grupie poniżej, m.in. funkcja `chisq.test()`
    + Test Cochrana-Mantela-Haenszel'a -- tylko jeżeli dwie zmienne są na skali porządkowej, funkcja `mantelhaen.test`
    + Uogólniony test Cochrana-Mantela-Haenszel'a -- stosujemy dla 3 zmiennych: zakładamy, że dwie zmienne są warunkowo niezależne, w warstwach utworzonych przez trzecią zmienną (zmienne mogą mieć charakter zarówno nominalny, jak i porządkowy) -- funkcja `mantelhaen.test`, `CMHtest() `z pakietu `vcdExtra`. Warto też zapoznać się z pakietem `coin`, który zawiera funkcje `cmh_test`.

Miary zależności:

+ współczynnik Yule’a -- dwie zmienne, mające porządek (tak, nie) i traktujemy symetryczne (zmiana porządku w tabeli nie wpływa na statystykę)
+ współczynnik C Pearsona -- dwie zmienne, nie mają porządku, traktujemy symetrycznie (odnosi się do tej i poniższych dwóch)
+ współczynnik V Cramera,
+ współczynnik T Czupurowa,
+ Goodman-Kruskal $\gamma$ -- dla zmiennych na skali porządkowej


#### Opis procedur testowych

1. Test $\chi^2$

```{r}
tab1 <- xtabs(~internet + plec, data = diagnoza2015)
chisq.test(tab1)
```

Wynik testu: występuje zależność miedzy korzystaniem z internetu, a płcią.

```{r}
tab2 <- xtabs(~internet + klm, data = diagnoza2015)
chisq.test(tab2)
```

Wynik testu: występuje zależność miedzy korzystaniem z internetu, a klasą wielkości miejscowości. Dodatkowo należy zauważyć, że statystyka jest ponad 5 razy większa od tej dla płci. 

```{r}
tab3 <- xtabs(~internet + woj, data = diagnoza2015)
chisq.test(tab3)
```

Wynik testu: występuje zależność miedzy korzystaniem z internetu, a województwem. 

2. Test t-Studenta (porównanie dwóch średnich)

```{r}
doubledecker(internet ~ wiek, diagnoza2015)
boxplot(wiek ~ internet, diagnoza2015)
var.test(wiek ~ internet, diagnoza2015)
t.test(wiek ~ internet, diagnoza2015)
wilcox.test(wiek ~ internet, diagnoza2015)
```

Porównanie median

```{r}
by(diagnoza2015$wiek,diagnoza2015$internet, median)
```

3.  Mantela-Haenszel'a $\chi^2$

Wykorzystujemy tylko wtedy jeżeli dane są na skali porządkowej! Hipoteza alternatywna jest liniowa zależność między

$$
Q_{MH} = (n-1) r^2 \sim \chi^2(1),
$$

gdzie $n$ to wielkość próby, a $r$ to współczynnik korelacji liniowej Pearsona. Statystykę tę możemy wykorzystać do weryfikacji następującego układu hipotez

+ $H_0$: brak związku między zmiennymi,
+ $H_1$: występuje związek między zmiennymi.

Możemy ją zbadać wykorzystując funkcję `vcdExtra::CMHtest()`

```{r}
doubledecker(internet~woj+plec, diagnoza2015)
tab3 <- xtabs(~internet+plec+woj,diagnoza2015)
mantelhaen.test(tab3)
woolf_test(tab3) ## czy można stosować test CMH ?
```


3. Uogólniony test Cochrana-Mantela-Haenszel'a

Zakładamy tabelę $2\times2\times K$ 

$$
Q_{CMH}=\frac{\left[ \sum_{k} (n_{11k}-\mu_{11k}) \right]^2}{\sum_{k}var(n_{11k})}
$$

gdzie 

$$
\mu_{11k} = E(n_{11k}) = n_{1+k}n_{+1k}/n_{++k},
$$

$$
var(n_{11k}) = \frac{n_{1+k}n_{2+k}n_{+1k}n_{+2k}}{n^2_{++k}(n_{++k}-1)}
$$

Obliczamy ja korzystając z funkcji `vcdExtra::CMHtest`

```{r}
diagnoza2015$wiek_4g <- cut(diagnoza2015$wiek,quantile(diagnoza2015$wiek))
tab4 <- xtabs(~internet+wiek_4g+woj,diagnoza2015)
mantelhaen.test(tab4)
```


### Regresja logistyczna

Funkcją logistyczną nazywamy funkcję postaci:

$$
f\left(x\right)=\frac{1}{1+e^{-x}}=\frac{e^{x}}{1+e^{x}}
$$


Funkcja logistyczna -- własności: 

+ $D=R$.
+ Funkcja logistyczna przyjmuje wartości z~przedziału $\left(0,1\right)$.
+ Funkcja logistyczna kształtem przypomina rozciągniętą literę S.
+ Funkcja logistyczna nadaje się do opisu takich zjawisk, które początkowo rosną coraz szybciej, następnie osiągają punkt przegięcia i~rosną coraz wolniej, zbliżając się asymptotycznie do poziomu nasycenia.
+ Popularność funkcji logistycznej wynika z~faktu, że przyjmuje ona wartości z~przedziału $\left(0,1\right)$. Może być ona zatem wykorzystana do modelowania prawdopodobieństwa.
+ Prawdopodobieństwo to może, na przykład, oznaczać ryzyko zachorowania lub szansę na wyzdrowienie, ryzyko niespłacenia kredytu bądź szansę jego spłaty, szansę zakupienia pewnego produktu bądź ryzyko zakupienia innego towaru etc.	    

Pojęcie szansy

+ Szansą nazywamy stosunek prawdopodobieństwa, że pewne zdarzenie wystąpi do prawdopodobieństwa, że to zdarzenie nie pojawi się.
+ Niech $A$ oznacza dowolne zdarzenie. Z~formalnego punktu widzenia szansę definiujemy w~następujący sposób:

$$
S\left(A\right)=\frac{P\left(A\right)}{1-P\left(A\right)}.
$$

Załóżmy, że prawdopodobieństwo zdarzenia, że klient supermarketu zakupi pewien produkt wynosi 0.1. Wówczas szansa zakupu tego produktu przez klienta supermarketu wynosi:

$$
S\left(A\right)=\frac{P\left(A\right)}{1-P\left(A\right)}=\frac{0.1}{1-0.1}=\frac{0.1}{0.9}=\frac{1}{9}.
$$

Oznacza to, że prawdopodobieństwo zakupu produktu jest równe jednej dziewiątej prawdopodobieństwa, że klient nie zakupi tego produktu. Innymi słowy możemy powiedzieć, że szansa zakupu danego produktu jest jak 1 do 9, tzn. na 10 klientów jeden z~nich zakupi dany towar.

Pojęcie ilorazu szans

+ Ilorazem szans dwóch zdarzeń $A$ i $B$ nazywamy stosunek szansy wystąpienia zdarzenia $A$ do szansy wystąpienia zdarzenia $B$.
+ Niech $A$ i $B$ oznaczają dowolne zdarzenia. Z~formalnego punktu widzenia iloraz szans definiujemy w~następujący sposób:
$$
OR_{AB}=\frac{S\left(A\right)}{S\left(B\right)}=\frac{P\left(A\right)}{1-P\left(A\right)}:\frac{P\left(B\right)}{1-P\left(B\right)}.
$$
+ Oznaczenie ilorazu szansy $OR_{AB}$ pochodzi od początkowych liter nazwy angielskiej (\textit{odds ratio}).
			
			
Pojęcie ilorazu szans -- przykład

Załóżmy, że poniższa tabela przedstawia informację na temat wyników badania dotyczącego kupna pewnego produktu. Niech $A$ oznacza zdarzenie polegające na kupnie tego produktu przez kobietę, a~$B$ zdarzenie polegające na kupnie tego produktu przez mężczyznę.

\begin{table}[htp]
\begin{tabular}{|l|l|l|}
\hline
\textbf{Płeć} & \multicolumn{2}{c|}{\textbf{Informacja o kupnie produktu}} \\ 
\cline{2-3}
\textbf{klienta} & \textbf{Produkt został kupiony} & \textbf{Produkt nie został kupiony} \\ 
\hline
\textbf{Kobieta} & \multicolumn{1}{c|}{243} & \multicolumn{1}{c|}{30} \\ 
\hline
\textbf{Mężczyzna} & \multicolumn{1}{c|}{48} & \multicolumn{1}{c|}{240} \\ 
\hline
\end{tabular}
\end{table}

Odpowiednie szansy zakupu tego produktu w~zależności od płci klienta wynoszą:

$$
\mathbf{S\left(A\right)}=\frac{243/273}{1-243/273}=\frac{243}{273}\cdot\frac{273}{30}=\frac{243}{30}=8.1\quad
\mathbf{S\left(B\right)}=\frac{48/288}{1-48/288}=\frac{48}{288}\cdot\frac{288}{240}=\frac{48}{240}=0.2
$$

Ilorazy szans zakupu tego produktu wynoszą odpowiednio:

$$
\mathbf{OR_{AB}}=8.1/0.2=40.5 \quad \mathbf{OR_{BA}}=0.2/8.1=0.0247
$$

Iloraz szans równy $40.5$ oznacza, że szansa zakupienia produktu przez kobietę jest czterdzieści razy większa niż szansa zakupu tego produktu przez mężczyznę.

+ W~przypadku gdyby $\mathbf{OR_{AB}}=1$ oznaczałoby to, że szansa zakupu danego produktu przez kobietę jest taka sama jak szansa zakupu produktu przez mężczyznę. Innymi słowy fakt zakupu produktu nie zależałby od płci kupującego klienta.
+ W~przypadku gdy $\mathbf{OR_{AB}}>1$ oznacza to, że szansa wystąpienia danego zdarzenia w~grupie $A$ jest większa niż w~grupie $B$. 
+ W~przypadku gdy $\mathbf{OR_{AB}}<1$ oznacza to, że szansa wystąpienia danego zdarzenia w~grupie $A$ jest mniejsza niż w~grupie $B$. W~naszym przykładzie $\mathbf{OR_{BA}}=0.2/8.1=0.0247 \cong 1/40$ co oznacza, że szansa zakupienia produktu przez mężczyznę jest czterdzieści razy mniejsza niż szansa zakupu tego produktu przez kobietę.
+ Większość pakietów statystycznych (w~tym i~SPSS) wyznacza $95\%$ przedział ufności dla ilorazu szans.

Model regresji logistycznej

+ Model regresji logistycznej z~formalnego punktu widzenia jest uogólnionym modelem liniowym (GLM), w~którym użyto tzw. transformacji logitowej.
+ Model regresji logistycznej wykorzystywany jest w~sytuacji gdy zmienna objaśniana $Y$ jest zmienną dychotomiczną tj. przyjmuje tylko dwie wartości.
+ Niech $Y$ oznacza zmienną dychotomiczną o~wartościach $1$ (z~reguły odnoszących się do zdarzeń pożądanych tj. wyzdrowienie, zakup produktu, spłata kredytu itd.) oraz $0$ (w~przeciwnym wypadku, na przykład choroba, nie kupienie produktu, brak spłaty kredytu itd.). Przy tak przyjętych oznaczeniach modelem regresji logistycznej jest:

$$
P=P\left(Y=1|X_{1},X_{2},\ldots,X_{k}\right)=\frac{e^{a_{0}+a_{1}X_{1}+\ldots+a_{k}X_{k}}}{1+e^{a_{0}+a_{1}X_{1}+\ldots+a_{k}X_{k}}}
$$
		$a_{i}, i=0,\ldots,k$	oznaczają współczynniki modelu regresji logistycznej, a~$x_{1},\ldots,x_{k}$ są zmiennymi objaśniającymi (niezależnymi), które mogą być zarówno ilościowe jak i jakościowe.
		
+ Z~formalnego punktu widzenia model regresji logistycznej wyraża warunkowe prawdopodobieństwo, że zmienna $Y$ przyjmie wartość równą $1$ dla określonych wartości zmiennych objaśniających $x_{1},\ldots,x_{k}$.

+ Parametry $a_{i}, i=0,\ldots,k$ modelu regresji logistycznej szacuje się metodą największej wiarogodności (MNW), której podstawy teoretyczne sformułował Fisher. 
		
			
Transformacja logitowa

+ Transformacją logitową (logitem) nazywamy przekształcenie prawdopodobieństwa w~logarytm naturalny ilorazu szans

$$
logit P=ln\frac{P}{1-P}
$$
+ W~przypadku modelu regresji logistycznej otrzymujemy:

$$
logit P=ln\frac{P}{1-P}=ln\left(P\right)-ln\left(1-P\right)=a_{0}+a_{1}x_{1}+\ldots+a_{k}x_{k}
$$ 

+ Powyższa równość nazywana jest logitową postacią modelu regresji logistycznej i~wyraża liniową funkcję regresji wielu zmiennych.

Interpretacja parametrów modelu regresji logistycznej

+ Parametry modelu regresji logistycznej mają bardzo prostą interpretację nawiązującą do ilorazu szans (OR).
+ Niech $A$ i~$B$ oznaczają dwie porównywane ze sobą grupy. Niech ponadto $X_{A_{i}}$ oznacza wartość $i$-tej zmiennej objaśniającej dla jednostki z~grupy $A$. Podobnie niech $X_{B_{i}}$ oznacza wartość $i$-tej zmiennej objaśniającej dla jednostki z~grupy $B$. 
+ Parametry modelu regresji logistycznej interpretujemy jako ilorazy szans:

$$
OR_{AB}=e^{\displaystyle\sum_{i=1}^{k}a_{i}\left(X_{A_{i}}-X_{B_{i}}\right)}
$$

Wybrane miary oceny jakości modelu regresji logistycznej

+ Praktyczne wykorzystanie regresji logistycznej, podobnie jak w przypadku regresji wielorakiej, powinno być poprzedzone analizą jakości uzyskanego modelu. Należy m.in. sprawdzić istotność uzyskanych parametrów czy dopasowanie modelu do danych empirycznych.
+ W pakietach statystycznych jedną z najczęściej wykorzystywanych miar oceny dopasowania modelu jest statystyka ilorazu największej wiarygodności postaci:
$$
LR=-2log\left(ML_{1}\right)-\left(-2log\left(ML_{2}\right)\right)=-2log\left(ML_{1}/ML_{2}\right). 
$$
+ Statystykę tę wykorzystuje się również do porównania dwóch modeli $M_{1}$ i $M_{2}$, przy czym zakłada się, że model $M_{1}$ jest podzbiorem modelu $M_{2}$ ($ML_{1}$ i $ML_{2}$ oznaczają odpowiednio wartość funkcji wiarygodności dla modelu $M_{1}$ i $M_{2}$ odpowiednio). 
+ Statystyka ta ma asymptotyczny rozkład chi-kwadrat z liczbą stopni swobody będącą różnicą między liczbą parametrów obydwu modeli.

+ Bardzo często pierwszym krokiem weryfikacji modelu jest przyjęcie założenia, że $M_{1}$ jest modelem, w którym jedynym parametrem jest wyraz wolny, a $M_{2}$ jest pełnym modelem zawierającym wszystkie zmienne objaśniające. 
+ W hipotezie zerowej zakłada się przy tym, że wszystkie zmienne nie są istotne w sensie statystycznym wobec hipotezy alternatywnej, że istnieją pewne zmienne, które są istotne. 
+ W przypadku odrzucenia hipotezy zerowej stwierdzamy, że model z pełnym zestawem zmiennych objaśniających poprawia dopasowanie do danych empirycznych i istotnie różni się od modelu z wyrazem wolnym.
+ W przypadku gdy jesteśmy zainteresowani sprawdzeniem czy dołączenie tylko jednej zmiennej poprawia model można wykorzystać statystykę Walda postaci:
$$
Wald_{i}=\left(\frac{b_{i}}{S_{b_{i}}}\right)^{2}
$$
+ Statystyka ta ma rozkład chi-kwadrat z jednym stopniem swobody ($S\left(b_{i}\right)$  oznacza średni błąd szacunku jaki popełniamy przy szacowaniu parametru przy i-tej zmiennej objaśniającej). 
+ W przypadku odrzucenia hipotezy zerowej stwierdzamy, że włączenie  i-tej zmiennej objaśniającej poprawia w sposób statystycznie istotny dopasowanie modelu do danych empirycznych. Innymi słowy warto zostawić taką zmienną w budowanym modelu regresji logistycznej. 

+ W klasycznym modelowaniu jedną z najczęściej wykorzystywanych miar dopasowania modelu do danych empirycznych jest współczynnik determinacji $R^{2}$.
+ Pewnym odpowiednikiem współczynnika determinacji w modelu regresji logistycznej jest zaproponowany przez McFaddena współczynnik, określany w literaturze mianem pseudo $R^{2}$.
$$
R^{2}_{McFadden}=1-\frac{Ln\left(L_{p}\right)}{Ln\left(L_{0}\right)}
$$ 
			$Ln\left(L_{p}\right)$ jest logarytmem naturalnym wartości funkcji wiarygodności dla pełnego modelu a $Ln\left(L_{0}\right)$ jest logarytmem naturalnym wartości funkcji wiarygodności dla modelu z wyrazem wolnym.  
+ Nie jest to jednak współczynnik unormowany w przedziale $[0,1]$ a jego maksymalna wartość jest mniejsza od 1. 
+ Generalnie, im większa jest wartość pseudo $R^{2}$ tym większa jest zdolność trafnego przewidywania badanej zmiennej z wykorzystaniem modelu regresji logistycznej.
+ Przyjmuje się przy tym, że wartości z przedziału 0,2-0,4 są w zupełności zadowalające.
+ W literaturze występuje również szereg innych miar wykorzystywanych do badania jakości dopasowania modelu regresji logistycznej.
+ Współczynnik Coxa-Snella
$$
R^{2}_{CS}=1-exp\left(\frac{-2Ln\left(L_{0}\right)+2Ln\left(L_{p}\right)}{n}\right)
$$
+ współczynnik Nagelkerkego
$$
R^{2}_{N}=1-\frac{1-e^{-\frac{2}{n}\left(Ln\left(L_{p}\right)-Ln\left(L_{0}\right)\right)}}{1-e^{-\frac{2}{n}Ln\left(L_{0}\right)}}
$$
+ współczynnik Cragga-Uhlera

$$
R^{2}_{CU}=\frac{1-\left(\frac{L_{0}}{L_{p}}\right)^{\frac{2}{n}}}{1-\left(L_{0}\right)^{\frac{2}{n}}}
$$
+ Podobnie jak w przypadku współczynnika McFaddena wyższe wartości współczynników Coxa-Snella, Negelkerke'a i Cragga-Uhlera są bardziej pożądane, gdyż świadczą o lepszym dopasowaniu oszacowanego modelu regresji logistycznej do danych empirycznych. 
+ Szczególnie użyteczny może być współczynnik Negelkerke'a, który jest unormowany i przyjmuje wartości z przedziału $[0,1]$. 
+ Przy interpretacji współczynników dopasowania w odniesieniu do modelu regresji logistycznej należy jednak zachować szczególną ostrożność. 
+ Jak pokazują praktyczne doświadczenia, uzyskanie miar dopasowania w przypadku modelu regresji logistycznej, na poziomie 0,2-0,4 można w wielu przypadkach uznać za zadowalające. 


## Dalsze kroki