# Multikolinearita
# priklad
rm(list=ls())
load("vydaje.RData")
attach(domacnosti)
domacnosti
alpha <- 0.05
n <- dim(domacnosti)[1]  # pocet mereni
k <- dim(domacnosti)[2]-1  # pocet promennych 
mod <- lm(vydaje ~ ., data = domacnosti) # klasicky model
(prehled <- summary(mod))
X <- model.matrix(mod)[,-1]  # matice s vysvetlujicimi promennymi
det(solve(t(X)%*%X))          # determinant (X'X)^(-1) je malinky => multikolinearita
library(car)
vif(mod)                      # VIF - variance inflarion factors
R <- cor(X)
R
W <- -(n-1-1/6*(2*k+7))*log(det(R)) # test H_0: R=I
W
qchisq(1-alpha,k*(k-1)/2)

d <- diag(solve(R))    # stejne jako VIF
Fj <- (n-k)/(k-1)*(d-1)
Fj
qf(1-alpha,k-1,n-k)

# Odstraneni multikolinearity
abs(cor(vydaje,X))   # 1. promenna ma nejvetsi korelaci s vydaji
mod1 <- lm(vydaje~clenu, data = domacnosti)
prehled1 <- summary(mod1)
prehled1
yhat1 <- mod1$fitted.values
SR1 <- sum((yhat1-mean(vydaje))^2)
SE1 <- sum(mod1$residuals^2)
(F1 <- (n-2)*SR1/SE1)  # stejne jako F-statistika
# nebo
ID1 <- prehled1$r.squared
ID1/(1-ID1)*(n-2)
qf(1-alpha,1,n-2) # F1 je vetsi => "clenu" ponechame v modelu

library(ggm)
#abs(pcor(c(1,3,2),var(domacnosti)))
abs(pcor(c("vydaje","deti","clenu"),var(domacnosti)))
abs(pcor(c("vydaje","vek","clenu"),var(domacnosti)))
abs(pcor(c("vydaje","prijem","clenu"),var(domacnosti)))
# promenna vek vyhrala => zahrneme ji do modelu
mod2 <- lm(vydaje~cbind(clenu,vek), data = domacnosti)
#mod2 <- update(mod1,.~.+vek,data = domacnosti)
prehled2 <- summary(mod2)
prehled2
yhat2 <- mod2$fitted.values
SR2 <- sum((yhat2-mean(vydaje))^2)
SE2 <- sum(mod2$residuals^2)
(F2 <- (n-3)*(SR2-SR1)/(SE2))
# nebo
ID2 <- prehled2$r.squared
(ID2-ID1)/(1-ID2)*(n-3)
qf(1-alpha,2,n-3) # F2 je jeste porad vetsi => "vek" ponechame v modelu

abs(pcor(c("vydaje","deti","clenu","vek"),var(domacnosti)))
abs(pcor(c("vydaje","prijem","clenu","vek"),var(domacnosti)))
# promenna deti vyhrala => zahrneme ji do modelu

mod3 <- lm(vydaje~cbind(clenu,vek,deti), data = domacnosti)
prehled3 <- summary(mod3)
prehled3
yhat3 <- mod3$fitted.values
SR3 <- sum((yhat3-mean(vydaje))^2)
SE3 <- sum(mod3$residuals^2)
(F3 <- (n-4)*(SR3-SR2)/(SE3))
# nebo
ID3 <- prehled3$r.squared
(ID3-ID2)/(1-ID3)*(n-4)
qf(1-alpha,3,n-4) # F3 uz je mensi => "deti" do modelu zahrnovat nebudeme

# NEBO
step(mod)

#
detach(domacnosti)
# ---------------------------------------------------------------------------