Generowanie macierzy o kolumnach z danego rozkładu

Przykład macierzy 10x2 z kolumnami z rozkładu normalnego:

X = matrix(0, nrow=10, ncol=2)
X[,1] = rnorm(10,0,1)
X[,2] = rnorm(10,0,1)

Generowanie wektora zer

Wektor zer o długości 10

rep(0, 10)

Generowanie punktów z rozkładu

Generowanie chmury punktów \(y_{i} = x_{i} + \varepsilon\), gdzie \(x_{i}\) to punkty z zakresu [0; 10] odległe o 0:1, natomiast \(\varepsilon_{i}\) to zmienne losowe z rozkładu normalnego o średniej \(\mu\) = 0 i odchyleniu standardowym = 3

xi <- seq(0, 10, by = 0.1)
yi <- xi + rnorm(length(xi), mean=0, sd=3)

Wgrywanie i wyświetalnie danych

M <- as.matrix(read.csv("obrazek.csv"))
image(M,  asp = TRUE, col = c("white", "black"), xaxt = "n", yaxt = "n")

par(mfrow = c(2, 2)) # c(liczba wierszy, liczba kolumn)

Rozkład SVD (Singular Value Decomposition)

\(X = U D V^T\)

svd_decomp = svd(X)

U = svd_decomp$u
D = diag(svd_decomp$d) # domyślnie D zwracane jest jako wektor
V = svd_decomp$v

all.equal(U %*% D %*% t(V), X)

Ten rozkład można wykorzystać np do kompresji obrazków, przez uwzględnienie jedynie części wartości osobliwych na diagonali. Poniższa funkcja właśnie w ten sposób konwertuje macierz pikseli \(M\).

generate_compressed <- function(M, compression_level=0.01) {
  svd_zebra <- svd(M)
  
  U <- svd_zebra$u
  V <- svd_zebra$v
  d <- svd_zebra$d
  
  number_of_observations <- floor(length(d)*compression_level)
  
  diag <- diag(c(d[1:number_of_observations], numeric(length(d) - number_of_observations)))
  
  compressed <- U %*% diag %*% t(V)
}

Rozkład spektralny

\(Σ = V D V^T\)

Σ = t(X) %*% X

eigen_decomp = eigen(Σ)
V = eigen_decomp$vectors
D = diag(eigen$values)

all.equal(V %*% D %*% t(V), Σ)

Rozkład QR

\(X = QR\), gdzie \(Q\) to macierz ortogonalna (\(Q^T = Q^{-1}\)), a \(R\) macierz trójkątna górna.

Głównie wykorzystywaliśmy go do sprawniejszego wyznaczania estymatora \(\hat \beta_{MNK}\) Metoda ta jest wygodna, gdy macierz \(X\) jest nieodwracalna, czyli posiada liniowo zależne kolumny.

QR <- qr(X)

R <- qr.R(QR)
Q <- qr.Q(QR)

# beta
beta <- backsolve(R, t(Q) %*% y)

Adjusted (skorygowane) \(R^2\)

\(R^2_{adj} = 1 - \frac{SSE/(n-|i|)}{SST/(n-1)} = 1 - \frac{\hat\sigma_i}{\hat\sigma^2_{null}}\) gdzie \(i\) to zmienne w modelu, model pełny \(|i|=p\)

nie ma uzasadnienia metodologicznego, ale jest używane w praktyce nie jest rosnące w zależności od liczby zmiennych tak jak \(R^2\), może zostać w ostateczności wykorzystane do selekcji zmiennych

Mozna to liczyc z definicji albo wziac z summary(mod)

# wzór na adj R^2
SST = sum((y - mean(y))^2)        
SSE = sum((y - mod$fitted.values)^2)
n=nrow(X)
p=ncol(X)
R2_adj <- 1-((SSE/(n-p))/(SST/(n-1)))

Stworzenie modelu oraz predykcja

modelinho = lm(y ~ x) # model liniowy

# sposób I
predict(modelinho, newdata = data.frame(x = 21.37))

# sposób II
coef(modelinho)[1] + coef(modelinho)[2]*21.37

Macierz eksperymentu

# Model liniowy gdzie zmienna Variable opisywana jest przez wszystkie pozostałe zmienne
modelinho = lm(Variable ~ ., data = dane)
# Macierz eksperymentu
X = model.matrix(modelinho)

Sposób bez modelu:

X # nasza macierz bez interceptu
intercept <- rep(1, nrow(X))
x_intercept <- cbind(intercept, X)
x_intercept <- as.matrix(X)

