1 #Author: Pierre Ratinaud
2 #Copyright (c) 2008-2011 Pierre Ratinaud
5 pp<-function(txt,val) {
9 MyChiSq<-function(x,sc,n){
11 E <- outer(sr, sc, "*")/n
12 STAT<-sum(((x - E)^2)/E)
16 MySpeedChi <- function(x,sc) {
18 E <- outer(sr, sc, "*")
19 STAT<-sum((x - E)^2/E)
23 find.max <- function(dtable, chitable, compte, rmax, maxinter, sc, TT) {
24 ln <- which(dtable==1, arr.ind=TRUE)
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))]
32 chitable[1,l]<-chitable[1,l]+1
33 chitable[2,l]<-chitable[2,l]-1
34 chi<-MyChiSq(chitable,sc,TT)
40 res <- list(maxinter=maxinter, rmax=rmax)
48 CHD<-function(data.in, x=9, mode.patate = FALSE, svd.method, libsvdc.path=NULL){
49 # sink('/home/pierre/workspace/iramuteq/dev/findchi2.txt')
51 row.names(dataori) <- rownames(data.in)
53 colnames(dtable) <- 1:ncol(dtable)
56 pp('ncol entree : ',ncol(dtable))
57 pp('nrow entree',nrow(dtable))
61 print('vire colonnes vides en entree')#FIXME : il ne doit pas y avoir de colonnes vides en entree !!
64 dtable<-dtable[,-which(sdt==0)]
65 print('vire lignes vides en entree')
68 rowelim<-as.integer(rownames(dtable)[which(sdt==0)])
69 print('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&')
71 print('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&')
72 dtable<-dtable[-which(sdt==0),]
77 listmere[[clnb]]<-mere
78 listmere[[clnb+1]]<-mere
79 list_fille[[mere]] <- c(clnb,clnb+1)
80 listcol[[clnb]]<-vector()
81 listcol[[clnb+1]]<-vector()
82 #extraction du premier facteur de l'afc
84 pp('taille dtable dans boucle (col/row)',c(ncol(dtable),nrow(dtable)))
85 afc<-boostana(dtable, nd=1, svd.method = svd.method, libsvdc.path=libsvdc.path)
86 pp('SV',afc$singular.values)
87 pp('V.P.', afc$eigen.values)
88 coordrow <- as.matrix(afc$row.scores[,1])
90 row.names(coordrow)<-rownames(dtable)
91 coordrow <- cbind(coordrow,1:nrow(dtable))
92 print('deb recherche meilleur partition')
93 ordert <- as.matrix(coordrow[order(coordrow[,1]),])
94 ordert <- cbind(ordert, 1:nrow(ordert))
95 ordert <- ordert[order(ordert[,2]),]
99 dtable <- dtable[order(ordert[,3]),]
100 sc <- colSums(dtable)
103 sc2 <- colSums(dtable) - sc1
104 chitable <- rbind(sc1, sc2)
109 inert <- find.max(dtable, chitable, compte, rmax, maxinter, sc, TT)
110 print('@@@@@@@@@@@@@@@@@@@@@@@@@@@@')
111 pp('max inter phase 1', inert$maxinter/TT)#max(listinter))
112 print('@@@@@@@@@@@@@@@@@@@@@@@@@@@@')
113 ordert <- ordert[order(ordert[,3]),]
114 listclasse<-ifelse(coordrowori<=ordert[(inert$rmax),1],clnb,clnb+1)
115 dtable <- dtable[order(ordert[,2]),]
118 #dtable<-cbind(dtable,'cl'= as.vector(cl))
120 N1<-length(listclasse[listclasse==clnb])
121 N2<-length(listclasse[listclasse==clnb+1])
124 ###################################################################
125 # reclassement des individus #
126 ###################################################################
132 ln <- which(dtable==1, arr.ind=TRUE)
134 lnz[1:nrow(dtable)] <- 0
136 for (k in 1:nrow(ln)) {lnz[[ln[k,1]]]<-append(lnz[[ln[k,1]]],ln[k,2])}
137 for (k in 1:nrow(dtable)) {lnz[[k]] <- lnz[[k]][-1]}
140 while (malcl!=0 & N1>=5 & N2>=5) {
142 listsub[[it]]<-vector()
143 txt <- paste('nombre iteration', it)
144 #pp('nombre iteration',it)
147 t1<-dtable[which(cl[,1]==clnb),]#[,-ncol(dtable)]
148 t2<-dtable[which(cl[,1]==clnb+1),]#[,-ncol(dtable)]
164 chtableori<-rbind(sc1,sc2)
166 interori<-MyChiSq(chtableori,sc,TT)/TT#chisq.test(chtableori)$statistic#/TT
167 txt <- paste(txt, ' - interori : ',interori)
168 #pp('interori',interori)
175 txt <- paste(txt, 'N1:', N1,'-N2:',N2)
181 if(cl[compte]==clnb){
182 chtable[1,l]<-chtable[1,l]-1
183 chtable[2,l]<-chtable[2,l]+1
185 chtable[1,l]<-chtable[1,l]+1
186 chtable[2,l]<-chtable[2,l]-1
188 interswitch<-MyChiSq(chtable,sc,TT)/TT#chisq.test(chtable)$statistic/TT
189 ws<-interori-interswitch
192 interori<-interswitch
193 if(cl[compte]==clnb){
197 listsub[[it]]<-append(listsub[[it]],compte)
202 listsub[[it]]<-append(listsub[[it]],compte)
204 vdelta<-append(vdelta,compte)
209 # for (val in vdelta) {
210 # if (cl[val]==clnb) {
212 # listsub[[it]]<-append(listsub[[it]],val)
215 # listsub[[it]]<-append(listsub[[it]],val)
218 print('###################################')
219 print('longueur < 0')
220 malcl<-length(vdelta)
221 if ((it>1)&&(!is.logical(listsub[[it]]))&&(!is.logical(listsub[[it-1]]))){
222 if (listsub[[it]]==listsub[[(it-1)]]){
227 print('###################################')
230 #dtable<-cbind(dtable,'cl'=as.vector(cl))
231 #dtable[,'cl'] <-as.vector(cl)
232 #######################################################################
233 # Fin reclassement des individus #
234 #######################################################################
235 # if (!(length(cl[cl==clnb])==1 || length(cl[cl==clnb+1])==1)) {
236 #t1<-dtable[dtable[,'cl']==clnb,][,-ncol(dtable)]
237 #t2<-dtable[dtable[,'cl']==clnb+1,][,-ncol(dtable)]
238 t1<-dtable[which(cl[,1]==clnb),]#[,-ncol(dtable)]
239 t2<-dtable[which(cl[,1]==clnb+1),]#[,-ncol(dtable)]
240 if (class(t1)=='numeric') {
247 if (class(t2)=='numeric') {
254 chtable<-rbind(sc1,sc2)
255 inter<-chisq.test(chtable)$statistic/TT
256 pp('last inter',inter)
257 print('=====================')
258 #calcul de la specificite des colonnes
259 mint<-min(nrowt1,nrowt2)
260 maxt<-max(nrowt1,nrowt2)
261 seuil<-round((1.9*(maxt/mint))+1.9,digit=6)
262 #sink('/home/pierre/workspace/iramuteq/dev/findchi2.txt')
263 # print('ATTENTION SEUIL 3,84')
281 #another try#########################################
282 one <- cbind(sc1,sc2)
283 cols <- c(length(which(cl==clnb)), length(which(cl==clnb+1)))
285 colss <- matrix(rep(cols,ncol(dtable)), ncol=2, byrow=TRUE)
287 rows <- cbind(rowSums(zero), rowSums(one))
289 for (m in 1:nrow(rows)) {
290 obs <- t(matrix(c(zero[m,],one[m,]),2,2))
291 E <- outer(rows[m,],cols,'*')/n
292 if ((min(obs[2,])==0) & (min(obs[1,])!=0)) {
294 } else if ((min(obs[1,])==0) & (min(obs[2,])!=0)) {
296 } else if (any(obs < 10)) {
297 chi <- sum((abs(obs - E) - 0.5)^2 / E)
299 chi <- sum(((obs - E)^2)/E)
305 if (obs[2,1] < E[2,1]) {
306 listcol[[clnb]]<-append(listcol[[clnb]],cn[m])
308 } else if (obs[2,2] < E[2,2]) {
309 listcol[[clnb+1]]<-append(listcol[[clnb+1]],cn[m])
314 ######################################################
315 print('resultats elim item')
316 pp(clnb+1,length(listcol[[clnb+1]]))
317 pp(clnb,length(listcol[[clnb]]))
319 pp('ncclnbp',ncclnbp)
320 listrownamedtable<-rownames(dtable)
321 listrownamedtable<-as.integer(listrownamedtable)
322 newcol<-vector(length=nrow(dataori))
323 #remplissage de la nouvelle colonne avec les nouvelles classes
326 newcol[listrownamedtable] <- cl[,1]
327 #recuperation de la classe precedante pour les cases vides
328 print('recuperation classes precedentes')
330 newcol[which(newcol==0)] <- dout[,ncol(dout)][which(newcol==0)]
332 if(!is.null(rowelim)) {
335 tailleclasse<-as.matrix(summary(as.factor(as.character(newcol))))
336 print('tailleclasse')
338 tailleclasse<-as.matrix(tailleclasse[!(rownames(tailleclasse)==0),])
339 plusgrand<-which.max(tailleclasse)
340 #???????????????????????????????????
341 #Si 2 classes ont des effectifs egaux, on prend la premiere de la liste...
342 if (length(plusgrand)>1) {
343 plusgrand<-plusgrand[1]
345 #????????????????????????????????????
347 #constuction du prochain tableau a analyser
348 print('construction tableau suivant')
349 dout<-cbind(dout,newcol)
350 classe<-as.integer(rownames(tailleclasse)[plusgrand])
351 dtable<-dataori[which(newcol==classe),]
352 row.names(dtable)<-rownames(dataori)[which(newcol==classe)]
353 colnames(dtable) <- 1:ncol(dtable)
355 listcolelim<-listcol[[as.integer(classe)]]
356 mother<-listmere[[as.integer(classe)]]
358 listcolelim<-append(listcolelim,listcol[[mother]])
359 mother<-listmere[[mother]]
362 listcolelim<-sort(unique(listcolelim))
363 pp('avant',ncol(dtable))
364 if (!is.logical(listcolelim)){
365 print('elimination colonne')
366 #dtable<-dtable[,-listcolelim]
367 dtable<-dtable[,!(colnames(dtable) %in% listcolelim)]
369 pp('apres',ncol(dtable))
370 #elimination des colonnes ne contenant que des 0
371 print('vire colonne inf 3 dans boucle')
374 dtable<-dtable[,-which(sdt<=3)]
376 #elimination des lignes ne contenant que des 0
377 print('vire ligne vide dans boucle')
378 if (ncol(dtable)==1) {
384 rowelim<-as.integer(rownames(dtable)[which(sdt==0)])
385 print('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&')
387 print('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&')
388 dtable<-dtable[-which(sdt==0),]
393 res <- list(n1 = dout, list_mere = listmere, list_fille = list_fille)