makenew
[iramuteq] / tabchi2mcnemar.py
1 #!/bin/env python
2 # -*- coding: utf-8 -*-
3 #Author: Pierre Ratinaud
4 #Copyright (c) 2016 Pierre Ratinaud
5 #License: GNU/GPL
6
7 import HTML
8 import os
9 import string
10 import wx
11 import os
12 import sys
13 import tempfile
14 from chemins import ffr,FFF
15 import wx.lib.sized_controls as sc
16 from time import sleep
17 from functions import exec_rcode, check_Rresult
18 from dialog import ChiDialog, PrefChi
19 from analysematrix import AnalyseMatrix
20                     
21 def make_res(line) :
22     if float(line[5]) <= 0.05 and line[6] != 'warning':
23         line.append('green')
24     elif float(line[5]) <= 0.05 and line[6] == 'warning':
25         line.append('blue')
26     else :
27         line.append('red')
28     return line
29
30 def clean_line(result) :
31     return [[val for val in line if val != '**'] for line in result]
32
33 def make_table(tabs, tab_title, res) :
34     return ['<br>'.join(['<font color=%s>%s</font>' % (res[i][-1], tab_title), HTML.table(tab)]) for i,tab in enumerate(tabs)]
35
36 def make_restab(res) :
37     return ['<br>'.join(['<font color=%s>%s</font>'% (line[-1], u'Résultats'),HTML.table([['chi', line[3]],['p', line[5]]])]) for i,line in enumerate(res)]
38
39 def make_htmlgraphs(graphs) :
40     return ['<img src=%s>' % os.path.basename(val) for val in graphs]
41
42 def make_link_list(res, text) :
43     return ['<a name=back_%i></a><a href=#%i><font color=%s>%s</font></a>' % (i, i, chi[-1], text[i]) for i, chi in enumerate(res)]
44
45 def make_title(res, text) :
46     return ['<a name=%i></a><br><font color=%s>%s</font><br><a href=#back_%i>retour</a><br>' % (i, val[-1], text[i], i) for i, val in enumerate(res)] 
47
48
49 chioption = { 'valobs' : True,
50               'valtheo' : True,
51               'resi' : False,
52               'contrib' : True,
53               'pourcent' : False,
54               'pourcentl' : True,
55               'pourcentc' : True,
56               'graph' : True,
57               'bw' : False,
58               }
59
60 class McNemar(AnalyseMatrix):
61     def doparametres(self, dlg = None):
62         if dlg is None :
63             return
64         dial = ChiDialog(self.parent, -1, u"Chi2 McNemar", chioption, self.tableau, size=(400, 350),
65                      style = wx.DEFAULT_DIALOG_STYLE
66                      )
67         dial.CenterOnParent()
68         val = dial.ShowModal()
69         if val==wx.ID_OK :     
70             self.colsel1 = dial.list_box_1.GetSelections()
71             self.colsel2 = dial.list_box_2.GetSelections()
72             if dial.chiopt :
73                 chioption['valobs'] = dial.dial.check1.GetValue()
74                 chioption['valtheo'] = dial.dial.check2.GetValue()
75                 chioption['resi'] = dial.dial.check3.GetValue()
76                 chioption['contrib'] = dial.dial.check4.GetValue()
77                 chioption['pourcent'] = dial.dial.check5.GetValue()
78                 chioption['pourcentl'] = dial.dial.check6.GetValue()
79                 chioption['pourcentc'] = dial.dial.check7.GetValue()
80                 chioption['graph'] = dial.dial.check8.GetValue()
81                 chioption['bw'] = dial.dial.checkbw.GetValue()
82                 dial.dial.Destroy()
83             dial.Destroy()
84             self.parametres.update(chioption)
85             self.chioption = chioption
86         else :
87             if dial.chiopt :
88                 dial.dial.Destroy()
89             dial.Destroy()
90             self.parametres = None
91             
92     def doanalyse(self):
93         self.count = 1
94         keepGoing = self.dlg.Update(self.count,u"Analyse dans R...")
95         self.OutFrame=tempfile.mktemp(dir=self.parent.TEMPDIR)
96         self.encode=self.parent.encode
97         self.TEMPDIR=self.parent.TEMPDIR
98         self.RPath=self.parent.PathPath.get('PATHS','rpath')
99         self.TextCroise=[]
100         for i in self.colsel1 :
101             for j in self.colsel2 :
102                 self.TextCroise.append(self.tableau.colnames[i] + ' / ' + self.tableau.colnames[j])
103         rchioption = {}
104         for val in self.chioption :
105             if self.chioption[val]:
106                 rchioption[val] = 'TRUE'
107             else :
108                 rchioption[val] = 'FALSE'
109         txt="""
110         source("%s")
111         """%ffr(self.parent.RscriptsPath['Rfunct'])
112         txt += """
113         source("%s")
114         """ % ffr(self.parent.RscriptsPath['Rgraph'])
115         txt += """
116         doobs <- %s
117         doexp <- %s
118         docontrib <- %s
119         doresi <- %s
120         dopr <- %s
121         doprl <- %s
122         doprc <- %s
123         dograph <- %s
124         bw <- %s
125         """ % (rchioption['valobs'], rchioption['valtheo'], rchioption['contrib'], rchioption['resi'], rchioption['pourcent'], rchioption['pourcentl'], rchioption['pourcentc'], rchioption['graph'], rchioption['bw'])
126         txt+="""
127         datadm <- read.csv2("%s", encoding="%s", header = TRUE, row.names = 1, sep='\\t', quote = '"', na.string = '')
128         listres<-list()
129         listcol<-list()
130         cont<-1
131         """%(ffr(self.tableau.parametres['csvfile']), self.tableau.parametres['syscoding'])
132         if len(self.colsel1)==1:
133             strsel1=str(self.colsel1).replace(',','')
134         else:
135             strsel1=str(self.colsel1)
136         if len(self.colsel2)==1:
137             strsel2=str(self.colsel2).replace(',','')
138         else:
139             strsel2=str(self.colsel2)
140         txt+="""
141         for (i in c%s) {""" % strsel1
142         txt+="""
143             for (j in c%s) {""" % strsel2
144         txt+="""
145                 tab<-table(datadm[,i+1],datadm[,j+1])
146                 if (min(dim(tab)) != 1) {
147                     chi <- mcnemar.test(tab)
148                     #chi<-chisq.test(tab)
149                     CS<-colSums(tab)
150                     RS<-rowSums(tab)
151                     GT<-sum(tab)
152                     chi$observed <- tab
153                     chi$expected <- tab
154                     chi$contrib <- tab
155                     chi$residuals <- tab
156                     #chi$contrib<-(tab-chi$expected)/sqrt(chi$expected * ((1 - RS/GT) %%*%% t(1 - CS/GT)))
157                     listres[[cont]]<-chi
158                     listcol[[cont]]<-ncol(tab)
159                     cont<-cont+1
160                 } else {
161                     chi <- list(observed = tab, residuals = tab, contrib = tab, statistic = 0, p.value = 1, expected = tab, message = 'pas de calcul')
162                     listres[[cont]] <- chi
163                     listcol[[cont]]<-ncol(tab)
164                     cont <- cont + 1
165                 }
166             }
167         }
168         maxcol<-max(unlist(listcol))+1
169         if (maxcol<7) {maxcol<-7}
170         frameout<-matrix('*',1,maxcol)
171         count<-0
172         for (chi in listres) {
173             if (min(chi$expected)<5) {
174                 att<-"warning"
175             } else {
176                 att<-""
177             }
178             if ('message' %%in%% attributes(chi)$names) {
179                 att <- "Ce chi2 n\'a pas été calculé"
180                 nom_colresi<-colnames(chi$observed)
181                 chi$prl <- chi$expected
182                 chi$prc <- chi$expected
183                 st <- sum(chi$observed)
184             } else {
185                 nom_colresi<-colnames(chi$observed)
186                 st <- sum(chi$observed)
187                 sc <- colSums(chi$observed)
188                 sr <- rowSums(chi$observed)
189                 chi$prl <- round((chi$observed/sr)*100,2)
190                 chi$prc <- t(round((t(chi$observed)/sc)*100,2))
191             }
192             fileout<-paste('histo_',count,sep='')
193             fileout<-paste(fileout,'.png',sep='')
194             count<-count+1
195             fileout<-file.path("%s",fileout)
196             if (max(nchar(colnames(chi$observed)))>15) {
197                 leg <- 1:length(colnames(chi$observed))
198             } else {
199                 leg <- colnames(chi$observed)
200             }
201             if (dograph) {
202                 width<-ncol(chi$observed)*100
203                 if (width < 350) {width <- 350}
204                 open_file_graph(fileout,width = width, height = 300)
205                 par(mar=c(0,0,0,0))
206                 layout(matrix(c(1,2),1,2, byrow=TRUE),widths=c(3,1))
207                 par(mar=c(2,2,1,0))
208                 par(cex=0.8)
209                 if (!bw) colors <- rainbow(length(rownames(chi$observed)))
210                 else colors <- gray.colors(length(rownames(chi$observed)))
211                 barplot(chi$prl,names.arg = leg, beside=TRUE,border=NA, col=colors)
212                 par(mar=c(0,0,0,0))
213                 par(cex=0.8)
214                 plot(0, axes = FALSE, pch = '')
215                 legend(x = 'center' , rownames(chi$observed), fill = colors)
216                 dev.off()
217             }
218             chi$prl <- cbind(chi$prl, total = rowSums(chi$prl))
219             chi$prc <- rbind(chi$prc, total = colSums(chi$prc))
220             chi$observed<-rbind(chi$observed,total=colSums(chi$observed))
221             chi$observed<-cbind(chi$observed,total=rowSums(chi$observed))
222             chi$pr <- round((chi$observed/st)*100,2)
223             chi$expected<-rbind(chi$expected,total=colSums(chi$expected))
224             chi$expected<-cbind(chi$expected,total=rowSums(chi$expected))
225             chi$expected<-round(chi$expected,digits=2)
226             chi$residuals<-round(chi$residuals,digits=2)
227             chi$contrib<-round(chi$contrib, digits=2)
228             nom_col<-colnames(chi$observed)
229            
230             if (ncol(chi$observed)<maxcol) {
231                 for (i in 1:(maxcol-ncol(chi$observed))) {
232                     chi$observed<-cbind(chi$observed,'**')
233                     chi$pr<-cbind(chi$pr,'**')
234                     chi$prl<-cbind(chi$prl,'**')
235                     chi$prc<-cbind(chi$prc,'**')
236                     chi$expected<-cbind(chi$expected,'**')
237                     chi$residuals<-cbind(chi$residuals,'**')
238                     chi$contrib<-cbind(chi$contrib,'**')
239                     nom_col<-append(nom_col,'**')
240                     nom_colresi<-append(nom_colresi,'**')
241                 }
242                 chi$residuals<-cbind(chi$residuals,'**')
243                 chi$contrib<-cbind(chi$contrib,'**')
244                 nom_colresi<-append(nom_colresi,'**')
245                 chi$prc<-cbind(chi$prc,'**')
246             } else if (ncol(chi$observed)==maxcol) {
247                 chi$residuals<-cbind(chi$residuals,'**')
248                 chi$contrib<-cbind(chi$contrib,'**')
249                 nom_colresi<-append(nom_colresi,'**')
250                 chi$prc<-cbind(chi$prc,'**')
251             }
252             if (doobs) {
253                 li<-matrix('*obs*',1,maxcol)
254                 frameout<-rbind(frameout,li)
255                 frameout<-rbind(frameout,nom_col)
256                 frameout<-rbind(frameout,chi$observed)
257             }
258             if (doexp) {
259                 li<-matrix('*exp*',1,maxcol)
260                 frameout<-rbind(frameout,li)
261                 frameout<-rbind(frameout,nom_col)
262                 frameout<-rbind(frameout,chi$expected)
263             }
264             if (doresi) {
265                 li<-matrix('*resi*',1,maxcol)
266                 frameout<-rbind(frameout,li)
267                 frameout<-rbind(frameout,nom_colresi)
268                 frameout<-rbind(frameout,chi$residuals)
269             }
270             if (docontrib) {
271                 li<-matrix('*contrib*',1,maxcol)
272                 frameout<-rbind(frameout,li)
273                 frameout<-rbind(frameout,nom_colresi)
274                 frameout<-rbind(frameout,chi$contrib)
275             }
276             if (dopr) {
277                 li<-matrix('*pr*', 1, maxcol)
278                 frameout<-rbind(frameout,li)
279                 frameout<-rbind(frameout,nom_col)
280                 frameout<-rbind(frameout,chi$pr)
281             }
282             if (doprl) {
283                 li<-matrix('*prl*', 1, maxcol)
284                 frameout<-rbind(frameout,li)
285                 frameout<-rbind(frameout,nom_col)
286                 frameout<-rbind(frameout,chi$prl)
287             }
288             if (doprc) {
289                 li<-matrix('*prc*', 1, maxcol)
290                 frameout<-rbind(frameout,li)
291                 frameout<-rbind(frameout,nom_colresi)
292                 frameout<-rbind(frameout,chi$prc)
293             }
294             res<-c('****','chi',chi$statistic,'p',chi$p.value,att,fileout)
295             frameout<-rbind(frameout,res)
296         }
297         li<-matrix('fin_analyse',1,maxcol)
298         frameout<-rbind(frameout,li)
299         write.csv2(frameout,file="%s")
300         """ % (ffr(self.parametres['pathout']),ffr(self.OutFrame))
301         tmpfile=tempfile.mktemp(dir=self.TEMPDIR)
302         print tmpfile
303         tmpscript=open(tmpfile,'w')
304         tmpscript.write(txt)
305         tmpscript.close()
306         pid = exec_rcode(self.RPath, tmpfile, wait = False)
307         while pid.poll() == None :
308             sleep(0.2)
309         check_Rresult(self.parent, pid)            
310         self.count += 1
311         keepGoing = self.dlg.Update(self.count,u"Ecriture des résultats")
312          
313         listfileout = self.dolayout(self.chioption)
314             #listfileout=dlg.ShowChi2(ColSel1,ColSel2)
315             #parent.FreqNum += 1
316             #parent.DictTab[u"Chi2_%s*"%parent.FreqNum]=listfileout
317 #             parent.newtab = wx.html.HtmlWindow(parent.nb, -1)
318 #             if "gtk2" in wx.PlatformInfo:
319 #                 parent.newtab.SetStandardFonts()
320 #             parent.newtab.LoadPage(listfileout[len(listfileout)-1])
321 #             parent.nb.AddPage(parent.newtab,u"Chi2_%s*"%parent.FreqNum)
322 #             parent.nb.SetSelection(parent.nb.GetPageCount()-1)
323 #             parent.ShowTab(wx.EVT_BUTTON)
324 #             parent.DisEnSaveTabAs(True)
325 #        self.count += 1
326 #        keepGoing = self.dlg.Update(self.count,u"Fini")
327
328     def dolayout(self, option):
329         ListFile=[False]
330         file=open(self.OutFrame,'rU')
331         content=file.readlines()
332         file.close()
333         lcont = [line.replace('"','').replace('\n','').split(';') for line in content]
334     
335         lcont.pop(0)
336         lcont.pop(0)
337         
338         allcoord = []
339         names = []
340
341         res = [chi for chi in lcont if chi[0]=='res']
342         res = [make_res(line) for line in res]
343         coord_res = [i for i,chi in enumerate(lcont) if chi[0]=='res']
344         if option['valobs']:
345             allcoord.append([i for i,chi in enumerate(lcont) if chi[1]=='*obs*'])
346             names.append(u'Valeurs observées')
347         if option['valtheo'] :
348             allcoord.append([i for i,chi in enumerate(lcont) if chi[1]=='*exp*'])
349             names.append(u'Valeurs théoriques')
350         if option['resi'] :
351             allcoord.append([i for i,chi in enumerate(lcont) if chi[1]=='*resi*'])
352             names.append(u'Residuals')
353         if option['contrib'] :
354             allcoord.append([i for i,chi in enumerate(lcont) if chi[1]=='*contrib*'])
355             names.append(u'Contributions a posteriori')
356         if option['pourcent'] : 
357             allcoord.append([i for i,chi in enumerate(lcont) if chi[1]=='*pr*'])
358             names.append(u'Pourcentages')
359         if option['pourcentl'] :
360             allcoord.append([i for i,chi in enumerate(lcont) if chi[1]=='*prl*'])
361             names.append(u'Pourcentages en ligne')
362         if option['pourcentc'] :
363             allcoord.append([i for i,chi in enumerate(lcont) if chi[1]=='*prc*'])
364             names.append(u'Pourcentages en colonne')
365     
366         allcoord.append(coord_res)
367         allhtml = [[clean_line(lcont[allcoord[i][j]+1:allcoord[i+1][j]]) for j, line in enumerate(allcoord[i])] for i, tab in enumerate(allcoord) if i!=len(allcoord)-1]
368         allhtml = [make_table(val,names[i],res) for i,val in enumerate(allhtml)]
369         links = make_link_list(res, self.TextCroise)
370         html_res = make_restab(res)
371         allhtml.insert(0,html_res)
372         titles = make_title(res, self.TextCroise)
373         allhtml.insert(0,titles)
374         
375         if option['graph'] :
376             graphs = [line[7] for line in res]
377             ListFile += graphs
378             html_graphs = make_htmlgraphs(graphs)
379             allhtml.append(html_graphs)
380
381         header=u"""
382         <html>\n
383         <meta http-equiv="content-Type" content="text/html; charset=%s" />\n
384         <body>\n
385         <h1>Test du Chi2</h1>\n
386         <br>
387         <table border=1><tr><td>
388         Légende : <br>
389         <font color=green>p &lt;= 0.05</font><br>
390         <font color=blue>p &lt;= 0.05 mais il y a des valeurs théoriques &lt; 5</font><br>
391         <font color=red>p &gt; 0.05</font>
392         </td></tr></table><br><br>
393         """%self.parent.SysEncoding
394
395         pretxt = '<br>\n'.join(links)+'<br><hr><br>\n'
396         txt = '<br><hr><br>\n'.join(['<br><br>'.join([tab[i] for tab in allhtml]) for i,val in enumerate(res)])
397         txt = header + pretxt + txt + '\n</body></html>'
398
399         fileout=os.path.join(self.parametres['pathout'],'resultats-chi2.html')
400         with open(fileout, 'w') as f :
401             f.write(txt)
402         ListFile.append(fileout)         
403         return ListFile