Interpretacja summary()

lm(formula = Price ~ ., data = dane)

Residuals:
     Min       1Q   Median       3Q      Max 
-12.7630  -4.0514   0.5389   2.3899  12.9855 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) 13.712572   9.514111   1.441   0.1677  
Bedroom     -7.756208   3.109374  -2.494   0.0232 *
Space        0.011626   0.008981   1.295   0.2128  
Room         5.097706   2.764303   1.844   0.0827 .
Lot          0.228063   0.195434   1.167   0.2593  
Tax          0.003374   0.006859   0.492   0.6291  
Bathroom     5.718372   4.276867   1.337   0.1988  
Garage       3.613603   2.064997   1.750   0.0982 .
Condition   -2.162027   4.137400  -0.523   0.6080  

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.337 on 17 degrees of freedom
Multiple R-squared:  0.7688,    Adjusted R-squared:   0.66 
F-statistic: 7.065 on 8 and 17 DF,  p-value: 0.0003757

• Residuals: rozkład rezyduów. Powinniśmy mieć medianę bliską zera, rozkład symetryczny \((Q3-Q2 \approx Q2-Q1)\). Średnia w modelu z interceptem zawsze równa 0
  

• Estimate: w MNK każda \(\beta_{i}\) ma rozkład normalny. “Estimate” to estymatory każdej z \(\beta_{i}\) (i tym samym średnie rozkładu)
  

• Std. Error: pierwiastek z wariancji rozkładu estymatora
  

• t value: wartość statystyki t, służącej do testowania dla każdego współczynnika hipotezy \(\mathbb{H}_{0}\), że dany współczynnik jest rowny 0
  

• Pr(>|t|): p-wartość testu; jeśli poniżej poziomu istotności, to możemy stwierdzić, że istotnie współczynnik jest różny od zera
  

• RSE: podniesiony do kwadratu daje estymator wariancji naszego błędu, jako że zakładamy, że \(\mathcal{E} \sim \mathbf{N}(0, \sigma^2)\)
  

• R^2: współczynnik determinancji (procent wyjaśnionej wariancji)
  

• Adjusted R^2: dostosowany \(R^2\), który jest odpowiednio pomniejszony ze względu na liczbę zmiennych w modelu (ponieważ zwiększenie ich liczby może jedynie zwiększać \(R^2\))
  

• F-statistic: wartość statystyki F, służącej do testowania \(\mathbb{H}_{0}\), że wszystkie współczynniki \(\beta_{i}\) są równe 0. Alternatywna hipoteza to taka, że przynajmniej jeden z nich jest niezerowy; de facto porównanie naszego modelu z modelem zerowym, w którym jest tylko intercept

Estymator \(\hat \beta_{MNK}\)

Szukamy takiego wektora \(\hat \beta_{MNK}\), który minimalizuje sumę kwadratów błędów: \[ Crit(b) = ||Y-Xb||^2 = (Y-Xb)^T(Y-Xb) = \sum_{i=1}^n |Y_i - \sum_{j=0}^{p-1} b_j x_{ij}|^2\]

Interpretacja geometryczna jest taka, że szukamy w hiperpłaszczyznie rozpiętej na kolumnach wektora \(X\) takiego wektora \(Xb_{MNK}\), który jest najbliżej wektora \(Y\). że jego odległość od wektora \(Y\) jest najmniejsza.

modelinho = lm(y ~ ., data = dane)
X = model.matrix(modelinho)
# X to macierz eksperymentu

solve(t(X) %*% X, t(X) %*% y) # parametry z definicji
coef(mod) # parametry z funkcji wbudowanej

Wartość \(\beta\) można wyznaczyć również za pomocą rozkładu \(QR\) macierzy \(X\).

X <- model.matrix(modelinho)
QR <- qr(X)

R <- qr.R(QR)
Q <- qr.Q(QR)


beta <- backsolve(R, t(Q) %*% y)

y_hat <- Q %*% t(Q) %*% y

Wyjaśnienie: \[ \beta = (X^TX)^{-1}X^Ty = (QR)^{-1}X^Ty = (R^TQ^TQR)^{-1}X^Ty = (R^TR)^{-1}R^TQ^Ty = R^{-1}Q^Ty \]

Dodatkowo \[ R^{-1} = (R^TR)^{-1}R^T \]

Kilka faktów

