Punkt 1 - Podział danych

Zamiast dzielić dane na dwie grupy, została przygotowana funkcja zwracająca indeksy danych spełniające warunki podziału. W ten sposób będzie można łatwo wybrać dane treningowe dane[trening,] oraz dane testowe dane[-trening].

trening = PodzielDane(dane, "STAN", podział, delta)

Frakcja w danych testowych wynosi 0.876, podczas gdy frakcja wszystkich danych wynosi 0.869.

Punkt 2 - Konstrukcja modelu regresji logistycznej

model = glm(STAN~., data = dane[trening,], family = binomial("logit"))

Do konstrukcji modelu została użyta standardowa funkcja glm.

Punkt 3 - Weryfikacja hipotezy \(H_0:\beta_2=...=\beta_p=0\)

Do zweryfikowania hipotezy zerowej \(H_0:\beta_2=...=\beta_p=0\) korzysta się ze statystyki dewiancji D, która ma rozkład \(\chi^2\) z \(k\) stopniami swobody. Do wyliczenia wartości dewiancji została użyta funkcja anova.

D = anova(model, update(model, ~1), test="Chisq")

Wartość dewiancji w tym przypadku wynosi -25.361, a prawdopodobieństwo testowe 0.003.Zatem należy odrzucić hipotezę zerową o braku liniowej zależności pomiędzy zmiennymi objaśniającymi, a zmienną objaśnianą.

Punkt 4 - Wyznaczenie współczynników regresji wraz z ich błędem standardowym oraz ich przedziałów ufności

summ(model, digits = 3)
Observations 186
Dependent variable STAN
Type Generalized linear model
Family binomial
Link logit
χ²(9) 25.361
Pseudo-R² (Cragg-Uhler) 0.242
Pseudo-R² (McFadden) 0.182
AIC 133.822
BIC 166.079
Est. S.E. z val. p
(Intercept) -8.085 3.444 -2.348 0.019
RMS10 -28.154 22.642 -1.243 0.214
RMS20 39.952 24.124 1.656 0.098
RMS30 -15.074 15.304 -0.985 0.325
PWD 6.908 2.097 3.294 0.001
A 0.001 0.011 0.120 0.904
DD -1.539 3.552 -0.433 0.665
YA 0.012 0.022 0.523 0.601
YDD 5.667 4.880 1.161 0.246
AR 0.232 1.572 0.148 0.883
Standard errors: MLE

Do wyznaczenia wartości estymatorów współczynników regresji \(\hat\beta_i\) oraz ich błędów standardowych została użyta funkcja summ. Jednak w przypadku regresji logistycznych, zamiast posługiwania się współczynnikami regresji \(\beta_i\), które trudno interpretować, lepiej jest przedstawić je w wartościach wykładniczych, które oznaczają zmianę szansy dla danego czynnika. Dla współczynników o wartościach mniejszych od 1 dany czynnik ma wpływ ograniczający, a dla wartości większej od 1 ma wpływ stymulujący.

summ(model, exp = TRUE, digits = 3)
Observations 186
Dependent variable STAN
Type Generalized linear model
Family binomial
Link logit
χ²(9) 25.361
Pseudo-R² (Cragg-Uhler) 0.242
Pseudo-R² (McFadden) 0.182
AIC 133.822
BIC 166.079
exp(Est.) 2.5% 97.5% z val. p
(Intercept) 0.000 0.000 0.263 -2.348 0.019
RMS10 0.000 0.000 11100082.412 -1.243 0.214
RMS20 224380767182166144.000 0.001 76801601961981899840046688488262880246.000 1.656 0.098
RMS30 0.000 0.000 3021361.624 -0.985 0.325
PWD 1000.494 16.402 61028.463 3.294 0.001
A 1.001 0.979 1.024 0.120 0.904
DD 0.215 0.000 226.416 -0.433 0.665
YA 1.012 0.968 1.057 0.523 0.601
YDD 289.038 0.020 4117572.940 1.161 0.246
AR 1.261 0.058 27.495 0.148 0.883
Standard errors: MLE

Tak samo w wartościach wykładniczych można przedstawić przedział ufności tych współczynników, co widać w powyższej tabeli. Aby zweryfikować hipotezę \(H_0:\beta_i = 0\) vs \(H_0:\beta_i \ne 0\) został zastosowany test Walda. Prawdopodobieństwo testowe dla tego testu jest również przedstawione w powyższych dwóch tabelach. Można z nich odczytać, że hipotezę zerową należy odrzucić tylko w przypadku wyrazu wolnego (Intercept) oraz w przypadku zmiennej PWD. Dla pozostałych predyktorów nie możemy, na przyjętym poziomie istotności \(\alpha\), odrzucić hipotezy o braku liniowej zależności pomiędzy nimi a zmienną objaśnianą.

Punkt 5 - Stworzenie macierzy konfuzji oraz wyznaczenie parametrów predykcyjnych modelu

