Standaryzacja

Standaryzacja polega na przekształceniu danych w taki sposób, aby ich średnia wynosiła 0, a odchylenie standardowe 1. Standaryzacja jest przydatna, gdy chcemy porównać zmienne o różnych jednostkach miary.

x <- (x - mean(x)) / sd(x)
X <- scale(X)

Usunięcie braków danych

d <- d[complete.cases(d), ]

Podział train/test

train_idx <- sample(1:n, 0.7*n, replace = F)
test_idx <- setdiff(1:n, train_idx)

TrainB <- df[train_idx,]
TestB <- df[test_idx,]

DANE TESTOWE CENTRUJEMY ŚREDNIMI Z DANYCH TRENINGOWYCH

Badanie istotności zmiennej z użyciem hipotezy liniowej

Założmy pewien podział naszych obserwacji na 3 kategorie. Chcemy sprawdzić, czy zmienna predyktorinho jest tak samo istotna dla każdej z kategorii. Chcemy przetestować, czy dla modelu gdzie rozbijemy kolumnę predyktorinho na 3 oddzielne kolumny przyjmujące wartości \(\neq 0\) tylko dla zmiennych przynależnych do danej kategorii, współczynniki przy tych kolumnach są równe.

i1 = ifelse(c$Group == 1, 1, 0)
i2 = ifelse(c$Group == 2, 1, 0)
i3 = ifelse(c$Group == 3, 1, 0)
i <- cbind(i1,i2,i3)
d <- data.frame(i1, i2, i3, i*c$predyktorinho, c$Y)
# z 2 kolumn zrobiliśmy 6
# kolumnę Group rozbiliśmy na 3 binarne kolumny dot. przynależności do grupy
# kolumnę z predyktorem rozbiliśmy na 3 kolumny, gdzie każda kolumna
# przyjmuje prawdziwe wartości dla danej grupy, 0 wp.p.
names(d)[4:7] <- c("s1", "s2", "s3", "Y")

Musimy wyrzucić z modelu Intercept jeżeli korzystamy z n kolumn opisujących przynależność do n grup (de facto wystarczyłoby n-1 kolumn).

modelinho <- lm(Y ~ .-1, data = d)
# testujemy hipotezę liniową, że współczynniki przy kolumnach s1, s2, s3 są wszystkie równe
lht(mod, c("s1 = s2", "s2 = s3"))
# p-value wyniosło 0.075, więc nie możemy odrzucić hipotezy że współczynniki są równe

Ridge: teoria

Ogólnie metoda Ridge (regresji grzbietowej) polega na regularyzacji funkcji straty poprzez dodanie do niej kary w normie L2.

Estymator regresji grzbietowej \(\hat{\beta}_{r}\) jest rozwiązaniem problemu optymalizacyjnego:

\[ \hat{\beta_{0,r}}, \hat{\beta}_{-0,r}^T = \arg \min_{b_0,b_{-0} } \left\{ \sum_{i=1}^n (y_i - \beta_0 - \sum_{j=1}^{p-1} x_{ij} \beta_j)^2 + \lambda \sum_{j=1}^{p-1} \beta_j^2 \right\} \]

\(\lambda > 0\) jest hiperparametrem, który na ogół jest dostrajany za pomocą crossvalidation.

Dla \(\lambda = 0\) otrzymujemy estymator najmniejszych kwadratów, zwiększanie \(\lambda\) powoduje zmniejszanie się wartości współczynników \(\beta_j\). Dla \(\lambda \to \infty\) współczynniki \(\beta_j\) dążą do zera - są ściągane wewnątrz kuli o promieniu \(\lambda\).

