progR

#D. NERINI, Nov. 2009
#Comparaison prevision avec Arbre de classification et analyse discriminante
#sur donnees etang de Berre disponible dans la rubrique precedente 'data mining"
#oxy : % oxygene dissous, a predire
#dens : densite de l eau
#moyve : moyenne du vt (m/s)
#maxve : maximum vent (m/s)
#deb : debit centrale (m3/s/j)
#debv : debit de la veille
#deb7 : debit cumule (m3/s/semaine)
#deb15 : debit cumule (m3/s/quinzaine)
##Ne pas hesiter a regarder l aide de R lorsque la fonction utilisee
#n est pas explicite
##Bien lire le script avant de le lancer pour le faire fonctionner par etape

#lecture des données apres stockage ou vous voulez
x=read.table("oxy.dat",header=T,row.names=1)

###Enlevement des données manquantes marquees par NA
ind=which(is.na(x[,1]))
x=x[-ind,]

library(tree)
#library(rpart) aussi a utiliser
#arbre de regression

#on fabrique une formule pour eviter de l ecrire a chaque fois
form0 = as.formula(paste("oxy~ ", "dens+moyve+maxve+deb+debv+deb7+deb15"))

#arbre de regression
xt=tree(form0,data=x)
plot(xt)
text(xt)

#prevision avec echantillon de construction + representation
xpred=predict(xt)
plot(x[,1],xpred,xlab="observed",ylab="predited")
abline(1,1,lty=2)

#analyse discriminante
library(MASS)
#on va creer des classes d oxygene
x11()
ind=which(x[,1]!=0)
plot(density(x[ind,1]))

#la densité nous informe d'une coupure bien sentie aux alentours de 50
## 2 classes
y=x
ind=which(x[,1]<50)
y[ind,1]="ano"
ind=which(x[,1]>=50)
y[ind,1]="oxy"

form1=as.formula(paste("as.factor(oxy)~ ","dens+moyve+maxve+deb+debv+deb7+deb15"))
ylda=lda(form1,data=y)
x11()
plot(ylda)

#prevision avec l echantillon de construction
ypred=predict(ylda)
ycla=y[,1]
misrate=1-sum(ypred$cla==ycla)/length(ycla)

###Avec un chantillon temoin : 1/3 des données
n=dim(y)[1]
nts=round(1/3*n,0)
ind=sample(1:n,nts,replace=F)
yts=y[ind,]
yls=y[-ind,]

###arbre classif
ytree=tree(form1,data=yls)
x11()
plot(ytree)
text(ytree)

##Prevision avec ech de construction
ypred=predict(ytree)
#regarder comment est fait ypred pr comprendre la suite

ind=apply(ypred,1,which.max)
ypred=colnames(ypred)[ind]

misrate=sum(yls[,1]!=ypred)/length(ypred)

##Avec echantillon test
ypredt=predict(ytree,yts)
ind=apply(ypredt,1,which.max)
ypredt=colnames(ypredt)[ind]

misrats=sum(yts[,1]!=ypredt)/length(ypredt)

##misrate<misrats evidemment

###Predicteurs agrégés (bagged predictors) et prevision
#Cette partie va mettre du temps a tourner car l operation de bagging est repetee
# 100 fois...
#Commencez par changer KK=10

KK=100
restot=matrix(-99,(KK-1),2)
n=nrow(y)
##cette fosi ci, 1/2 des données pr echantillon test
##ceci va permettre de montrer l amelioration de la convergence

nts=round(1/2*n,0)
ind=sample(1:n,nts,replace=F)
yts=y[ind,]
yls=y[-ind,]
nls=nrow(yls)

for(K in 2:KK){
    restree=matrix(-99,nts,K)
    reslda=matrix(-99,nts,K)
for (k in 1:K){
    indB=sample(1:nls,nls,replace=T)
    yB=yls[indB,]
#Avec un arbre
    ytree=tree(form1,data=yB)
    ytreepred=predict(ytree,yts)
    ind=apply(ytreepred,1,which.max)
    ytreepred=colnames(ytreepred)[ind]
    restree[,k]=ytreepred
#Avec discriminante
    ylda=lda(form1,data=yB)
    yldapred=predict(ylda,yts)
    ycla=as.character(yldapred$class)
    reslda[,k]=ycla
    }
##Avant de continuer on REGARDE le contenu de restree et rescla

#apply permet de renvoyer pour les lignes de restree une table de frequence
#sortie sous forme de liste
   pourtree=apply(restree,1,table)
#A partir de la liste precedente, je detecte pour chaque element de la liste
#le max
   pourtree=lapply(pourtree,which.max)
#je recupere le nom du max, hehe
   pourtree=lapply(pourtree,names)
#je mets la liste sous forme d un vecteur
   pourtree=unlist(pourtree)
#meme chose pour lda
   pourlda=apply(reslda,1,table)
   pourlda=lapply(pourlda,which.max)
   pourlda=lapply(pourlda,names)
   pourlda=unlist(pourlda)
#Je calcule maintenant les taux de mauvaise classif dans les deux cas
#avec ech temoin
   mcr.tree=sum(yts[,1]!=pourtree)/length(pourtree)
   mcr.lda=sum(yts[,1]!=pourlda)/length(pourlda)

####ceci est repete pour differentes valeurs de K (de 2 a KK) et stocké dans restot sous la forme
#d une matrice de dimension K-1 lignes et 2 colonnes
    restot[K-1,]=cbind(mcr.tree,mcr.lda)
}

##le graphique ci-dessous represente l evolution du taux de mauvaise classification pr arbre et lda
## qd le nombre K de predicteurs agreges augmente
#l arbre de regression donne une meilleure prevision dans ce cas et les taux de mauvaise classification
#se stabilisent au cours des iterations, ceci de maniere plus rapide pour lda

x11()
matplot(2:KK,restot, type="b",col=c(1,2), lty=1,pch=1,ylab="Misclassification rate",xlab="num. of aggregated predictors")