Do wyznaczenia macierzy konfuzji oraz wszelkich parametrów predykcyjnych modelu, została wykorzystana własna funkcja LogParam.

Param0.5 = LogParam(model, dane[-trening,], 0.5)
Param0.7 = LogParam(model, dane[-trening,], 0.7)

Macierz konfuzji (ang. Classification Table) dla poziomu odcięcia \(\pi_0=0.5\).

kable(Param0.5 %>% select(TP:FN), format = 'html') %>%
  kable_styling(bootstrap_options = "striped", full_width = F)
TP TN FP FN
65 1 11 4

Parametry predykcyjne modelu dla poziomu odcięcia \(\pi_0=0.5\).

kable(Param0.5 %>% select(Acc:Err), digits = rep(3, 8), format = 'html') %>%
  kable_styling(bootstrap_options = "striped", full_width = F)
Acc Ppv Npv Tnr Fnr Tpr Fpr Err
0.815 0.855 0.2 0.333 0.058 0.942 0.917 0.185

gdzie:

Macierz konfuzji (ang. Classification Table) dla poziomu odcięcia \(\pi_0=0.7\).

kable(Param0.7 %>% select(TP:FN), format = 'html') %>%
  kable_styling(bootstrap_options = "striped", full_width = F)
TP TN FP FN
61 3 9 8

Parametry predykcyjne modelu dla poziomu odcięcia \(\pi_0=0.7\).

kable(Param0.7 %>% select(Acc:Err), digits = rep(3, 8), format = 'html') %>%
  kable_styling(bootstrap_options = "striped", full_width = F)
Acc Ppv Npv Tnr Fnr Tpr Fpr Err
0.79 0.871 0.273 0.667 0.116 0.884 0.75 0.21

Punkt 6 - Optymalny próg odcięcia

W celu znalezienia optymalnego progu odcięcia, można posłużyć się wykresem odsetek prawidłowych klasyfikacji w funkcji progu odcięcia. Pomocny jest tutaj także wykres ROC.

DaneTest = dane[-trening,]
ROCpred = prediction(predict(model, DaneTest, type = "response"),
                     DaneTest$STAN)
ROCperf = performance(ROCpred, 'tpr', 'fpr')
plot(ROCperf, colorize = TRUE, text.adj = c(-0.2, 1.7),
     main = "Wykres ROC")

ROCperf = performance(ROCpred, 'acc')
plot(ROCperf, 
     main="Odsetek prawidłowych klasyfikacji w funkcji progu odcięcia",
     xlab="Próg odcięcia",
     ylab="Odsetek prawidłowych klasyfikacji")

cost.perf1 = performance(ROCpred, "cost", cost.fp = 1, cost.fn = 10)
cost.perf2 = performance(ROCpred, "cost", cost.fp = 10, cost.fn = 1)

Analizując powyższe wykresy należy zauważyć, że im niższy próg odcięcia, tym wyższa wartość współczynnika Acc. Dodatkowo można to potwierdzić wyznaczając optymalny próg odcięcia dla minimalnego kosztu błędu. Przy założeniu, że koszt błędu FN jest dziesięciokrotnie większy od błędu FP, optymalny próg odcięcia wynosi 0.331. Jeżeli natomiast odwrotnie - koszt błędu FP jest dziesięciokrotnie większy od błędu FN, to optymalny próg odcięcia wynosi 0.942.

Funkcje dodatkowe

Poniżej został umieszczony kod funkcji dodatkowych.

PodzielDane = function(dane, y, podział = 0.8, delta = 0.025){
  frak0 = sum(dane[y])/nrow(dane)
  repeat{
    wybrane = sample(seq(nrow(dane)), podział*nrow(dane))
    frak = sum(dane[wybrane,y])/length(wybrane)
    if (abs(frak0 - frak) <= delta) break
  }
  wybrane
}


LogParam = function(model, DaneTest, poziom = 0.5){
  df = tibble(test = DaneTest[[1]],
              p = predict(model, DaneTest, type = "response"),
              resp = ifelse(p>poziom, 1, 0))
  wynik = tibble(n = nrow(df),
                 P = sum(df$test == 1),
                 N = sum(df$test == 0),
                 TP = sum(df$test == 1 & df$resp == 1),
                 TN = sum(df$test == 0 & df$resp == 0),
                 FP = sum(df$test == 0 & df$resp == 1),
                 FN = sum(df$test == 1 & df$resp == 0),
                 Accu = (TP + TN)/(P+N),   #overall proportion of correct classifications
                 Ppv = TP/(TP + FP),
                 Npv = TN/(TN + FN),
                 Tnr = FN/N,  #specifity
                 Fnr = FN/P,  #false negative rate
                 Tpr = TP/P,  #sensivity
                 Fpr = FP/N,  #false positive rate
                 Err = (FP+FN)/(P+N))
  wynik
}