iramuteq.org Git - iramuteq/blob - Rscripts/CHD.R

   1 #Author: Pierre Ratinaud
   2 #Copyright (c) 2008-2011 Pierre Ratinaud
   3 #License: GNU/GPL
   4
   5 pp<-function(txt,val) {
   6         d<-paste(txt,' : ')
   7         print(paste(d,val))
   8 }
   9 MyChiSq<-function(x,sc,n){
  10         sr<-rowSums(x)
  11         E <- outer(sr, sc, "*")/n
  12         STAT<-sum(((x - E)^2)/E)
  13         STAT
  14 }
  15
  16 MySpeedChi <- function(x,sc) {
  17     sr <-rowSums(x)
  18     E <- outer(sr, sc, "*")
  19         STAT<-sum((x - E)^2/E)
  20         STAT
  21 }
  22
  23 find.max <- function(dtable, chitable, compte, rmax, maxinter, sc, TT) {
  24     ln <- which(dtable==1, arr.ind=TRUE)
  25     lo <- list()
  26     lo[1:nrow(dtable)] <- 0
  27     for (k in 1:nrow(ln)) {lo[[ln[k,1]]]<-append(lo[[ln[k,1]]],ln[k,2])}
  28     for (k in 1:nrow(dtable)) {lo[[k]] <- lo[[k]][-1]}
  29         ## lo<-lo[-c(1,length(lo))]
  30         ## for (l in lo) {
  31         ##     compte <- compte + 1
  32         ##     chitable[1,l]<-chitable[1,l]+1
  33         ##     chitable[2,l]<-chitable[2,l]-1
  34         ##     chi<-MyChiSq(chitable,sc,TT)
  35                 ## if (chi>maxinter) {
  36                 ##     maxinter<-chi
  37                 ##     rmax<-compte
  38                 ## }
  39     #}
  40         lo<-lo[-c(1)]
  41         for (l in lo) {
  42                 chi<-MyChiSq(chitable,sc,TT)
  43                 if (chi>maxinter) {
  44                         maxinter<-chi
  45                         rmax<-compte
  46                 }
  47                 compte <- compte + 1
  48                 chitable[1,l]<-chitable[1,l]+1
  49                 chitable[2,l]<-chitable[2,l]-1
  50         }
  51     res <- list(maxinter=maxinter, rmax=rmax)
  52     res
  53 }
  54
  55
  56
  57
  58
  59 CHD<-function(data.in, x=9, mode.patate = FALSE, svd.method, libsvdc.path=NULL){
  60 #       sink('/home/pierre/workspace/iramuteq/dev/findchi2.txt')
  61         dataori <- data.in
  62     row.names(dataori) <- rownames(data.in)
  63         dtable <- data.in
  64         colnames(dtable) <- 1:ncol(dtable)
  65     dout <- NULL
  66         rowelim<-NULL
  67         pp('ncol entree : ',ncol(dtable))
  68         pp('nrow entree',nrow(dtable))
  69         listcol <- list()
  70         listmere <- list()
  71         list_fille <- list()
  72         print('vire colonnes vides en entree')#FIXME : il ne doit pas y avoir de colonnes vides en entree !!
  73         sdt<-colSums(dtable)
  74         if (min(sdt)==0)
  75                 dtable<-dtable[,-which(sdt==0)]
  76     print('vire lignes vides en entree')
  77     sdt<-rowSums(dtable)
  78         if (min(sdt)==0) {
  79         rowelim<-as.integer(rownames(dtable)[which(sdt==0)])
  80         print('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&')
  81         print(rowelim)
  82         print('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&')
  83                 dtable<-dtable[-which(sdt==0),]
  84         }
  85         mere<-1
  86         for (i in 1:x) {
  87                 clnb<-(i*2)
  88                 listmere[[clnb]]<-mere
  89                 listmere[[clnb+1]]<-mere
  90                 list_fille[[mere]] <- c(clnb,clnb+1)
  91                 listcol[[clnb]]<-vector()
  92                 listcol[[clnb+1]]<-vector()
  93                 #extraction du premier facteur de l'afc
  94                 print('afc')
  95                 pp('taille dtable dans boucle (col/row)',c(ncol(dtable),nrow(dtable)))
  96                 afc<-boostana(dtable, nd=1, svd.method = svd.method, libsvdc.path=libsvdc.path)
  97                 pp('SV',afc$singular.values)
  98                 pp('V.P.', afc$eigen.values)
  99                 coordrow <- as.matrix(afc$row.scores[,1])
 100                 coordrowori<-coordrow
 101                 row.names(coordrow)<-rownames(dtable)
 102         coordrow <- cbind(coordrow,1:nrow(dtable))
 103                 print('deb recherche meilleur partition')
 104         ordert <- as.matrix(coordrow[order(coordrow[,1]),])
 105         ordert <- cbind(ordert, 1:nrow(ordert))
 106         ordert <- ordert[order(ordert[,2]),]
 107
 108                 listinter<-vector()
 109                 listlim<-vector()
 110         dtable <- dtable[order(ordert[,3]),]
 111         sc <- colSums(dtable)
 112         TT <- sum(sc)
 113         sc1 <- dtable[1,]
 114         sc2 <- colSums(dtable) - sc1
 115         chitable <- rbind(sc1, sc2)
 116         compte <- 1
 117         maxinter <- 0
 118         rmax <- NULL
 119
 120         inert <- find.max(dtable, chitable, compte, rmax, maxinter, sc, TT)
 121         print('@@@@@@@@@@@@@@@@@@@@@@@@@@@@')
 122                 pp('max inter phase 1', inert$maxinter/TT)#max(listinter))
 123                 print('@@@@@@@@@@@@@@@@@@@@@@@@@@@@')
 124         ordert <- ordert[order(ordert[,3]),]
 125                 listclasse<-ifelse(coordrowori<=ordert[(inert$rmax),1],clnb,clnb+1)
 126             dtable <- dtable[order(ordert[,2]),]
 127                 cl<-listclasse
 128                 pp('TT',TT)
 129                 #dtable<-cbind(dtable,'cl'= as.vector(cl))
 130
 131                 N1<-length(listclasse[listclasse==clnb])
 132                 N2<-length(listclasse[listclasse==clnb+1])
 133                 pp('N1',N1)
 134                 pp('N2',N2)
 135 ###################################################################
 136 #                  reclassement des individus                     #
 137 ###################################################################
 138         if (!mode.patate) {
 139                 malcl<-1000000000000
 140                 it<-0
 141                 listsub<-list()
 142                 #in boucle
 143             ln <- which(dtable==1, arr.ind=TRUE)
 144             lnz <- list()
 145             lnz[1:nrow(dtable)] <- 0
 146
 147             for (k in 1:nrow(ln)) {lnz[[ln[k,1]]]<-append(lnz[[ln[k,1]]],ln[k,2])}
 148             for (k in 1:nrow(dtable)) {lnz[[k]] <- lnz[[k]][-1]}
 149                 TT<-sum(dtable)
 150
 151                 while (malcl!=0 & N1>=5 & N2>=5) {
 152                         it<-it+1
 153                         listsub[[it]]<-vector()
 154                 txt <- paste('nombre iteration', it)
 155                         #pp('nombre iteration',it)
 156                         vdelta<-vector()
 157                         #dtable[,'cl']<-cl
 158                         t1<-dtable[which(cl[,1]==clnb),]#[,-ncol(dtable)]
 159                         t2<-dtable[which(cl[,1]==clnb+1),]#[,-ncol(dtable)]
 160                         ncolt<-ncol(t1)
 161                         #pp('ncolt',ncolt)
 162
 163                 if (N1 != 1) {
 164                     sc1<-colSums(t1)
 165                 } else {
 166                     sc1 <- t1
 167                 }
 168                 if (N2 != 1) {
 169                             sc2<-colSums(t2)
 170                 } else {
 171                     sc2 <- t2
 172                 }
 173
 174                 sc<-sc1+sc2
 175                         chtableori<-rbind(sc1,sc2)
 176                         chtable<-chtableori
 177                         interori<-MyChiSq(chtableori,sc,TT)/TT#chisq.test(chtableori)$statistic#/TT
 178                         txt <- paste(txt, ' - interori : ',interori)
 179                 #pp('interori',interori)
 180
 181                         N1<-nrow(t1)
 182                         N2<-nrow(t2)
 183
 184                         #pp('N1',N1)
 185                         #pp('N2',N2)
 186                         txt <- paste(txt, 'N1:', N1,'-N2:',N2)
 187                 print(txt)
 188                 compte <- 0
 189                         for (l in lnz){
 190                         chi.in<-chtable
 191                         compte <- compte + 1
 192                                         if(cl[compte]==clnb){
 193                                                 chtable[1,l]<-chtable[1,l]-1
 194                                                 chtable[2,l]<-chtable[2,l]+1
 195                                         }else{
 196                                                 chtable[1,l]<-chtable[1,l]+1
 197                                                 chtable[2,l]<-chtable[2,l]-1
 198                                         }
 199                                         interswitch<-MyChiSq(chtable,sc,TT)/TT#chisq.test(chtable)$statistic/TT
 200                                         ws<-interori-interswitch
 201
 202                                         if (ws<0){
 203                                                 interori<-interswitch
 204                                                 if(cl[compte]==clnb){
 205                                                         #sc1<-chtable[1,]
 206                                                         #sc2<-chtable[2,]
 207                                                         cl[compte]<-clnb+1
 208                                                         listsub[[it]]<-append(listsub[[it]],compte)
 209                                                 } else {
 210                                                         #sc1<-chtable[1,]
 211                                                         #sc2<-chtable[2,]
 212                                                         cl[compte]<-clnb
 213                                                         listsub[[it]]<-append(listsub[[it]],compte)
 214                                                 }
 215                                                 vdelta<-append(vdelta,compte)
 216                         } else {
 217                             chtable<-chi.in
 218                                         }
 219                    }
 220     #                   for (val in vdelta) {
 221     #                           if (cl[val]==clnb) {
 222     #                                   cl[val]<-clnb+1
 223     #                                   listsub[[it]]<-append(listsub[[it]],val)
 224     #                                   }else {
 225     #                                   cl[val]<-clnb
 226     #                                   listsub[[it]]<-append(listsub[[it]],val)
 227     #                           }
 228     #                   }
 229                         print('###################################')
 230                         print('longueur < 0')
 231                         malcl<-length(vdelta)
 232                         if ((it>1)&&(!is.logical(listsub[[it]]))&&(!is.logical(listsub[[it-1]]))){
 233                                 if (listsub[[it]]==listsub[[(it-1)]]){
 234                                         malcl<-0
 235                                 }
 236                         }
 237                         print(malcl)
 238                         print('###################################')
 239                 }
 240         }
 241                 #dtable<-cbind(dtable,'cl'=as.vector(cl))
 242         #dtable[,'cl'] <-as.vector(cl)
 243 #######################################################################
 244 #                 Fin reclassement des individus                      #
 245 #######################################################################
 246 #               if (!(length(cl[cl==clnb])==1 || length(cl[cl==clnb+1])==1)) {
 247                         #t1<-dtable[dtable[,'cl']==clnb,][,-ncol(dtable)]
 248                         #t2<-dtable[dtable[,'cl']==clnb+1,][,-ncol(dtable)]
 249                     t1<-dtable[which(cl[,1]==clnb),]#[,-ncol(dtable)]
 250                         t2<-dtable[which(cl[,1]==clnb+1),]#[,-ncol(dtable)]
 251             if (class(t1)=='numeric') {
 252                 sc1 <- as.vector(t1)
 253                 nrowt1 <- 1
 254             } else {
 255                 sc1 <- colSums(t1)
 256                 nrowt1 <- nrow(t1)
 257             }
 258             if (class(t2)=='numeric') {
 259                 sc2 <- as.vector(t2)
 260                 nrowt2 <- 1
 261             } else {
 262                 sc2 <- colSums(t2)
 263                 nrowt2 <- nrow(t2)
 264             }
 265             chtable<-rbind(sc1,sc2)
 266                         inter<-chisq.test(chtable)$statistic/TT
 267                         pp('last inter',inter)
 268                         print('=====================')
 269                         #calcul de la specificite des colonnes
 270                         mint<-min(nrowt1,nrowt2)
 271                         maxt<-max(nrowt1,nrowt2)
 272                         seuil<-round((1.9*(maxt/mint))+1.9,digit=6)
 273                         #sink('/home/pierre/workspace/iramuteq/dev/findchi2.txt')
 274 #                       print('ATTENTION SEUIL 3,84')
 275 #                       seuil<-3.84
 276                         pp('seuil',seuil)
 277                         sominf<-0
 278                         nv<-0
 279                         nz<-0
 280                         ncclnb<-0
 281                         ncclnbp<-0
 282             NN1<-0
 283             NN2<-0
 284             maxchip<-0
 285             nbzeroun<-0
 286             res1<-0
 287             res2<-0
 288             nbseuil<-0
 289             nbexe<-0
 290             nbcontrib<-0
 291                         cn<-colnames(dtable)
 292             #another try#########################################
 293             one <- cbind(sc1,sc2)
 294             cols <- c(length(which(cl==clnb)), length(which(cl==clnb+1)))
 295             print(cols)
 296             colss <- matrix(rep(cols,ncol(dtable)), ncol=2, byrow=TRUE)
 297             zero <- colss - one
 298             rows <- cbind(rowSums(zero), rowSums(one))
 299             n <- sum(cols)
 300             for (m in 1:nrow(rows)) {
 301                 obs <- t(matrix(c(zero[m,],one[m,]),2,2))
 302                 E <- outer(rows[m,],cols,'*')/n
 303                 if ((min(obs[2,])==0) & (min(obs[1,])!=0)) {
 304                     chi <- seuil + 1
 305                 } else if ((min(obs[1,])==0) & (min(obs[2,])!=0)) {
 306                     chi <- seuil - 1
 307                 } else if (any(obs < 10)) {
 308                     chi <- sum((abs(obs - E) - 0.5)^2 / E)
 309                 } else {
 310                     chi <- sum(((obs - E)^2)/E)
 311                 }
 312                 if (is.na(chi)) {
 313                     chi <- 0
 314                 }
 315                 if (chi > seuil) {
 316                     if (obs[2,1] < E[2,1]) {
 317                         listcol[[clnb]]<-append(listcol[[clnb]],cn[m])
 318                         ncclnb<-ncclnb+1
 319                     } else if (obs[2,2] < E[2,2]) {
 320                                                 listcol[[clnb+1]]<-append(listcol[[clnb+1]],cn[m])
 321                                                 ncclnbp<-ncclnbp+1
 322                     }
 323                 }
 324             }
 325             ######################################################
 326                         print('resultats elim item')
 327                         pp(clnb+1,length(listcol[[clnb+1]]))
 328                         pp(clnb,length(listcol[[clnb]]))
 329                         pp('ncclnb',ncclnb)
 330                         pp('ncclnbp',ncclnbp)
 331                         listrownamedtable<-rownames(dtable)
 332                         listrownamedtable<-as.integer(listrownamedtable)
 333                         newcol<-vector(length=nrow(dataori))
 334                         #remplissage de la nouvelle colonne avec les nouvelles classes
 335                         print('remplissage')
 336 #                       num<-0
 337             newcol[listrownamedtable] <- cl[,1]
 338                         #recuperation de la classe precedante pour les cases vides
 339                         print('recuperation classes precedentes')
 340             if (i!=1) {
 341                 newcol[which(newcol==0)] <- dout[,ncol(dout)][which(newcol==0)]
 342             }
 343             if(!is.null(rowelim)) {
 344                 newcol[rowelim] <- 0
 345             }
 346                         tailleclasse<-as.matrix(summary(as.factor(as.character(newcol))))
 347                         print('tailleclasse')
 348                         print(tailleclasse)
 349                         tailleclasse<-as.matrix(tailleclasse[!(rownames(tailleclasse)==0),])
 350                         plusgrand<-which.max(tailleclasse)
 351                         #???????????????????????????????????
 352                         #Si 2 classes ont des effectifs egaux, on prend la premiere de la liste...
 353                         if (length(plusgrand)>1) {
 354                                 plusgrand<-plusgrand[1]
 355                         }
 356                         #????????????????????????????????????
 357
 358                         #constuction du prochain tableau a analyser
 359                         print('construction tableau suivant')
 360             dout<-cbind(dout,newcol)
 361                         classe<-as.integer(rownames(tailleclasse)[plusgrand])
 362                         dtable<-dataori[which(newcol==classe),]
 363                         row.names(dtable)<-rownames(dataori)[which(newcol==classe)]
 364             colnames(dtable) <- 1:ncol(dtable)
 365                         mere<-classe
 366                         listcolelim<-listcol[[as.integer(classe)]]
 367                         mother<-listmere[[as.integer(classe)]]
 368                         while (mother!=1) {
 369                                 listcolelim<-append(listcolelim,listcol[[mother]])
 370                                 mother<-listmere[[mother]]
 371                         }
 372
 373                         listcolelim<-sort(unique(listcolelim))
 374                         pp('avant',ncol(dtable))
 375                         if (!is.logical(listcolelim)){
 376                                 print('elimination colonne')
 377                                 #dtable<-dtable[,-listcolelim]
 378                 dtable<-dtable[,!(colnames(dtable) %in% listcolelim)]
 379                         }
 380                         pp('apres',ncol(dtable))
 381                         #elimination des colonnes ne contenant que des 0
 382                         print('vire colonne inf 3 dans boucle')
 383                         sdt<-colSums(dtable)
 384                         if (min(sdt)<=3)
 385                                 dtable<-dtable[,-which(sdt<=3)]
 386
 387                         #elimination des lignes ne contenant que des 0
 388                         print('vire ligne vide dans boucle')
 389                         if (ncol(dtable)==1) {
 390                                 sdt<-dtable[,1]
 391                         } else {
 392                                 sdt<-rowSums(dtable)
 393                         }
 394                         if (min(sdt)==0) {
 395                                 rowelim<-as.integer(rownames(dtable)[which(sdt==0)])
 396                                 print('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&')
 397                                 print(rowelim)
 398                                 print('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&')
 399                                 dtable<-dtable[-which(sdt==0),]
 400                         }
 401 #               }
 402         }
 403 #       sink()
 404         res <- list(n1 = dout, list_mere = listmere, list_fille = list_fille)
 405         res
 406 }