rm(list=ls()) #czyszczenie obszaru roboczego
set.seed(2020) #aby wszyscy otrzymali takie same wyniki


#przykładowe sytuacje:

#1. Zwyczajne obserwacje
n=250
x=runif(n)
y=1+2*x+0.1*rnorm(n)
x[125]=0.5
y[125]=1+2*0.5+0.1


plot(y~x, col="grey", pch=16, xlim=c(0,1),ylim=c(0.5,5))
abline(lm(y~x), col="blue")
points(y[125]~x[125], col="red", pch=16)
r1<-coef(lm(y~x)) #zapisujemy współczynniki regresji
r1 

#2. obserwacja odstająca (ale nie wpływowa)
x2<-x
y2<-y
y2[125]=y[125]+2

plot(y2~x2, col="grey", pch=16, xlim=c(0,1),ylim=c(0.5,5))
abline(lm(y2~x2), col="blue")
points(y2[125]~x2[125], col="red", pch=16, cex=1.3)
lm(y2~x2)
r2<-coef(lm(y2~x2)) #zapisujemy współczynniki regresji
r2 

summary(lm(y2~x2)) #wyniki regresji ze wszystkimi obserwacjami
summary(lm(y2[-125]~x2[-125])) #wyniki regresji z usuniętą obserwacją odstającą

#3. obserwacja odstająca i wpływowa
x3<-x
y3<-y
x3[125]=3
y3[125]=0.2

plot(y3~x3, col="grey", pch=16, xlim=c(0,3),ylim=c(0,7))
abline(lm(y3~x3), col="blue")
points(y3[125]~x3[125], col="red", pch=16, cex=1.3)
lm(y3~x3)
r3<-coef(lm(y3~x3)) #zapisujemy współczynniki regresji
r3 

summary(lm(y3~x3)) #wyniki regresji ze wszystkimi obserwacjami
summary(lm(y3[-125]~x3[-125])) #wyniki regresji z usuniętą obserwacją odstającą

#4. obserwacja wpływowa, ale nie odstająca 
x4<-x
y4<-y
x4[125]=3
y4[125]=1+2*3+0.1

plot(y4~x4, col="grey", pch=16, xlim=c(0,3),ylim=c(0,7))
abline(lm(y4~x4), col="blue")
points(y4[125]~x4[125], col="red", pch=16)
lm(y4~x4)
r4<-coef(lm(y4~x4)) #zapisujemy współczynniki regresji
r4 

Reg<-rbind(r1[2],r2[2],r3[2],r4[2])
rownames(Reg)=c("zwykła", "odstająca", "odst. i wpływowa", "tylko wpływowa")
Reg

#Badanie próby
#zaczynamy od bibliotek
if(!require("car")) {install.packages("car"); library("car")} #przykładowe dane i wykresy z pakietu car
if(!require("olsrr")) {install.packages("olsrr"); library("olsrr")} #miary wykrywające obserwacje odstające



fit <- lm(mpg~disp+hp+wt, data=mtcars) #od czego zależy zużycie paliwa - disp = obj. skokowa; hp = moc; wt = masa; 
summary(fit)

#Wykres kwantyl-kwantyl: porównujemy rozkład reszt z rozkładem normalnym
qqnorm(fit$residuals)
qqline(fit$residuals, col="red")
#Wersja alternatywna, gdzie mamy wyznaczone przedziały ufności na podstawie rozkładu t-studenta
qqPlot(fit, main="QQ Plot")

#Wykresy dźwigni
leveragePlots(fit) 


#Test Bonferroniego - formalny test obserwacji odstających
#oparty na resztach studentyzowanych usuniętych; p-value skorygowane o liczbę testów (N)
outlierTest(fit)


#Studentized Residual Plot
ols_plot_resid_stud(fit)

#Studentized Residuals vs Leverage Plot
ols_plot_resid_lev(fit)


#dystans Cook'a (D-Cook's distance)
ols_plot_cooksd_bar(fit)
ols_plot_cooksd_chart(fit)
#wartość D-Cooka
cooks.distance(fit)

#DFFITS
ols_plot_dffits(fit)
dffits(fit)

#DFBETAS (wersja standaryzowana - taka jest rysowana przez ols_plot_dfbetas)
ols_plot_dfbetas(fit)
dfbetas(fit)

#Jeszcze raz to samo za pomocą funkcji plot
plot(fit) 
#Klikamy enter i otrzymujemy po kolei następujace wykresy:
# 1 Residuals vs Fitted
# 2 Q-Q plot
# 3 Scale - location
# 4 Residuals vs Leverage


#Obserwacje odstające czy błędna specyfikacja?
#wszystkie trzy największe reszty studentyzowane są DODATNIE (model niedoszacowuje)
sort(rstudent(fit), decreasing=TRUE)[1:3]
#i leżą na krańcach rozkładu masy - to objaw zależności wypukłej, a nie błędnych danych
mtcars[c("Toyota Corolla","Fiat 128","Chrysler Imperial"), c("mpg","wt")]
range(mtcars$wt)
#sprawdzenie: logarytm masy zamiast poziomu
fit_log <- lm(mpg~disp+hp+log(wt), data=mtcars)
c(liniowy=summary(fit)$r.squared, logarytm=summary(fit_log)$r.squared)
outlierTest(fit_log)


#Estymacja uwzględniająca występowanie outlierów
if(!require("robustbase")) {install.packages("robustbase"); library("robustbase")} #pakiet zawierający LTS

# 1. Least Trimmed Squares (LTS)
fit2<-ltsReg(mpg~disp+hp+wt, data=mtcars, alpha=0.9) #alpha - jaki % obserwacji ma pozostać w próbie
summary(fit2)

#2. Least Absolute Deviations (LAD) estimator
if(!require("L1pack")) {install.packages("L1pack"); library("L1pack")}
fit3<-lad(mpg~disp+hp+wt, data=mtcars)
summary(fit3)

#3. Quantile regression (regresja kwantylowa)
if(!require("quantreg")) {install.packages("quantreg"); library("quantreg")}
fit4<-rq(mpg~disp+hp+wt, data=mtcars, tau = seq(0.1, 0.9, by=0.1)) #tau - dla których kwantyli szacujemy regresję
summary(fit4)
plot(summary(fit4))

#Podsumowanie
Reg_coef<-cbind(coef(fit),coef(fit2),coef(fit3))
colnames(Reg_coef)<-c("OLS","LTS","LAD")
Reg_coef<-cbind(Reg_coef, coef(fit4))
Reg_coef #porównaj tau=0.5 z LAD



#####################################
# Praca domowa
#Dla każdej z par x i y (numer 2, 3, 4) narysuj wykres qqplot oraz policz wartości odległości Cooka, DFFITS i DFBETAS:
#przeprowadź regresję kwantylową (kwantyle 10, 20, ..., 90) dla pary x3 i y3, 
#a następnie zapisz wartości współczynników i porównaj je z wynikami KMNK