library (car)
library (nlme)		# gls
library (lmtest)	#	D-W test
library (ggm)

summary (mtcars)

alpha <- 0.05
n <- dim(mtcars)[1]  # pocet mereni
k <- dim(mtcars)[2]-1  # pocet promennych 
mod <- lm(mpg ~ ., data = mtcars) # klasicky model
prehled <- summary(mod)
prehled
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
y<-mtcars$mpg
abs(cor(y,X))   # 5. promenna ma nejvetsi korelaci se spotrebou
mod1 <- lm(mpg~wt, data = mtcars)
prehled1 <- summary(mod1)
prehled1

ID1 <- prehled1$r.squared
ID1/(1-ID1)*(n-2)
qf(1-alpha,1,n-2) # F1 je vetsi => "wt" ponechame v modelu

library(ggm)
abs(pcor(c("mpg","cyl","wt"),var(mtcars)))
abs(pcor(c("mpg","disp","wt"),var(mtcars)))
abs(pcor(c("mpg","hp","wt"),var(mtcars)))
abs(pcor(c("mpg","drat","wt"),var(mtcars)))
abs(pcor(c("mpg","qsec","wt"),var(mtcars)))
abs(pcor(c("mpg","vs","wt"),var(mtcars)))
abs(pcor(c("mpg","am","wt"),var(mtcars)))
abs(pcor(c("mpg","gear","wt"),var(mtcars)))
abs(pcor(c("mpg","carb","wt"),var(mtcars)))
# promenna "cyl" vyhrala => zahrneme ji do modelu
mod2 <- update(mod1,.~. + cyl)
prehled2 <- summary(mod2)
prehled2

ID2 <- prehled2$r.squared
(ID2-ID1)/(1-ID2)*(n-3)
qf(1-alpha,2,n-3) # F2 je jeste porad vetsi => "cyl" ponechame v modelu

abs(pcor(c("mpg","hp","wt","cyl"),var(mtcars)))
abs(pcor(c("mpg","disp","wt","cyl"),var(mtcars)))
abs(pcor(c("mpg","drat","wt","cyl"),var(mtcars)))
abs(pcor(c("mpg","qsec","wt","cyl"),var(mtcars)))
abs(pcor(c("mpg","vs","wt","cyl"),var(mtcars)))
abs(pcor(c("mpg","am","wt","cyl"),var(mtcars)))
abs(pcor(c("mpg","gear","wt","cyl"),var(mtcars)))
abs(pcor(c("mpg","carb","wt","cyl"),var(mtcars)))
# promenna "hp" tesne vyhrala => zahrneme ji do modelu

mod3 <- update(mod2,.~. + hp)
prehled3 <- summary(mod3)
prehled3
ID3 <- prehled3$r.squared
(ID3-ID2)/(1-ID3)*(n-4)
qf(1-alpha,3,n-4) # F3 je mensi => "hp" do modelu nezahrneme

# NEBO
step(mod) #!!! uplne jiny vysledek
#
# Vitez
mod <- mod2 # vitezny model
prehled <- summary(mod)
prehled
# overeni multikolinearity
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)
# Multikolinearita jiz neni prokazana

# Testovani normality
# graficky
res <- mod2$residuals
qqnorm(res)
qqline(res)
# KS test
ks.test(res,"pnorm")

# Shapiro-Wilk
shapiro.test(res)
# Testy sice normalitu zamitly, ale celkem tesne. I z obrazku je patrne, ze muzeme residua povazovat za normalni.