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)
44 CHD<-function(data.in, x=9, libsvdc=FALSE, libsvdc.path=NULL){
45 # sink('/home/pierre/workspace/iramuteq/dev/findchi2.txt')
47 row.names(dataori) <- rownames(data.in)
49 colnames(dtable) <- 1:ncol(dtable)
52 pp('ncol entree : ',ncol(dtable))
53 pp('nrow entree',nrow(dtable))
57 print('vire colonnes vides en entree')#FIXME : il ne doit pas y avoir de colonnes vides en entree !!
60 dtable<-dtable[,-which(sdt==0)]
61 print('vire lignes vides en entree')
64 rowelim<-as.integer(rownames(dtable)[which(sdt==0)])
65 print('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&')
67 print('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&')
68 dtable<-dtable[-which(sdt==0),]
73 listmere[[clnb]]<-mere
74 listmere[[clnb+1]]<-mere
75 list_fille[[mere]] <- c(clnb,clnb+1)
76 listcol[[clnb]]<-vector()
77 listcol[[clnb+1]]<-vector()
78 #extraction du premier facteur de l'afc
80 pp('taille dtable dans boucle (col/row)',c(ncol(dtable),nrow(dtable)))
81 afc<-boostana(dtable, nd=1, libsvdc=libsvdc, libsvdc.path=libsvdc.path)
82 pp('SV',afc$singular.values)
83 pp('V.P.', afc$eigen.values)
84 coordrow <- as.matrix(afc$row.scores[,1])
86 row.names(coordrow)<-rownames(dtable)
87 coordrow <- cbind(coordrow,1:nrow(dtable))
88 print('deb recherche meilleur partition')
89 ordert <- as.matrix(coordrow[order(coordrow[,1]),])
90 ordert <- cbind(ordert, 1:nrow(ordert))
91 ordert <- ordert[order(ordert[,2]),]
95 dtable <- dtable[order(ordert[,3]),]
99 sc2 <- colSums(dtable) - sc1
100 chitable <- rbind(sc1, sc2)
105 inert <- find.max(dtable, chitable, compte, rmax, maxinter, sc, TT)
106 print('@@@@@@@@@@@@@@@@@@@@@@@@@@@@')
107 pp('max inter phase 1', inert$maxinter/TT)#max(listinter))
108 print('@@@@@@@@@@@@@@@@@@@@@@@@@@@@')
109 ordert <- ordert[order(ordert[,3]),]
110 listclasse<-ifelse(coordrowori<=ordert[(inert$rmax),1],clnb,clnb+1)
111 dtable <- dtable[order(ordert[,2]),]
114 #dtable<-cbind(dtable,'cl'= as.vector(cl))
116 N1<-length(listclasse[listclasse==clnb])
117 N2<-length(listclasse[listclasse==clnb+1])
120 ###################################################################
121 # reclassement des individus #
122 ###################################################################
127 ln <- which(dtable==1, arr.ind=TRUE)
129 lnz[1:nrow(dtable)] <- 0
131 for (k in 1:nrow(ln)) {lnz[[ln[k,1]]]<-append(lnz[[ln[k,1]]],ln[k,2])}
132 for (k in 1:nrow(dtable)) {lnz[[k]] <- lnz[[k]][-1]}
135 while (malcl!=0 & N1>=5 & N2>=5) {
137 listsub[[it]]<-vector()
138 txt <- paste('nombre iteration', it)
139 #pp('nombre iteration',it)
142 t1<-dtable[which(cl[,1]==clnb),]#[,-ncol(dtable)]
143 t2<-dtable[which(cl[,1]==clnb+1),]#[,-ncol(dtable)]
159 chtableori<-rbind(sc1,sc2)
161 interori<-MyChiSq(chtableori,sc,TT)/TT#chisq.test(chtableori)$statistic#/TT
162 txt <- paste(txt, ' - interori : ',interori)
163 #pp('interori',interori)
170 txt <- paste(txt, 'N1:', N1,'-N2:',N2)
176 if(cl[compte]==clnb){
177 chtable[1,l]<-chtable[1,l]-1
178 chtable[2,l]<-chtable[2,l]+1
180 chtable[1,l]<-chtable[1,l]+1
181 chtable[2,l]<-chtable[2,l]-1
183 interswitch<-MyChiSq(chtable,sc,TT)/TT#chisq.test(chtable)$statistic/TT
184 ws<-interori-interswitch
187 interori<-interswitch
188 if(cl[compte]==clnb){
192 listsub[[it]]<-append(listsub[[it]],compte)
197 listsub[[it]]<-append(listsub[[it]],compte)
199 vdelta<-append(vdelta,compte)
204 # for (val in vdelta) {
205 # if (cl[val]==clnb) {
207 # listsub[[it]]<-append(listsub[[it]],val)
210 # listsub[[it]]<-append(listsub[[it]],val)
213 print('###################################')
214 print('longueur < 0')
215 malcl<-length(vdelta)
216 if ((it>1)&&(!is.logical(listsub[[it]]))&&(!is.logical(listsub[[it-1]]))){
217 if (listsub[[it]]==listsub[[(it-1)]]){
222 print('###################################')
224 #dtable<-cbind(dtable,'cl'=as.vector(cl))
225 #dtable[,'cl'] <-as.vector(cl)
226 #######################################################################
227 # Fin reclassement des individus #
228 #######################################################################
229 # if (!(length(cl[cl==clnb])==1 || length(cl[cl==clnb+1])==1)) {
230 #t1<-dtable[dtable[,'cl']==clnb,][,-ncol(dtable)]
231 #t2<-dtable[dtable[,'cl']==clnb+1,][,-ncol(dtable)]
232 t1<-dtable[which(cl[,1]==clnb),]#[,-ncol(dtable)]
233 t2<-dtable[which(cl[,1]==clnb+1),]#[,-ncol(dtable)]
234 if (class(t1)=='numeric') {
241 if (class(t2)=='numeric') {
248 chtable<-rbind(sc1,sc2)
249 inter<-chisq.test(chtable)$statistic/TT
250 pp('last inter',inter)
251 print('=====================')
252 #calcul de la specificite des colonnes
253 mint<-min(nrowt1,nrowt2)
254 maxt<-max(nrowt1,nrowt2)
255 seuil<-round((1.9*(maxt/mint))+1.9,digit=6)
256 #sink('/home/pierre/workspace/iramuteq/dev/findchi2.txt')
257 # print('ATTENTION SEUIL 3,84')
275 #another try#########################################
276 one <- cbind(sc1,sc2)
277 cols <- c(length(which(cl==clnb)), length(which(cl==clnb+1)))
279 colss <- matrix(rep(cols,ncol(dtable)), ncol=2, byrow=TRUE)
281 rows <- cbind(rowSums(zero), rowSums(one))
283 for (m in 1:nrow(rows)) {
284 obs <- t(matrix(c(zero[m,],one[m,]),2,2))
285 E <- outer(rows[m,],cols,'*')/n
286 if ((min(obs[2,])==0) & (min(obs[1,])!=0)) {
288 } else if ((min(obs[1,])==0) & (min(obs[2,])!=0)) {
290 } else if (any(obs < 10)) {
291 chi <- sum((abs(obs - E) - 0.5)^2 / E)
293 chi <- sum(((obs - E)^2)/E)
299 if (obs[2,1] < E[2,1]) {
300 listcol[[clnb]]<-append(listcol[[clnb]],cn[m])
302 } else if (obs[2,2] < E[2,2]) {
303 listcol[[clnb+1]]<-append(listcol[[clnb+1]],cn[m])
308 ######################################################
309 print('resultats elim item')
310 pp(clnb+1,length(listcol[[clnb+1]]))
311 pp(clnb,length(listcol[[clnb]]))
313 pp('ncclnbp',ncclnbp)
314 listrownamedtable<-rownames(dtable)
315 listrownamedtable<-as.integer(listrownamedtable)
316 newcol<-vector(length=nrow(dataori))
317 #remplissage de la nouvelle colonne avec les nouvelles classes
320 newcol[listrownamedtable] <- cl[,1]
321 #recuperation de la classe precedante pour les cases vides
322 print('recuperation classes precedentes')
324 newcol[which(newcol==0)] <- dout[,ncol(dout)][which(newcol==0)]
326 if(!is.null(rowelim)) {
329 tailleclasse<-as.matrix(summary(as.factor(as.character(newcol))))
330 print('tailleclasse')
332 tailleclasse<-as.matrix(tailleclasse[!(rownames(tailleclasse)==0),])
333 plusgrand<-which.max(tailleclasse)
334 #???????????????????????????????????
335 #Si 2 classes ont des effectifs egaux, on prend la premiere de la liste...
336 if (length(plusgrand)>1) {
337 plusgrand<-plusgrand[1]
339 #????????????????????????????????????
341 #constuction du prochain tableau a analyser
342 print('construction tableau suivant')
343 dout<-cbind(dout,newcol)
344 classe<-as.integer(rownames(tailleclasse)[plusgrand])
345 dtable<-dataori[which(newcol==classe),]
346 row.names(dtable)<-rownames(dataori)[which(newcol==classe)]
347 colnames(dtable) <- 1:ncol(dtable)
349 listcolelim<-listcol[[as.integer(classe)]]
350 mother<-listmere[[as.integer(classe)]]
352 listcolelim<-append(listcolelim,listcol[[mother]])
353 mother<-listmere[[mother]]
356 listcolelim<-sort(unique(listcolelim))
357 pp('avant',ncol(dtable))
358 if (!is.logical(listcolelim)){
359 print('elimination colonne')
360 #dtable<-dtable[,-listcolelim]
361 dtable<-dtable[,!(colnames(dtable) %in% listcolelim)]
363 pp('apres',ncol(dtable))
364 #elimination des colonnes ne contenant que des 0
365 print('vire colonne inf 3 dans boucle')
368 dtable<-dtable[,-which(sdt<=3)]
370 #elimination des lignes ne contenant que des 0
371 print('vire ligne vide dans boucle')
372 if (ncol(dtable)==1) {
378 rowelim<-as.integer(rownames(dtable)[which(sdt==0)])
379 print('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&')
381 print('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&')
382 dtable<-dtable[-which(sdt==0),]
387 res <- list(n1 = dout, list_mere = listmere, list_fille = list_fille)