Data <- read.table("ozone_01.txt",header=TRUE);
 # Data <- read.table("essence_bon.dat",header=TRUE);

 # Fit the regression model with complete variables
 mod.lm <- lm(maxO3~ T9 + T12 + T15+ Ne9 +Ne12 +Ne15 + Vx9 + Vx12 + Vx15 ,data= Data); 
 #modr.lm <- lm(Consommation ~  Permis+ Revenu+ Route+ Population ,data= Data); 
 #modc.lm <- lm(Consommation ~   Population + Permis+ Revenu+ Route+ Permis_capita+ Consom_capita ,data= Data); 



############################################
# Diagnostic of multicollinearity      #
########################################################
# The 'Paires' make scatter plots between each two variables.

#pairs(Consommation ~  Permis+ Revenu+ Route+ Population ,data= Data)
# Get the multicolinrearatiy 1 pic
#pairs(maxO3~ T9 + T12 + T15+ Ne9 +Ne12 +Ne15 + Vx9 ,data= Data)
# For example
pairs(maxO3~ T9 + T12 + Vx9 ,data= Data)


# Inflation factor of the variance

library(car) # The function VIF is in the Car package

# Which one is greater than 10
# If VIF of the jth variable is greater than 10, 
# it shows that a colinearity between jth and some of others
vif(mod.lm)# See the influation factor of the variance

#################################
# Check postulates( Hypothesis) #
#################################

modele2.lm<-  lm(maxO3~ T9 + T12 + Ne9 +Ne12 + Vx9 + Vx12 + Vx15 ,data= Data);
ri<-rstudent(modele2.lm) # To caculate the Rstudent
plot(predict(modele2.lm),ri)

objects(modele2.lm)
modele2.lm$rank

# The verification of the normality graphs

par(mfrow=c(1,3))
hist(ri)
boxplot(ri)
title('Moustache diagram of residuals')
qqnorm(ri)

#Test of normality

shapiro.test(ri)   #H0: W=1, We accept the normality distribution

ks.test(ri,'pnorm')    #H0 D=0, We accept the normality distribution


#To verify the with sesiduals graph


leveragePlots(modele2.lm)

# To verify the presence of the influentes / aberrantes values

par(mfrow=c(1,1))
influencePlot(modele2.lm,scale=1,id.n = 1)

dcook<-cooks.distance(modele2.lm)

plot(dcook)
 
4/(length(Data[,1])-modele2.lm$rank)

a=which(dcook>4/(length(Data[,1])-modele2.lm$rank))
dcook[a]

DFBETA<-dfbeta(modele2.lm)#>2/sqrt(n)
DFBETA
2/sqrt(length(Data[,1]))
a1=which(DFBETA[,2]>2/length(Data[,1]))
DFBETA[a1,2]

DFFITS<-dffits(modele2.lm)#Not belang to [-2,2]
DFFITS

hii<-hatvalues(modele2.lm)
hii

COVRATIO<-covratio(modele2.lm)
COVRATIO