Na ogół zakładamy, że nasz model regresji liniowej ma postać: \[ Y = X\beta + \varepsilon \] gdzie \(\varepsilon \sim \mathbf{N}(0, \sigma^2I)\)
Estymator MNK jest nieobciążony, czyli \(E(\hat \beta_{MNK}) = \beta\)
Wariancja MNK wynosi \(\Sigma_{\hat \beta_{MNK}} = \sigma^2(X'X)^{-1}\)

Wariancja \(\hat \beta\)

Wariancja \(\hat \beta\) to \(\sigma^2(X'X)^{-1}\)
Aby wyznaczyć tę wartość należy użyć wzoru: \(\Sigma_{HA} = H\Sigma_A H'\)

xtx_inv <- solve(t(X) %*% X) # funkcja solve odwraca macierz
xtx_inv <- solve(R) %*% t(solve(R)) # alternatywny sposób wyznaczenia (X'X)^-1
sigma2 <- SSE / (n-p)

beta_var <- sigma2 * xtx_inv

Wariancja \(\hat \beta_i\)

beta_var_i = beta_var[i,i]
# lub
beta_var_i = sigma2 * diag(xtx_inv)[i]

Na ogół nie znamy rzeczywistej wariancji \(\sigma^2\), dlatego zamiast niej stosujemy estymator \(s^2 = \hat \sigma^2 = \frac{SSE}{degf}\) Korzystamy tutaj z błędu standardowego \[ SE_{hat \beta_i}^2 = s^2 (X'X)^{-1}_{ii} = s^2 (1-h_{ii})1\]

Test istotności pojedynczego współczynnika

Czyli t-statystyka oraz jej p-value dla zmiennej objaśniającej

Interpretacja hipotezy zerowej: predyktor o indeksie i nie ma istotnego wpływu na \(Y\) jeżeli w modelu są pozostałe predyktory.

Przy założeniu hipotezy zerowej poniższa statystyka ma rozkład t-Studenta z \(n-p\) stopniami swobody: statystyka t = \(\frac{\beta}{SE(\beta)}\)

p-value: prawdopodobieństwo, że statystyka przyjmie wartość bardziej sprzyjającą odrzuceniu. W przypadku statystyki t, p-value = P(T > |t|)

A p-value measures the probability of obtaining the observed results, assuming that the null hypothesis is true

Czyli istotne wzory to: \[ t = \frac{\hat \beta}{SE(\hat \beta)} \] \[ SE(\hat \beta) = \sqrt{\frac{SSE}{degf} (X'X)^{-1}}\] \[ p = 2 * min[P(T > t), P(T<t)] \]

Liczenie t-stat i p-value z definicji:

QR <- qr(X)
R <- qr.R(QR)
Q <- qr.Q(QR)
beta <- backsolve(R, t(Q) %*% Y)

xtx_inv <- solve(R)%*%t(solve(R))

SSE <- sum((Y - mod$fitted.values)^2)
df <- length(Y) - ncol(X)

sigma2 <- SSE / df

se_beta <- sqrt(sigma2 * diag(xtx_inv))

t <- beta / se_beta
p_val <- 2 * pmin(pt(t, nrow(X) - ncol(X)), 1 - pt(t, nrow(X) - ncol(X)))

Z funkcji wbudowanych:

t_builtin <- summary(mod)$coefficients[,3]

Rezydua studentyzowane

i-te rezyduum dzielimy przez pierwiastek z jego błąd standardowy: \(|r_i| = |\frac{e_i}{\sqrt{SE(e_i)}| = |\frac{e_i}{\sqrt{s^2(1-h_{ii})}}|\)

xtx_inv <- solve(t(X) %*% X)
sigma2 <- SSE/(n - p)

r <- e / sqrt(sigma2 * (1 - diag(xtx_inv)))

# funkcja wbudowana
r <- rstandard(modelinho)

Rezydua modyfikowane studentyzowane

Zdarza się tak, że duża wartość jednego rezyduum może znacząco wpłynąć na wartość \(s^2\), dlatego stosuje się rezydua modyfikowane studentyzowane. Chcąc ocenić, czy i−ta obserwacja jest odstająca, usuwamy ją ze zbioru danych

Rezyduum modyfikowane można wyznaczyć bez koniecznosci ponownego dopasowania modelu: \[ e_{i,-i} = \frac{e_i}{1-h_{ii}} \] \[ Var(e_{i,-i}) = \frac{\sigma^2}{1-h_{ii}} \] \[t_i = r_i\sqrt{\frac{n-p-1}{n-p-r_{i}^2}}\] Więcej wzorów i dokładniejsze objaśnienie tutaj

t <- rstudent(modelinho)

Studentyzowane rezyduum modyfikowane mają rozkład t-Studenta z \(n-p-1\) stopniami swobody.

Wykres indeksowy rezyduów

Według wykładu najważniejsza metoda sprawdzania dopasowania modelu, powinien być sprawdzany na początku.

e <- model$residuals
plot(1:length(e), e)
abline(0,0)

Na tym wykresie dla modelu dobrze dopasowanego powinna być widoczna chmura punktów rozrzucona losowo wokół prostej \(y=0\).

Wykres częściowej regresji

Jeżeli model jest źle wyspecyfikowany, to na wykresie częściowej regresji powinna być widoczna zależność nieliniowa między rezyduami a zmienną objaśniającą.
Wykres cześciowej regresji jest dwuwymiarowym wykresem rozproszenia \(e_{Y|X_{-j}}\) względem \(e_{x_j|x_{-j}}\).

model_y = lm(y~.-j, data=d)
model_j = lm(j~.-y, data=d) # zaleznosc j od innych zmiennych objasniajacych (nie wliczamy w to y!)

plot(model_j$residuals, model_y$residuals) # tutaj argumenty podawane są (x,y) 
abline(lm(model_y$residuals~model_j$residuals)) # tutaj (y~x)
# na osi y rezydua modelu gdzie Y objaśniany jest wszystkimi zmiennymi oprócz Xj

# weyfikacja z funkcjami wbudowanymi
car::avPlot(modelinho) # oryginalny model, y opisywane wszystkimi zmiennymi

Na tym wykresie prosta MNK ma dokładne nachylenie równe \(\beta_{j,MNK}\). Na tej podstawie również można wnioskować na temat istotności danej zmiennej.

Interpretacja: Ten wykres mówi nam, jaka jest relacja pomiedzy zmienną objasniającą j a zmienną objasnianą y, jeśli pozbędziemy się wpływu wszystkich pozostałych zmiennych objasniających.
chcemy sprawdzić, czy zmienna j wnosi cos do modelu. Trenujemy dwa modele - oba zależą od wszystkich zmiennych objaśniających oprócz j. W pierwszym objaśnianą jest oryginalny y, w drugim j. Potem bierzemy rezydua, bo rezydua to ta część zmienności, która nie jest wyjaśniana przez zmienne objaśniające. Po wyplotowaniu tego mamy wplyw zmiennej j na y przy założeniu, że wartości pozostałych zmiennych objaśniających są ustalone. Jeśli prosta jest bliska poziomej (nachylenie bliskie 0) to znaczy że zmienna j jest nieistotna

Wykres częściowych rezyduów

Na tym wykresie można lepiej zaobserwować nieliniowości względem badanej zmiennej/skupienia od wykresu częściowej regresji.

\[ y-\sum_{i\neq j}\hat\beta_i x_i = x_j \hat \beta_j + e \text{ versus }x_j\]

# funkcje wbudowane
library(faraway)
prplot(model, i) # i to numer zmiennej ktora badamy

# z definicji
plot(model$res + d$i[-outliers]*model$coefficients[i] ~ d$i[-outliers])

Na tym wykresie można też czasami zaobserwować klastry, które są grupami o różnych cechach. Tzw. paradoks Simpsona - zależność może zmieniać swój charakter w zależności od tego, czy rozważamy podzbiory danych, czy całość.

Metoda Boxa-Coxa

Ta metoda pozwala nam znaleźć taką trasformację \(h(Y)\), ze \(h(Y)\) jest generowane przez model liniowy. Przeszukujemy rodzinę parametryczną funkcji \((g_{\lambda}(.))_{\lambda \geq 0}\), zdefiniowanej jako:

\[ g_{\lambda} (y) = \frac{y^{\lambda}-1}{\lambda} \] dla\(\lambda >0\) oraz równa \(log(y)\) dla \(\lambda = 0\).

library(MASS)

boxcox1 = boxcox(modelinho, lambda = seq(0,1, length.out=100)) 
# dopasowujemy 100 modeli dla lambda # rysuje przy okazji wykres
# estymatora największej wiarogodności + przedział ufności

lambda = boxcox1$x[which(boxcox1$y == max(boxcox1$y))]

mlambda <- lm((y^lambda)~x, data=d)

Estymator Metodą Ważonych Najmniejszych Kwadratów (MWNK)

Jeżeli \(\varepsilon \sim \mathbf{N}(0, \sigma^2 W)\), gdzie \(W = diag(sigma_1^2, sigma_2^2, \dots, sigma_n^2)\) \[ Crit(b) = \sum_{i=1}^n \frac{1}{\sigma_i^2} |Y_i - x_ib|^2\] Sens ważenia jest taki, że obserwacje, które mają większą wariancje powinny mieć mniejszy wpływ na estymator \(\hat \beta\). \[ \hat{\beta}_{MWNK} = (X'WX)^{-1}X'WY \]

rez <- residuals(mod)

W <- 1 / lm(abs(rez) ~ mod$fitted.values)$fitted.values^2 # wagi

# wizualizacja żeby sprawdzić jak wygląda rozkład rezyduów
plot(mod$fitted.values, abs(rez))

# nowy model, pozbyliśmy się heteroskedastyczności
mod_weighted <- lm(y ~ x, weights = w)

# możemy teraz zobaczyć, że nie ma już heteroskedastyczności
plot(1:length(rez), rstudent(mod_weighted))

Metoda Uogólnionych Najmniejszych Kwadratów (MUNK)

W tym przypadku rozważamy ogólny przypadek macierzy kowariancji błędów, redukując ten problem do takiego, który już znamy i potrafimy rozwiązać przy pomocy MNK.

Jeżeli \(\Sigma_{\varepsilon} = \sigma^2 \Sigma\) to \(\Sigma^{-1/2} \varepsilon \sim \mathbf{N}(0, \sigma^2 I)\) Wobec tego jeżeli weźmiemy \(Y_1 = \Sigma^{-1/2} Y\) oraz \(X_1 = \Sigma^{-1/2} X\) to otrzymamy model liniowy z błędami o stałej wariancji, który możemy rozwiązać za pomocą MNK.

\(\beta_{MUNK} = (X'\Sigma^{-1}X) ^{-1} X'\Sigma^{-1}Y\)

Łatwo sprawdzić, że \(\Sigma_{\hat \beta_{MUNK}} = \sigma^2 (X'\Sigma^{-1}X)^{-1}\)

Rezydua studentyzowane lub modyfikowane studentyzowane lub modyfikowane

Można np. przeanalizować wykres indeksowy rezyduów. Wykres kwantylowy też pozwala wykryć nietypowe rezydua.

Wzory podane w zakładce Rezydua

Heurystyczna metoda \(r_i \geq 2\) lub \(t_i \geq 2\)

res <- rstudent(modelinho) # rezydua studentyzowane modyfikowane
res2 <- rstandard(modelinho) # rezydua studentyzowane
outliers = which(abs(res) >= 2) # tutaj dostajemy obserwacje odstające
model_res = lm(Y~., data=d, subset = -outliers)

Rezydua modyfikowane mają rozkład tstudenta z \(n-p-1\) stopniami swobody. Wobec tego można zastosować test tstudenta, aby sprawdzić czy rezydua są odstające. Uwaga! jeżeli sprawdzamy czy rezyduum o najwiekszej wartosci \(t\) jest odstające na zadanym poziomie istotności \(\alpha\), to powinniśmy sprawdzić czy \(t > t_{n-p-1, 1-\alpha/n}\). Jest to problem wielokrotnych powtórzeń. Tę operację nazywamy poprawką Bonferroniego.

Jeżeli wybieramy, np. pierwszą obserwację (wybor obserwacji jest niezalezny od danych) to możemy zastosować zwykły test tstudenta.

REZYDUA MODYFIKOWANE NIESTUDENTYZOWANE

X<- model.matrix(mod)
res <- mod$residuals
res_mod <- res/(1-hatvalues(mod))

odległość Cooka

## Odległość Cooka
Oceniamy wpływ usunięcia i -tej obserwacji na prognozę wszystkich obserwacji.
Obserwacja jest wpływowa jeżeli jej usunięcie znacząco zmienia parametry modelu.

$$D_i = \frac{\|X b_{-i} - \hat{Y}\|}{ps^2}$$
$$D_i = \frac{r_i^2 h_{ii}}{p(1-h_{ii})}$$

```{}
# funkcja wbudowana
cooks.distance(modelinho) 
plot(modelinho, 4) # wykres cooka

bigcook <- which(cooks.distance(modelinho)>=1) # jedna z heurystyk

mod2 <- lm(Savings ~ .-Country, data = s, subset=-bigcook)