boissons<-read.table("/Users/carraro/Documents/Laurent/enseignement/axe_MSA/06-07/boissons.txt",dec=".",header=TRUE)
mod<-lm(temps~nb+distance,data=boissons)
summary(mod)

#résidus bruts
layout(matrix(1:4,2,2))
plot(mod$residuals, xlab="index",ylab="résidus")
plot(mod$residuals~mod$fitted.values, xlab="réponse estimée",ylab="résidus")
plot(mod$residuals~boissons$nb, xlab="nombre",ylab="résidus")
plot(mod$residuals~boissons$dist, xlab="distance",ylab="résidus")

#résidus standardisés
std<-rstandard(mod)
layout(matrix(1:4,2,2))
plot(std, xlab="index",ylab="résidus standardisés")
plot(std~mod$fitted.values, xlab="réponse estimée",ylab="résidus standardisés")
plot(std~boissons$nb, xlab="nombre",ylab="résidus standardisés")
plot(std~boissons$dist, xlab="distance",ylab="résidus standardisés")

#résidus studentisés
stud<-rstudent(mod)
layout(matrix(1:4,2,2))
plot(stud, xlab="index",ylab="résidus studentisés")
q<-qt(0.975,df=24)
abline(a=q,b=0,lty="dashed")
abline(a=-q,b=0,lty="dashed")
plot(stud~mod$fitted.values, xlab="réponse estimée",ylab="résidus studentisés")
abline(a=q,b=0,lty="dashed")
abline(a=-q,b=0,lty="dashed")
plot(stud~boissons$nb, xlab="nombre",ylab="résidus studentisés")
abline(a=q,b=0,lty="dashed")
abline(a=-q,b=0,lty="dashed")
plot(stud~boissons$dist, xlab="distance",ylab="résidus studentisés")
abline(a=q,b=0,lty="dashed")
abline(a=-q,b=0,lty="dashed")

############################
#on enlève l'observation n°9
boissons2<-boissons[c(1:8,10:25),]
mod2<-lm(temps~nb+distance,data=boissons2)
summary(mod2)

#résidus studentisés
stud2<-rstudent(mod2)
layout(matrix(1:4,2,2))
plot(stud2, xlab="index",ylab="résidus studentisés")
q<-qt(0.975,df=23)
abline(a=q,b=0,lty="dashed")
abline(a=-q,b=0,lty="dashed")
plot(stud2~mod2$fitted.values, xlab="réponse estimée",ylab="résidus studentisés")
abline(a=q,b=0,lty="dashed")
abline(a=-q,b=0,lty="dashed")
plot(stud2~boissons2$nb, xlab="nombre",ylab="résidus studentisés")
abline(a=q,b=0,lty="dashed")
abline(a=-q,b=0,lty="dashed")
plot(stud2~boissons2$dist, xlab="distance",ylab="résidus studentisés")
abline(a=q,b=0,lty="dashed")
abline(a=-q,b=0,lty="dashed")


layout(matrix(1,1,1))
qqnorm(stud2,ylab="résidus studentisés")
qqline(stud2)

#validation
ks.test(stud, "pt", 24)
ks.test(stud2, "pt", 23)
ks.test(stud,"pnorm",0,1)
ks.test(stud2,"pnorm",0,1)
shapiro.test(stud)
shapiro.test(stud2)


#leviers
lev<-hatvalues(mod)
plot(lev,type="h",xlab="index",ylab="levier",ylim=c(0,1))
points(lev,pch=19,col="blue")

#distance de Cook
cook<-cooks.distance(mod)
plot(cook,type="h",xlab="index",ylab="distance de Cook")
points(cook,pch=19,col="blue")

#########################################
#intervalles de confiance et de prévision
data_pluie<-read.table("/Users/carraro/Documents/Laurent/enseignement/axe_MSA/06-07/donnees_pluie.csv",sep=";",header=TRUE, dec=",")
pluie1<-lm(rendement~pluie,data=data_pluie)
pluie2<-lm(rendement~pluie+I(pluie^2),data=data_pluie)
newdata<-data.frame(pluie=seq(0,0.5,length=500))

#modèle de degré 1
p.conf<-predict(pluie1,newdata=newdata,interval="confidence") 
p.pred<-predict(pluie1,newdata=newdata,interval="prediction") 
plot(data_pluie,main="intervalles de confiance et de prévision pour le modèle de degré 1")
lines(newdata$pluie,p.conf[,2],lty=2,col=2)
lines(newdata$pluie,p.conf[,3],lty=2,col=2)
lines(newdata$pluie,p.conf[,1],col=1)
lines(newdata$pluie,p.pred[,2],lty=3,col=3)
lines(newdata$pluie,p.pred[,3],lty=3,col=3)
legend("topleft",legend=c("intervalle de confiance","intervalle de prévision"),lty=c(2,3),col=c(2,3))

#modèle de degré 2
p.conf<-predict(pluie2,newdata=newdata,interval="confidence") 
p.pred<-predict(pluie2,newdata=newdata,interval="prediction") 
plot(data_pluie,main="intervalles de confiance et de prévision pour le modèle de degré 2")
lines(newdata$pluie,p.conf[,2],lty=2,col=2)
lines(newdata$pluie,p.conf[,3],lty=2,col=2)
lines(newdata$pluie,p.conf[,1],col=1)
lines(newdata$pluie,p.pred[,2],lty=3,col=3)
lines(newdata$pluie,p.pred[,3],lty=3,col=3)
legend("topleft",legend=c("intervalle de confiance","intervalle de prévision"),lty=c(2,3),col=c(2,3))

#modèle de degré 7
pluie7<-lm(rendement~poly(pluie,degree=7),data=data_pluie)
p.conf<-predict(pluie7,newdata=newdata,interval="confidence") 
p.pred<-predict(pluie7,newdata=newdata,interval="prediction") 
plot(data_pluie,ylim=c(0,1),main="intervalles de confiance et de prévision pour le modèle de degré 7")
lines(newdata$pluie,p.conf[,2],lty=2,col=2)
lines(newdata$pluie,p.conf[,3],lty=2,col=2)
lines(newdata$pluie,p.conf[,1],col=1)
lines(newdata$pluie,p.pred[,2],lty=3,col=3)
lines(newdata$pluie,p.pred[,3],lty=3,col=3)
legend("topleft",legend=c("intervalle de confiance","intervalle de prévision"),lty=c(2,3),col=c(2,3))