Estymator Ridge jest lepszy, jeżeli zmienne są ze sobą skorelowane. Postać estymatora ridge wygląda (dla problemu bez wyrazu wolnego) następująco: \[ \hat{\beta}_{r} = (X^T X + \lambda I)^{-1} X^T Y \] Widać zatem, że jeżeli \(X'X\) jest słabo uwarunkowana (tzn jej najmniejsza wartość własna jest bliska zero), to estymator MNK postaci \((X^T X)^{-1} X^T Y\) jest źle uwarunkowany (bierzemy odwrotność czegoś co jest prawie osobliwe). Jeżeli wprowadzamy regularyzacje, to do każdej z wartości własnych dodajemy \(\lambda\), co powoduje, że estymator jest lepiej uwarunkowany - macierz \((X^T X + \lambda I)\) jest zawsze odwracalna.

Wartość wyrazu wolnego w modelu regresji liniowej to: \[ \hat{\beta}_{0,r} = \bar{y} - \sum_{j=1}^{p-1} \hat{\beta}_{j,r} \bar{x}_j = \bar{y} - \bar{x}'_{-0} \hat{\beta}_{-0,r} \]

Aby znaleźć \(\hat{\beta}_{-0, r}\) należy zastosować algorytm:
1 Centrujemy \(Y\) i standaryzujemy predyktory \(X\) ( tzn od \(Y\) odejmujemy średnią elementwise, każdą kolumnę \(X\) odejmujemy średnią i dzielimy przez normę po odjęciu). Wtedy \(\hat{\beta}_{0,r} = 0\).
2. Rozwiązujemy problem bez wyrazu wolnego
3. skalujemy na powrót do oryginalnej skali
4. liczymy \(\hat{\beta}_{0,r}\) ze wzoru wyżej

Użyteczny może być też tutaj rozkład SVD macierzy \(X = UDV^T\): \[ (X'X)^{=1} = VD^{-2}V' \]

\[ \hat{Y}_r = U(D^2 + \lambda I)^{-1} D^2 U^T Y \]

Ogólnie \(\hat{Y}_{MNK}\) to rzut Y na przestrzeń generowaną przez kolumny macierzy U (czyli kolumny X, ale zortogonalizowane), \(\hat{Y}_r\) to rzut Y na przestrzeń generowaną przez kolumny macierzy U, ale ściągnięte do środka układu współrzędnych.

Ridge jest obciążony, ale macierz kowariancji jest mniejsza niż MNK: \[ \Sigma_r \leq \Sigma_{MNK} \] (macierz można łatwo policzyć ze wzoru \(\Sigma_{AY} = A\Sigma_YA'\))

Własności Ridge

1. Własność przesunięcia \[ Y_i = Y_i + c => \hat{\beta}_{-0,i} = \hat{\beta}_{-0,i}, \hat{\beta}_0 = \hat{\beta}_0 + c \]
2. Własność skalowania

\[ X_i = cX_i => \hat{\beta}_i = \frac{\hat{\beta}_i}{c} \]

żeby osiągnąć 1. kara z regularyzacji nie jest dla wyrazu wolnego, aby osiągnąć 2. stosujemy normalizacje

Lasso: teoria

W przypadku lasso kara w regularyzacji jest w normie L1, a nie L2. Estymator lasso \(\hat{\beta}_l\) jest rozwiązaniem problemu optymalizacyjnego: \[ \hat{\beta}_{0,l}, \hat{\beta}_{-0,l}^T = \arg \min_{b_0,b_{-0} } \left\{ \sum_{i=1}^n (y_i - \beta_0 - \sum_{j=1}^{p-1} x_{ij} \beta_j)^2 + 2\lambda \sum_{j=1}^{p-1} |\beta_j| \right\} \]

Podobnie jak w ridge centrujemy \(Y\) i standaryzujemy kolumny \(X\), więc na ogół rozpatrujemy lasso bez wyrazu wolnego. Wtedy estymator lasso ma postać: \[ \hat{\beta}_{-0,l} = \arg \min_{b_{-0} } \left\{ (Y- X\beta)'(Y-X\beta) + 2\lambda |\beta|_1 \right\} \]

Lasso działa troche jak selektor zmiennych - zwiększając \(\lambda\) zeruje kolejne współczynniki \(\beta_j\).

Optymalizacja nieróżniczkowalnych wypukłych funkcji - do teorii

Funkcja celu lasso jest wypukła, ale nieróżniczkowalna. Aby poradzić sobie z tym problemem wprowadzamy pojęcie subgradientu. Subgradient \(F: \mathbb{R}^p \rightarrow \mathbb{R}\) w punkcie \(x_0\) to dowolny wektor \(g \in \mathbb{R}^p\) taki, że: \[ F(x) \geq F(x_0) + g^T(x-x_0) \] Subróżniczka \(\partial F(x_0)\) to zbiór wszystkich subgradientów w punkcie \(x_0\), innymi słowy jest to zbiór wszystkich współczynników \(\alpha\) takich, że prosta przechodząca przez punkt \((x_0, F(x_0))\) i o kierunku \(\alpha\) leży poniżej wykresu funkcji \(F\) w sąsiedztwie \(x_0\).

Subróżniczka \(F(x) = |x|\) w 0 to zbiór \(\{a: |a \leq 1, ax |x|\}\) Wobec tego \[ (\partial |x|_1)_j = \begin{cases} -1 & x_j < 0 \\ [-1,1] & x_j = 0 \\ 1 & x_j > 0 \end{cases} \] co można również zapisać jako: \[ \partial |x|_1 = \{g: |g|_{\infty} \leq 1, g^Tx = |x|_1\} \] Minimalizacja funkcji celu lassso jest wobec tego osiągalna wtedy i tylko wtedy gdy \[ 0 \in \partial \left\{ \sum_{i=1}^n (y_i - \sum_{j=1}^{p-1} x_{ij} \beta_j)^2 + 2\lambda \sum_{j=1}^{p-1} |\beta_j| \right\} = \partial( (Y- X\beta)'(Y-X\beta) + 2\lambda |\beta|_1) = \{2X'X\beta - 2X'Y + 2\lambda g: |g|_{\infty} \leq 1, g^Tx = |x|_1\} \] łatwo sprawdzić, że istotnie 0 należy do tego zbioru, wystarczy tylko wziąć odpowiednie \(g\).

Wybór współczynników

Lasso, w związku z karą w normie L1, ma ciekawą własność zerowania współczynników dla coraz większych wartości parametru \(\lambda\). Tę własność można zobaczyć na wykresie współczynników od wartości \(log(\lambda)\).

plot(m, xvar="lambda", label=TRUE)
coefplot::coefpath(m) # interaktywny wykres

coef(m, s=0.1) # dla lambda = 0.1

m.cv <- glmnet::cv.glmnet(X, y, alpha=1)

plot(m.cv)

m.cv$lambda.min
m.cv$lambda.1se

lasso_1se <- glmnet::glmnet(X, y, alpha=1, lambda=m.cv$lambda.1se)
# ewentualnie
predict(m.cv, type = "coef", s = m.cv$lambda.1se)

PCA

Principal Component Analysis - Analiza Głównych Składowych, polega na wyszukiwaniu kierunków które dadzą nam największą wariancje.

Dla kierunku \(||a||_2 = 1\) maksymalizujemy wariancję \(a^TX\): \[ \text{Var}(a^TX) = a^T \Sigma_X a \] czyli rozwiązujemy zadanie: \[ \arg\max_{a: ||a||_2 = 1} a^T \Sigma_X a = v_1 \] gdzie \(v_1\) to wektor własny macierzy $_X $ odpowiadający największej wartości własnej.

To be precise, w przypadku gdy nie mamy dostępu do macierzy kowariancji \(\Sigma_X\) to maksymalizujemy wariancję \(a^TS\), gszie \(S\) to próbkowa macierz kowariancji.

\[ S = \frac{1}{n-1} \sum_{i=1}^n (x_i - \bar{x})(x_i - \bar{x})^T \] Dla macierzy \(X\) z centrowanymi kolumnami, czyli \(\bar{x} = 0\) mamy: \[ S = \frac{1}{n-1} X^TX \]

Następne kierunki wyznaczamy w analogiczny sposób, ale tak aby były ortogonalne do poprzednich.

Aby łatwiej wyznaczyć kierunki główne z PCA wystarczy wziąć kolumny macierzy \(V\) z rozkładu spektralnego macierzy \(X\), wtedy \(X'X = VD^2V'\).

PCR

Principal Component Regression - Regresja z wykorzystaniem PCR, polega na zastąpieniu zmiennych pierwotnych \(r\) kierunkami głównymi wyznaczonymi przez PCA, a następnie przeprowadzeniu regresji liniowej na tych zmiennych.

model_pcr <- pcr(y ~., data = data, scale = T, validation = "CV")

validationplot(model_pcr, val.type = "RMSEP")
validationplot(model_pcr, val.type = "R2")

pcr_pred <- predict(model_pcr, test, ncomp = Variable)
# Variable to liczba kierunków głównych

model_pcr <- pcr(medv ~ ., data = train, scale = TRUE, validation = "CV")

# X to ramka danych zmiennych objaśniających, Y to kolumna zmiennej objaśnianej
X <- scale(X)

svd <- svd(X)
V <- svd$v

Z <- X %*% V # składowe główne
pcr_data <- data.frame(y=Y, Z)

m_pcr <- lm(y ~ ., data = pcr_data)

# jeszcze bardziej ręcznie - liczymy współczynniki
p <- ncol(pcr_data)
theta <- rep(0, p)
theta[1] <- mean(Y)

for (i in 2:p) {
    theta[i] <- pcr_data$medv %*% pcr_data[,i] / (pcr_data[,i] %*% pcr_data[,i])
}

# ostatecznie jeszcze trzeba wybrać ile współczynników używamy

# współczynniki przy beta - orginalny estymator
r <- 13

beta_pcr <- rep(0, p-1)
for (i in 1:r) {
  beta_pcr = beta_pcr + theta[i+1]*V[,i]
  
}
y_pred <- X %*% t(t(beta_pcr)) + theta[1]

PLSR

PCR nie uwzględnia zmiennej celu. PLSR już tak. Dla scentrowanej macierzy \(X\) wykonuje podobne operacje co PCR, tylko że nie rzutujemy na przestrzeń kolumn \(X\).

mod_plsr <- plsr(y ~., data = data, scale = T, validation = "CV")

validationplot(mod_plsr, val.type = "MSEP")
validationplot(mod_plsr, val.type = "R2")

plsr_pred <- predict(model_plsr, test, ncomp = Variable)
# Variable to liczba kierunków głównych

# y to wektor, X to macierz

X <- scale(X)

n <- nrow(X)
e <- rep(1, n)
r <- 10 # liczba krokow

y_plsr <- mean(y) * e # y startowy
X_plsr <- X
fi <- rep(0, p)
z <- matrix(0, n, r)
theta <- rep(0, r)

for (m in 1:r){
  
  fi <- as.vector((t(y) %*% X_plsr) / diag((t(X_plsr) %*% X_plsr)))
  z[, m] <- fi %*% t(X_plsr)
  
  theta[m] <- t(y) %*% z[, m] / (t(z[, m]) %*% z[, m])
  y_plsr <- y_plsr - theta[m] * z[, m]
  
  X_plsr <- X_plsr - z[, m]%*% t((t(X_plsr) %*% z[, m]) / c(t(z[, m]) %*% z[, m]))
  X_plsr <- scale(X_plsr, center = FALSE, scale = TRUE)
}


pred <- e*mean(y) + z %*% (theta)

Wiarogodność w modelu logistycznym

prawdopodobieństwo w modelu logistycznym wyraża się jako \[ f_i(Y_i|x_i) = \pi_i^{y_i} (1-\pi_i)^{1-y_i} \] gdzie \(\pi_i = P(Y_i = 1 | X_i = x_i) = \sigma(\beta'x_i) = \sigma(\beta_0 + \beta_1 x_{i1} + \dots + \beta_p x_{ip})\)

\[ L(\beta) = f(Y_1, Y_2, \dots, Y_n | X_1 = x_1, \dots, X_n=x_n) = \prod_{i=1}^n p_i^{y_i} (1-p_i)^{1-y_i} \]

Funkcja log-wiarogodności: \[ \mathbf{L}(\beta) = \sum_{i=1}^n y_i \beta'x_i + \sum_{i=1} \log (1 + exp(\beta'x_i)) \] jest to funkcja wklęsła parametru \(\beta\). \[ \partial \log L = X'(Y- \pi) \] gdzie \(\pi = (\pi_1, \dots, \pi_n)'\).

Liniowa separowalność klas - efekt Haucka-Donnnera

Wykres pozwalający dla prostego przypadku zidentyfikować liniową separowalność klas.

ggplot(data=data,
       aes(
         x=var_x,
         y=var_y,
         col=category
       )) + geom_point()

W wypadku idealnej separowalności klas, czyli gdy istnieje hiperpłaszczyzna rozdzielająca klasy, to mamy do czynienia z efektem Haucka-Donnnera. Algorytm odpowiedzialny za stworzenie modelu regresji liniowej nie zbiega, tzn. cwel próbuje stworzyć funkcję bliską funkcji indykatorowej, zatem współczynniki \(\beta\) stają się kurewsko duże. Dymy polegają na tym, że \(SE(\beta)\) dąży do \(\infty\) dużo szybciej, co za tym idzie statystyka testowa Walda postaci \(\frac{\beta}{SE(\beta)}\) dąży do 0. Wtenczas p-value jest bliskie 1, czyli \(\textbf{model leci w Variablea}\), że zmienne są zupełnie nieistotne (mimo że są istotne).

Szansa (odds)

jest to \(\frac{p}{1-p}\), gdzie \(p\) to prawdopodobieństwo sukcesu.

Test Walda

można znaleźć w summary od modelu, jest to z value, pwartość oznaczamy tu jako Pr(>|z|). Test Walda (z-value) to wartość estymatora przez błąd standardowy

Test dewiancji

testowana jest funkcja log-wiarogodności ( bo estymator \(\beta\) jest konstruowany jako log-wiarogności):

\(H_0\): mniejszy model \(\omega\) jest adekwatny \(H_1\): większy model \(\Omega\) jest adekwatny, a mniejszy nieadekwatny

\[ D_{\omega. \Omega} = 2\log \frac{L(\beta_\Omega)}{L(\beta_\omega)} \geq 0 \]

Jeżeli \(H_0\) jest prawdziwa, to \(D_{\omega. \Omega}\) ma rozkład \(\chi^2\) z \(df = p - q\) stopniami swobody, gdzie \(p\) to liczba parametrów w modelu, a \(q\) to liczba parametrów w modelu H\(H_0\).

Żeby sprawdzić czy zbiór zmiennych jest istotny porównujemy model vs model tylko z interceptem (null) Istotność pojedynczej zmiennej przy założeniu istnienia pozostałych sprawdzamy porównując model bez tej zmiennej do takiego z.

Residual deviance to \(D_{\omega_{n}. \omega}\) gdzie \(\omega_{n}\) to model nasycony - czyli taki co jest adekwatny i pełny, \(\omega\) to nasz model co go testujemy, a null deviance to \(D_{\omega_n. \omega_0}\) gdzie \(\omega_0\) to model z samym interceptem.

residual deviance używamy tylko do danych grupowych, gdy liczba obserwacji dla ustalonej wartości zmiennej to co najmniej 5

dane grupowe to takie:

dane <- data.frame(
  dead = c(2, 1, 1),
  alive = c(2, 3, 1),
  V1 = c(1, 5, 6),
  V2 = c(2, 3, 4)
)

dane indywidualne:

dane <- data.frame(
  y <- c(1, 0, 1, 0, 1, 0, 1, 0, 1, 0),
  V1 <- c(1, 1, 1, 2, 5, 5, 5, 5, 6, 6),
  V2 <- c(2, 2, 2, 1, 3, 3, 3, 3, 4, 4)
)

Porównanie modeli regresji logistycznej

Porównujemy dwa modele regresji logistycznej, podobnie jak z modelami liniowymi oraz hipotezą liniową. \(H_0\) jest taka, że model mniejszy jest lepszy.

SA.logit <- glm(Y ~ ., data = data, family = "binomial")
SA.logit_null <- glm(Y ~ 1, data = data, family = "binomial") # model pusty

anova(SA.logit_null, SA.logit, test = "Chisq") # tu od razu p-wartosc
pv = 1 - pchisq(123.97, df = 9)

Log-wiarogodność oraz dewiancja z definicji

sum <- 0
n <- nrow(d[train, ])

log_wiarogodnosc <- function (y, model) {
  sum <- 0
  n <- nrow(d[train, ])
  for (i in 1:n) {
    sum <- sum + c(y[i]) * c(coef(model)) %*% model.matrix(model)[i, ] - log(1 + exp(coef(model) %*% model.matrix(model)[i, ]))
  }
  return(sum)
}

aic <- function (y, model) {
  return(-2*log_wiarogodnosc(y, model) + 2*length(coef(model)))
}

deviance_test <- function (y, model1, model2) {
  return(2*(log_wiarogodnosc(y, model1) - log_wiarogodnosc(y, model2)))
}