#####################################
#regression - option maths applis
#EMSE
#L. Carraro (Telecom Saint-Etienne)
#decembre 09
#####################################

load("/Users/carraro/Documents/Laurent/enseignement/axe_MSA/06-07/cartons/cartons.RData")
dechetsA<-dechets[dechets$prestataire=="A",]
dechetsB<-dechets[dechets$prestataire=="B",]

#interactions
coplot(log10(cartons)~log10(CA)|log10(surface),panel=panel.smooth,data=dechets)
coplot(log10(cartons)~log10(CA)|categorie,panel=panel.smooth,data=dechets)
coplot(log10(cartons)~log10(CA)|prestataire,panel=panel.smooth,data=dechets)

#données pour prestataire B
coplot(log10(cartons)~log10(CA)|log10(surface),panel=panel.smooth,data=dechetsB)
coplot(log10(cartons)~log10(CA)|categorie,panel=panel.smooth,data=dechetsB)
coplot(log10(cartons)~log10(CA)|categorie+log10(surface),panel=panel.smooth,data=dechetsB)

mod_lB1<-lm(log10(cartons)~log10(CA)+log10(surface)+categorie,data=dechetsB)
summary(mod_lB1)
mod_lB2<-lm(log10(cartons)~log10(CA),data=dechetsB)
summary(mod_lB2)

#sélections automatiques de modèles
mod_l0<-lm(log10(cartons)~1,data=dechets)
#forward
addterm(mod_l0,test="F",scope=cartons~(log10(CA)+log10(surface)+prestataire+categorie)^2,sorted=TRUE)
mod_l1<-update(mod_l0,.~.+log10(CA))
addterm(mod_l1,test="F",scope=cartons~(log10(CA)+log10(surface)+prestataire+categorie)^2,sorted=TRUE)
mod_l2<-update(mod_l1,.~.+prestataire)
addterm(mod_l2,test="F",scope=cartons~(log10(CA)+log10(surface)+prestataire+categorie)^2,sorted=TRUE)
mod_l3<-update(mod_l2,.~.+log10(CA):prestataire)
addterm(mod_l3,test="F",scope=cartons~(log10(CA)+log10(surface)+prestataire+categorie)^2,sorted=TRUE)
#backward
#départ = modèle linéaire sans interaction
mod_lc<-lm(log10(cartons)~log10(CA)+log10(surface)+prestataire+categorie,data=dechets)
dropterm(mod_lc,test="F",sorted=TRUE)
mod_lc1<-update(mod_lc,.~.-categorie)
dropterm(mod_lc1,test="F",sorted=TRUE)
mod_lc2<-update(mod_lc1,.~.-log10(surface))
dropterm(mod_lc2,test="F",sorted=TRUE)
#départ = modèle linéaire avec interactions
mod_lc<-lm(log10(cartons)~(log10(CA)+log10(surface)+prestataire+categorie)^2,data=dechets)
dropterm(mod_lc,test="F",sorted=TRUE)
mod_lc1<-update(mod_lc,.~.-log10(surface):prestataire)
dropterm(mod_lc1,test="F",sorted=TRUE)
mod_lc2<-update(mod_lc1,.~.-log10(CA):prestataire)
dropterm(mod_lc2,test="F",sorted=TRUE)
mod_lc3<-update(mod_lc2,.~.-log10(CA):categorie)
dropterm(mod_lc3,test="F",sorted=TRUE)
#avec AIC et BIC
#AIC
st<-step(mod_l0, ,direction = "both",scope=list(lower=cartons~1,upper=cartons~(log10(CA)+log10(surface)+prestataire+categorie)^2),trace = 1)
#BIC
step(mod_l0, ,direction = "both",k=log(137),scope=list(lower=cartons~1,upper=cartons~(log10(CA)+log10(surface)+prestataire+categorie)^2),trace = 1)

#étude des données de pluie : exemple 2
data_pluie<-read.table("/Users/carraro/Documents/Laurent/enseignement/axe_MSA/06-07/donnees_pluie.csv",sep=";",header=TRUE, dec=",")
plot(data_pluie,xlab="pluie en m",ylab="rendement sans unité")
pluie1<-lm(rendement~pluie,data=data_pluie)
summary(pluie1)
pluie2<-lm(rendement~pluie+I(pluie^2),data=data_pluie)
summary(pluie2)
pluie3<-lm(rendement~pluie+I(pluie^2)+I(pluie^3),data=data_pluie)
summary(pluie3)

#résidus pour les cartons
#prestataire B
stud<-studres(mod_lB2)
plot(stud,main="résidus studentisés",sub="modèle linéaire log10(cartons)~log10(CA) - prestataire B",ylab="résidus")

layout(matrix(1:3,3,1))
plot(stud~log10(CA),data=dechetsB,main="résidus studentisés",ylab="résidus")
plot(stud~log10(surface),data=dechetsB,main="résidus studentisés",ylab="résidus")
boxplot(stud~categorie,data=dechetsB,horizontal=TRUE,main="résidus studentisés")

layout(matrix(1,1,1))
qqnorm(stud)
qqline(stud)

layout(matrix(1:2,1,2))
qqnorm(stud[log10(dechetsB$CA)<4.4],main=c("Droite de Henri","résidus tels que log10(CA) < 4.4)"))
qqline(stud[log10(dechetsB$CA)<4.4])
qqnorm(stud[log10(dechetsB$CA)>4.4],main=c("Droite de Henri","résidus tels que log10(CA) > 4.4)"))
qqline(stud[log10(dechetsB$CA)>4.4])

#on enlève l'observation n°12 dont le résidu dépasse 4
dechetsB2<-dechetsB[c(1:11,13:35),]
mod_lB3<-lm(log10(cartons)~log10(CA),data=dechetsB2)
summary(mod_lB3)

stud2<-studres(mod_lB3)
plot(stud2,main=c("résidus studentisés - prestataire B sans observation n°12","modèle linéaire log10(cartons)~log10(CA)"),ylab="résidus")

layout(matrix(1,1,1))
qqnorm(stud2)
qqline(stud2)

layout(matrix(1:2,1,2))
qqnorm(stud2[log10(dechetsB2$CA)<4.4],main=c("Droite de Henri","résidus tels que log10(CA) < 4.4)"))
qqline(stud2[log10(dechetsB2$CA)<4.4])
qqnorm(stud2[log10(dechetsB2$CA)>4.4],main=c("Droite de Henri","résidus tels que log10(CA) > 4.4)"))
qqline(stud2[log10(dechetsB2$CA)>4.4])

#tests d'adéquation
ks.test(stud, "pt",df=33)
ks.test(stud, "pnorm",0,1)
shapiro.test(stud)
ks.test(stud2, "pt",df=32)
ks.test(stud2, "pnorm",0,1)
shapiro.test(stud2)