code stacie

# predire : log de la constante d'affinite des molecules pour la HSA
# 1666 descripteurs
# 94 individus

# elminer variables inutiles
# justifier methodes choisies
# validation croisee pour determiner les parametres a choisir
# avantages/inconvenients des methodes utilisees
# proposer perspectives

#----------------------------- Importer les donnees --------------------------#

d = read.table('X.txt', header=T)
d = na.omit(d)

cte = read.table('Y.txt', header=T)

#---------------------------- Nettoyer les donnees ---------------------------#

# on enleve les descripteurs qui ont une variance nulle
# car pas informatifs
d = d[,sapply(d,sd)!=0]
# il reste 1492 descripteurs

d.cte = cbind(d,cte) # ajout de la colonne 'x'
# 1493 colonnes, 1493eme colonne contient le log des constantes d'affinite 'x'

#----- elimination des variables correlees

elimcor_sansY<-function(X,s=0.95)
{
  #X matrice contenant les variables a grouper
  #Y vecteur contenant les groupes a predire
  #s valeur seuil de correlation
  correl=cor(X)
  stop=F
  possetap=1:ncol(X)
  groupes=as.list(1:ncol(X))

  while (stop==F)
  {
    ##regroupement des var pour lesquelles |corr|>0.95
    gplist<-list(NULL)
    possglob=1:ncol(correl)
    for (i in 1:(ncol(correl)))
    {
      poss=possglob[-i]
      gplist[[i]]=c(i,poss[abs(correl[i,poss])>s])
    }
    ##on trie les groupes du plus gros au plus petit
    gplisteff=unlist(lapply(gplist,length))
    if (any(gplisteff>1))
    {
      gplistfin=gplist[gplisteff>1]
      gplistuniq=unlist(gplist[gplisteff==1])
      gpsel=NULL
      ##on selectionne dans chaque groupe une variable au hasard
      for (i in 1:length(gplistfin))
      {
        selloc=min(gplistfin[[i]])
        gploc=groupes[[possetap[selloc]]]
        for (j in 1:length(gplistfin[[i]]))
        {
          gploc=c(gploc,groupes[[possetap[gplistfin[[i]][j]]]])                 }
        groupes[[possetap[selloc]]]=unique(gploc)
        gpsel=c(gpsel,selloc)
      }
      possetap=possetap[c(gplistuniq,unique(gpsel))]
      correl=cor(X[,possetap])
    }
    else stop=T
  }
  #groupeseff=unlist(lapply(groupes,length))
  #groupes=groupes[groupeseff>1]
  return(list(possetap=possetap,groupes=groupes))
}

# elimcor sans le Y sur la matrice normalisee
d.norm = scale(d.cte, center=T, scale=T)
d.elimcor = d.norm[,elimcor_sansY(d.norm[,-1493],0.95)$possetap]
# reste 661 descripteurs