rm(list=ls()) #clear the workspace
set.seed(2020) #so that everyone gets the same results


#example situations:

#1. Ordinary observations
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)) #store the regression coefficients
r1 

#2. an outlier (but not an influential observation)
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)) #store the regression coefficients
r2 

summary(lm(y2~x2)) #regression results with all observations
summary(lm(y2[-125]~x2[-125])) #regression results with the outlier removed

#3. an observation that is both an outlier and influential
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)) #store the regression coefficients
r3 

summary(lm(y3~x3)) #regression results with all observations
summary(lm(y3[-125]~x3[-125])) #regression results with the outlier removed

#4. an influential observation that is not an outlier 
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)) #store the regression coefficients
r4 

Reg<-rbind(r1[2],r2[2],r3[2],r4[2])
rownames(Reg)=c("ordinary", "outlier", "outlier & influential", "influential only")
Reg

#Analysing a sample
#we start with the libraries
#note: require() alone does not load the package after installing it - hence library() in the braces
if(!require("car")) {install.packages("car"); library("car")} #example data and plots from the car package
if(!require("olsrr")) {install.packages("olsrr"); library("olsrr")} #outlier detection measures



fit <- lm(mpg~disp+hp+wt, data=mtcars) #what drives fuel consumption - disp = displacement; hp = horsepower; wt = weight; 
summary(fit)

#Quantile-quantile plot: we compare the distribution of the residuals with the normal distribution
qqnorm(fit$residuals)
qqline(fit$residuals, col="red")
#Alternative version, with confidence bands based on the t-distribution
qqPlot(fit, main="QQ Plot")

#Leverage plots
leveragePlots(fit) 


#Bonferroni test - a formal test for outliers
#based on studentized deleted residuals; the p-value is adjusted for the number of tests (N)
outlierTest(fit)


#Studentized Residual Plot
ols_plot_resid_stud(fit)

#Studentized Residuals vs Leverage Plot
ols_plot_resid_lev(fit)


#Cook's distance (Cook's D)
ols_plot_cooksd_bar(fit)
ols_plot_cooksd_chart(fit)
#Cook's D values
cooks.distance(fit)

#DFFITS
ols_plot_dffits(fit)
dffits(fit)

#DFBETAS (the standardised version - this is what ols_plot_dfbetas draws)
ols_plot_dfbetas(fit)
dfbetas(fit)

#The same once again, using the plot function
plot(fit) 
#Press enter and you get the following plots in turn:
# 1 Residuals vs Fitted
# 2 Q-Q plot
# 3 Scale - location
# 4 Residuals vs Leverage


#Outliers or a misspecified model?
#all three largest studentized residuals are POSITIVE (the model under-predicts)
sort(rstudent(fit), decreasing=TRUE)[1:3]
#and they sit at the extremes of the weight distribution - a symptom of a convex relationship, not of bad data
mtcars[c("Toyota Corolla","Fiat 128","Chrysler Imperial"), c("mpg","wt")]
range(mtcars$wt)
#check: log of weight instead of its level
fit_log <- lm(mpg~disp+hp+log(wt), data=mtcars)
c(linear=summary(fit)$r.squared, log=summary(fit_log)$r.squared)
outlierTest(fit_log)


#Estimation that accounts for outliers
if(!require("robustbase")) {install.packages("robustbase"); library("robustbase")} #package containing LTS

# 1. Least Trimmed Squares (LTS)
fit2<-ltsReg(mpg~disp+hp+wt, data=mtcars, alpha=0.9) #alpha - what % of observations is to remain in the sample
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
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 - the quantiles for which we estimate the regression
summary(fit4)
plot(summary(fit4))

#Summary
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 #compare tau=0.5 with LAD



#####################################
# Homework
#For each of the x and y pairs (number 2, 3, 4) draw a qqplot and compute Cook's distance, DFFITS and DFBETAS:
#run a quantile regression (quantiles 10, 20, ..., 90) for the pair x3 and y3, 
#then store the coefficient values and compare them with the OLS results