font size...
[iramuteq] / Rscripts / CHDKmeans.R
1 #Author: Pierre Ratinaud
2 #Copyright (c) 2008 Pierre Ratinaud
3 #Lisense: GNU/GPL
4
5
6
7 #library(ca)
8 #library(MASS)
9 #source('/home/pierre/workspace/iramuteq/Rscripts/afc.R')
10 #data<-read.table('output/corpus_bin.csv',header=TRUE,sep='\t')
11 #source('/home/pierre/workspace/iramuteq/Rscripts/anacor.R')
12 #library(ade4)
13 pp<-function(txt,val) {
14         d<-paste(txt,' : ')
15         print(paste(d,val))
16 }
17 MyChiSq<-function(x,sc,n){
18         sr<-rowSums(x)
19         E <- outer(sr, sc, "*")/n
20         STAT<-sum((abs(x - E))^2/E)
21         STAT
22 }
23
24 CHD<-function(data,x=9){
25 #       sink('/home/pierre/workspace/iramuteq/dev/findchi2.txt')
26         dataori=as.matrix(data)
27         row.names(dataori)=rownames(data)
28         dtable=as.matrix(data)
29         row.names(dtable)=rownames(data)
30         rowelim<-NULL
31         pp('ncol entree : ',ncol(dtable))
32         pp('nrow entree',nrow(dtable))
33         listcol <- list()
34         listmere <- list()
35         list_fille <- list()
36         print('vire colonnes vides en entree')#FIXME : il ne doit pas y avoir de colonnes vides en entree !!
37         sdt<-colSums(dtable)
38         if (min(sdt)==0)
39                 dtable<-dtable[,-which(sdt==0)]
40         mere<-1
41         for (i in 1:x) {
42                 clnb<-(i*2)
43                 listmere[[clnb]]<-mere
44                 listmere[[clnb+1]]<-mere
45                 list_fille[[mere]] <- c(clnb,clnb+1)
46                 listcol[[clnb]]<-vector()
47                 listcol[[clnb+1]]<-vector()
48                 #extraction du premier facteur de l'afc
49                 print('afc')
50                 #afc<-ca(dtable,nd=1)
51                 #afc<-corresp(dtable,nd=1)
52                 #afc<-fca(dtable)
53                 pp('taille dtable dans boucle (col/row)',c(ncol(dtable),nrow(dtable)))
54 #               afc<-boostana(dtable,nd=1)
55 #               #afc<-dudi.coa(dtable,scannf=FALSE,nf=1)
56 #               pp('SV',afc$singular.values)
57 #               pp('V.P.', afc$eigen.values)
58 #               #pp('V.P.',afc$eig[1])
59 #               #pp('V.P.',afc$factexpl[1])
60 #               #print(afc$chisq)
61 #               afcchi<-chisq.test(dtable)
62 #               Tinert<-afcchi$statistic/sum(dtable)
63 #               pp('inertie totale',Tinert)
64 #               #coordonnees des colonnes sur le premier facteur
65 #               #coordrow=afc$rowcoord
66 #               coordrow=as.matrix(afc$row.scores)
67 #               #coordrow<-as.matrix(afc$rproj[,1])
68 #               #coordrow<-as.matrix(afc$li)
69 #               coordrowori<-coordrow
70 #               #row.names(coordrow)<-afc$rownames
71 #               row.names(coordrow)<-rownames(dtable)
72 #               #classement en fonction de la position sur le premier facteur
73 #               print('deb recherche meilleur partition')
74 #               ordertable<-cbind(dtable,coordrow)
75 #               coordrow<-as.matrix(coordrow[order(coordrow[,1]),])
76 #               ordertable<-ordertable[order(ordertable[,ncol(ordertable)]),][,-ncol(ordertable)]
77 #               listinter<-vector()
78 #               listlim<-vector()
79                 TT=sum(dtable)
80 #               sc<-colSums(dtable)
81 #               sc1<-ordertable[1,]
82 #               sc2<-colSums(dtable)-ordertable[1,]
83 #               chitable<-rbind(sc1,sc2)
84 #               #listinter<-vector()
85 #               maxinter<-0
86 #               rmax<-NULL
87 ##################################################################
88 #               for (l in 2:(nrow(dtable)-2)){
89 #                       chitable[1,]<-chitable[1,]+ordertable[l,]
90 #                       chitable[2,]<-chitable[2,]-ordertable[l,]
91 ##                      sc1<-sc1+ordertable[l,]
92 ##                      sc2<-sc2-ordertable[l,]
93 ##                      chitable<-rbind(sc1,sc2)
94 #                       chi<-MyChiSq(chitable,sc,TT)
95 #                       if (chi>maxinter) {
96 #                               maxinter<-chi
97 #                               rmax<-l
98 #                       }
99 #               #       listinter<-append(listinter,MyChiSq(chitable,sc,TT))#,chisq.test(chitable)$statistic/TT)
100 #               }
101 #               #plot(listinter)
102 #               print('@@@@@@@@@@@@@@@@@@@@@@@@@@@@')
103 #               pp('max inter phase 1', maxinter/TT)#max(listinter))
104 #               print('@@@@@@@@@@@@@@@@@@@@@@@@@@@@')
105 #               #maxinter<-which.max(listinter)+1
106 #
107 #               listclasse<-ifelse(coordrowori<=coordrow[(rmax),1],clnb,clnb+1)
108 #               
109 #               cl<-listclasse
110 #
111 #               pp('TT',TT)
112 #               #dtable<-cbind(dtable,'cl'= as.vector(cl))
113 #
114 #               N1<-length(listclasse[listclasse==clnb])
115 #               N2<-length(listclasse[listclasse==clnb+1])
116 #               pp('N1',N1)
117 #               pp('N2',N2)
118 ####################################################################
119 ##                  reclassement des individus                     #
120 ####################################################################
121 #               malcl<-1000000000000
122 #               it<-0
123 #               listsub<-list()
124 #               #in boucle
125 #
126 #               while (malcl!=0 & N1>=5 & N2>=5) {
127 #                       it<-it+1
128 #                       listsub[[it]]<-vector()
129 #                       pp('nombre iteration',it)
130 #                       vdelta<-vector()
131 #                       #dtable[,'cl']<-cl
132 #                       t1<-dtable[cl==clnb,]#[,-ncol(dtable)]
133 #                       t2<-dtable[cl==clnb+1,]#[,-ncol(dtable)]
134 #                       ncolt<-ncol(t1)
135 #                       pp('ncolt',ncolt)
136 #                       sc1<-colSums(t1)
137 #                       sc2<-colSums(t2)
138 #                       sc<-sc1+sc2
139 #                       TT<-sum(dtable)
140 #                       chtableori<-rbind(sc1,sc2)
141 #                       interori<-MyChiSq(chtableori,sc,TT)/TT#chisq.test(chtableori)$statistic#/TT
142 #                       chtable<-chtableori
143 #                       pp('interori',interori)
144 #                       N1<-nrow(t1)
145 #                       N2<-nrow(t2)
146 #                       pp('N1',N1)
147 #                       pp('N2',N2)
148 #                       
149 #                       for (l in 1:nrow(dtable)){
150 #                                       if(cl[l]==clnb){
151 #                                               chtable[1,]<-sc1-dtable[l,]
152 #                                               chtable[2,]<-sc2+dtable[l,]
153 #                                       }else{
154 #                                               chtable[1,]<-sc1+dtable[l,]
155 #                                               chtable[2,]<-sc2-dtable[l,]
156 #                                       }
157 #                                       interswitch<-MyChiSq(chtable,sc,TT)/TT#chisq.test(chtable)$statistic/TT
158 #                                       ws<-interori-interswitch
159 #
160 #                                       if (ws<0){
161 #                                               interori<-interswitch
162 #                                               if(cl[l]==clnb){
163 #                                                       sc1<-chtable[1,]
164 #                                                       sc2<-chtable[2,]
165 #                                                       cl[l]<-clnb+1
166 #                                                       listsub[[it]]<-append(listsub[[it]],l)
167 #                                               }else{
168 #                                                       sc1<-chtable[1,]
169 #                                                       sc2<-chtable[2,]
170 #                                                       cl[l]<-clnb
171 #                                                       listsub[[it]]<-append(listsub[[it]],l)
172 #                                               }
173 #                                               vdelta<-append(vdelta,l)
174 #                                       }
175 #                               }
176 ##                      for (val in vdelta) {
177 ##                              if (cl[val]==clnb) {
178 ##                                      cl[val]<-clnb+1
179 ##                                      listsub[[it]]<-append(listsub[[it]],val)
180 ##                                      }else {
181 ##                                      cl[val]<-clnb
182 ##                                      listsub[[it]]<-append(listsub[[it]],val)
183 ##                              }
184 ##                      }
185 #                       print('###################################')
186 #                       print('longueur < 0')
187 #                       malcl<-length(vdelta)
188 #                       if ((it>1)&&(!is.logical(listsub[[it]]))){
189 #                               if (listsub[[it]]==listsub[[(it-1)]]){
190 #                                       malcl<-0
191 #                               }
192 #                       }
193 #                       print(malcl)
194 #                       print('###################################')
195 #               }
196 ##########################################################################
197 # kmeans
198 #        clori<-kmeans(dtable,2)$cluster
199 #        cl<-ifelse(clori==1,clnb,clnb+1)
200 #pam
201         library(cluster)
202         clori<-pam(dist(dtable,method='binary'),2,diss=TRUE)$cluster
203         cl<-ifelse(clori==1,clnb,clnb+1)
204         
205             dtable<-cbind(dtable,'cl'=as.vector(cl))
206
207 #######################################################################
208 #                 Fin reclassement des individus                      #
209 #######################################################################
210 #               if (!(length(cl[cl==clnb])==1 || length(cl[cl==clnb+1])==1)) {
211                         t1<-dtable[dtable[,'cl']==clnb,][,-ncol(dtable)]
212                         t2<-dtable[dtable[,'cl']==clnb+1,][,-ncol(dtable)]
213                 
214                         chtable<-rbind(colSums(t1),colSums(t2))
215                         inter<-chisq.test(chtable)$statistic/TT
216                         pp('last inter',inter)
217                         print('=====================')
218                         
219                         #calcul de la specificite des colonnes
220                         mint<-min(nrow(t1),nrow(t2))
221                         maxt<-max(nrow(t1),nrow(t2))
222                         seuil<-round((1.9*(maxt/mint))+1.9,digit=6)
223                         #sink('/home/pierre/workspace/iramuteq/dev/findchi2.txt')
224 #                       print('ATTENTION SEUIL 3,84')
225 #                       seuil<-3.84
226                         pp('seuil',seuil)
227                         sominf<-0
228                         nv<-0
229                         nz<-0
230                         ncclnb<-0
231                         ncclnbp<-0
232             NN1<-0
233             NN2<-0
234             maxchip<-0
235             nbzeroun<-0
236             res1<-0
237             res2<-0
238             nbseuil<-0
239             nbexe<-0
240             nbcontrib<-0
241                         cn<-colnames(dtable)[1:(ncol(dtable)-1)]
242                         for (k in 1:(ncol(dtable)-1)) {
243                                 #print(k)
244                                 tab<-table(dtable[,k],dtable[,ncol(dtable)])
245                                 #print(cn[k])
246                                 #print (tab)
247                                 #chidemoi<-(sum(t)*(((t[1,1]*t[2,2])-(t[1,2]*t[2,1]))^2))/((rowSums(t)[1])*(rowSums(t)[2])*(colSums(t)[1])*(colSums(t)[2]))
248                                 if (rownames(tab)[1]==1 & nrow(tab)==1) {
249                                         tab<-rbind(c(0,0),tab[1])
250                                         print('tableau vide dessus')
251                                 }
252                                 if (min(tab)<=10){
253                                         chi<-chisq.test(tab,correct=TRUE)
254                                 }else{
255                                         chi<-chisq.test(tab,correct=FALSE)
256                                 }
257                                 if ((min(tab[2,])==1) & (min(tab[1,])!=0)){
258                                         chi$statistic<-seuil+1
259                                 #    print('min tab 2 == 0')    
260                                 }
261 #                if(((tab[2,1]>tab[1,1]) & (tab[2,2]>tab[1,2]))){
262 #                                       nz<-nz+1
263 #                                       chi$statistic<-seuil-1
264 #                                       sominf<-sominf-1
265 #                               }
266                 if ((min(tab[1,])==0) & (min(tab[2,])!=0)){
267                                         chi$statistic<-seuil-1
268                 #    print('min tab 1 == 0')
269                                 }
270                                 if (is.na(chi$statistic)) {
271                                         chi$statistic<-0
272                                         print('NA dans chi$statistic')
273                                 }
274 #                               print('#####################')
275                                 if (chi$statistic>seuil){
276                                 i#      print('sup seuil')
277                                         if (tab[2,1]<chi$expected[2,1]){
278                                                 listcol[[clnb]]<-append(listcol[[clnb]],cn[k])
279                                                 ncclnb<-ncclnb+1
280 #                                               print('@@@@@@@@specifique de la classe@@@@@@')
281 #                                               print(clnb)
282                                         } else if (tab[2,2]<chi$expected[2,2]){
283                                                 listcol[[clnb+1]]<-append(listcol[[clnb+1]],cn[k])
284                                                 ncclnbp<-ncclnbp+1
285 #                                               print('@@@@@@@@specifique de la classe@@@@@@')
286 #                                               print(clnb+1)
287                                         }
288                                 }
289 #                               if (chi$statistic<seuil){
290 #                                       print('inf seuil')
291 #                                       print(tab)
292 #                                       print(chi$expected)
293 #                                       print(chi$statistic)
294 #                                       sominf<-sominf+1
295 #                               }
296                                 
297 #                               print('tableau')
298 #                               print(tab)
299 #                               chi<-chisq.test(tab,correct=FALSE)
300 #                               print('sans correction')
301 #                               print(chi$statistic)
302 #                               print(chi$expected)
303 #                               print(chi$residuals)
304 #                chit<-((tab-chi$expected)^2)/chi$expected
305 #                print(mean(chi$residuals))
306 #                print(sd(as.vector(chi$residuals)))
307                 #print(chi$residuals/sd(as.vector(chi$residuals)))
308                 #print(abs(chi$residuals/sd(as.vector(chi$residuals))))
309                 #print(qnorm(abs(chi$residuals/sd(as.vector(chi$residuals)))))
310 #                difftab<-tab-chi$expected
311  #               print('contribution')
312                 #print((difftab-mean(difftab))/sd(as.vector(difftab)))
313                 #print(abs((chit-mean(chit))/sd(as.vector(chit))))
314 #                tabin1<-matrix(1,2,2)
315 #                tabin1<-tabin1-(colSums(tab)/sum(tab))
316 #                tabin2<-matrix(1,2,2)
317 #                tabin2<-(tabin2-(rowSums(tab)/sum(tab)))
318 #                tabinv<-tabin1*tabin2
319 #                CS<-colSums(tab)
320 #                RS<-rowSums(tab)
321 #                GT<-sum(tab)
322 #                print(tabinv)
323 #                contrib<-(tab-chi$expected)/sqrt(chi$expected * ((1 - RS/GT) %*% t(1 - CS/GT)))
324 #                print(chit)
325 #                print('TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT')
326 #                print(contrib)
327 #                if (((which.max(chit)==2) || (which.max(chit)==4)) & chi$statistic>=seuil) {
328 #                    if (which.max(chit)==2) NN1<-NN1+1
329 #                    if (which.max(chit)==4) NN2<-NN2+1
330 #                }
331 #                if (max(abs(contrib[2,])>=1.96)) {
332 #                    nbcontrib<-nbcontrib+1
333 #                }
334 #                if (max(chit[2,]>=seuil)) maxchip<-maxchip+1
335 #                if (chi$statistic >=seuil) nbseuil<-nbseuil+1
336 #                if ((min(tab[2,])==0) & (chi$statistic<seuil)) nbzeroun<-nbzeroun+1
337 #                if ((((which.max(chi$residual)==2) || (which.max(chi$residual)==4) || (which.min(chi$residual)==2) || (which.min(chi$residual)==4))) & chi$statistic>=seuil) {
338 #                    if ((which.max(chi$residual)==2) || (which.min(chi$residual)==4)) res1<-res1+1
339 #                    if ((which.max(chi$residual)==4) || (which.min(chi$residual)==2)) res2<-res2+1
340 #                } else if (!((which.max(chi$residual)==2) || (which.max(chi$residual)==4) || (which.min(chi$residual)==2) || (which.min(chi$residual)==4)) & chi$statistic>=seuil) {
341 #                    print('AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA')
342 #                    print('AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA')
343 #                }
344 #                if (length(unique(as.vector(chi$residual)))==2) {
345 #                    nbexe<-nbexe+1
346 #                }
347
348 #                               print('avec correction')
349 #                if (min(tab)<5) chi<-chisq.test(tab,correct=TRUE)
350 #                               print(chi$statistic)
351 #                #if (chi$statistic >=seuil) nbseuil<-nbseuil+1
352 #                               print(chi$expected)
353 #                               print(chi$residuals)
354 #                chit<-((abs(tab-chi$expected)-0.5)^2)/chi$expected
355 #                print(chit)
356 #                               print('#####################')
357                         }
358 #            pp('NN1',NN1)
359 #            pp('NN2',NN2)
360 #            pp('maxchip',maxchip)
361 #            pp('nbzeroun',nbzeroun)
362 #            pp('res1',res1)
363 #            pp('res2',res2)
364 #            pp('nbseuil',nbseuil)
365 #            pp('nbexe',nbexe)
366 #            pp('nbcontrib', nbcontrib)
367                         print('resultats elim item')
368                         pp(clnb+1,length(listcol[[clnb+1]]))
369                         pp(clnb,length(listcol[[clnb]]))
370 #                       pp('inf seuil',sominf)
371 #                       pp('nv',nv)
372                         pp('ncclnb',ncclnb)
373                         pp('ncclnbp',ncclnbp)
374 #                       pp('nz',nz)
375                         #sink()
376                         #lignes concernees
377                         listrownamedtable<-rownames(dtable)
378                         listrownamedtable<-as.integer(listrownamedtable)
379                         newcol<-vector(length=nrow(dataori))
380                         #remplissage de la nouvelle colonne avec les nouvelles classes
381                         print('remplissage')
382                         num<-0
383                         for (ligne in listrownamedtable) {
384                                 num<-num+1
385                                 newcol[ligne]<-as.integer(as.vector(dtable[,'cl'][num])[1])
386                         }
387                         #recuperation de la classe precedante pour les cases vides
388                         print('recuperation classes precedentes')
389                         matori<-as.matrix(dataori)
390                         if (i!=1) {
391                                 #    options(warn=-1)
392                                 for (ligne in 1:length(newcol)) {
393                                         #       print(newcol[ligne])
394                                         if (newcol[ligne]==0) { # ce test renvoie un warning
395                                                 if (ligne %in% rowelim){
396                                                         newcol[ligne]<-0
397                                                 } else {
398                                                         newcol[ligne]<-matori[ligne,length(matori[1,])]
399                                                 }
400                                                 
401                                         }
402                                 }
403                                 #   options(warn=0)
404                         }
405                         dataori<-cbind(dataori,newcol)
406                         
407                         tailleclasse<-as.matrix(summary(as.factor(as.character(dataori[,ncol(dataori)]))))
408                         #tailleclasse<-as.integer(tailleclasse)
409                         print('tailleclasse')
410                         print(tailleclasse)
411                         tailleclasse<-as.matrix(tailleclasse[!(rownames(tailleclasse)==0),])
412                         plusgrand<-which.max(tailleclasse)
413                         #???????????????????????????????????
414                         #Si 2 classes ont des effectifs egaux, on prend la premiere de la liste...
415                         if (length(plusgrand)>1) {
416                                 plusgrand<-plusgrand[1]
417                         }
418                         #????????????????????????????????????
419                         
420                         #constuction du prochain tableau a analyser
421                         print('construction tableau suivant')
422                         classe<-as.integer(rownames(tailleclasse)[plusgrand])
423                         dtable<-dataori[dataori[,ncol(dataori)]==classe,]
424                         row.names(dtable)<-rownames(dataori[dataori[,ncol(dataori)]==classe,])
425                         dtable<-dtable[,1:(ncol(dtable)-i)]
426                         mere<-classe
427                         listcolelim<-listcol[[as.integer(classe)]]
428                         mother<-listmere[[as.integer(classe)]]
429                         while (mother!=1) {
430                                 listcolelim<-append(listcolelim,listcol[[mother]])
431                                 mother<-listmere[[mother]]
432                         }
433                         
434                         listcolelim<-sort(unique(listcolelim))
435                         pp('avant',ncol(dtable))
436                         if (!is.logical(listcolelim)){
437                                 print('elimination colonne')
438                                 #dtable<-dtable[,-listcolelim]
439                                 dtable<-dtable[,!(colnames(dtable) %in% listcolelim)]
440                         }
441                         pp('apres',ncol(dtable))
442                         #elimination des colonnes ne contenant que des 0
443                         print('vire colonne inf 3 dans boucle')
444                         sdt<-colSums(dtable)
445                         if (min(sdt)<=3)
446                                 dtable<-dtable[,-which(sdt<=3)]
447         
448                         #elimination des lignes ne contenant que des 0
449                         print('vire ligne vide dans boucle')
450                         if (ncol(dtable)==1) {
451                                 sdt<-dtable[,1]
452                         } else {
453                                 sdt<-rowSums(dtable)
454                         }
455                         if (min(sdt)==0) {
456                                 rowelim<-as.integer(rownames(dtable)[which(sdt==0)])
457                                 print('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&')
458                                 print(rowelim)
459                                 print('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&')
460                                 dtable<-dtable[-which(sdt==0),]
461                         }
462 #               }
463         }
464 #       sink()
465         res <- list(n1 = dataori[,(ncol(dataori)-x+1):ncol(dataori)], list_mere = listmere, list_fille = list_fille)
466         res
467 }
468
469 #dataout<-CHD(data,9)
470
471 #igraph pour graph de relation