Mercurial > repos > devteam > rcve

comparison rcve.py @ 0:7740956d197b draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Imported from capsule None
author devteam
date Tue, 01 Apr 2014 10:52:35 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:7740956d197b
1 #!/usr/bin/env python
2
3 import sys
4 from rpy import *
5 import numpy
6
7 def stop_err(msg):
8 sys.stderr.write(msg)
9 sys.exit()
10
11
12 def sscombs(s):
13 if len(s) == 1:
14 return [s]
15 else:
16 ssc = sscombs(s[1:])
17 return [s[0]] + [s[0]+comb for comb in ssc] + ssc
18
19
20 infile = sys.argv[1]
21 y_col = int(sys.argv[2])-1
22 x_cols = sys.argv[3].split(',')
23 outfile = sys.argv[4]
24
25 print "Predictor columns: %s; Response column: %d" % ( x_cols, y_col+1 )
26 fout = open(outfile,'w')
27
28 for i, line in enumerate( file ( infile )):
29 line = line.rstrip('\r\n')
30 if len( line )>0 and not line.startswith( '#' ):
31 elems = line.split( '\t' )
32 break
33 if i == 30:
34 break # Hopefully we'll never get here...
35
36 if len( elems )<1:
37 stop_err( "The data in your input dataset is either missing or not formatted properly." )
38
39 y_vals = []
40 x_vals = []
41
42 for k, col in enumerate(x_cols):
43 x_cols[k] = int(col)-1
44 x_vals.append([])
45 """
46 try:
47 float( elems[x_cols[k]] )
48 except:
49 try:
50 msg = "This operation cannot be performed on non-numeric column %d containing value '%s'." % ( col, elems[x_cols[k]] )
51 except:
52 msg = "This operation cannot be performed on non-numeric data."
53 stop_err( msg )
54 """
55 NA = 'NA'
56 for ind, line in enumerate( file( infile )):
57 if line and not line.startswith( '#' ):
58 try:
59 fields = line.split("\t")
60 try:
61 yval = float(fields[y_col])
62 except Exception, ey:
63 yval = r('NA')
64 #print >>sys.stderr, "ey = %s" %ey
65 y_vals.append(yval)
66 for k, col in enumerate(x_cols):
67 try:
68 xval = float(fields[col])
69 except Exception, ex:
70 xval = r('NA')
71 #print >>sys.stderr, "ex = %s" %ex
72 x_vals[k].append(xval)
73 except:
74 pass
75
76 x_vals1 = numpy.asarray(x_vals).transpose()
77 dat = r.list( x=array(x_vals1), y=y_vals )
78
79 set_default_mode(NO_CONVERSION)
80 try:
81 full = r.lm( r("y ~ x"), data=r.na_exclude(dat) ) #full model includes all the predictor variables specified by the user
82 except RException, rex:
83 stop_err("Error performing linear regression on the input data.\nEither the response column or one of the predictor columns contain no numeric values.")
84 set_default_mode(BASIC_CONVERSION)
85
86 summary = r.summary(full)
87 fullr2 = summary.get('r.squared','NA')
88
89 if fullr2 == 'NA':
90 stop_err("Error in linear regression")
91
92 if len(x_vals) < 10:
93 s = ""
94 for ch in range(len(x_vals)):
95 s += str(ch)
96 else:
97 stop_err("This tool only works with less than 10 predictors.")
98
99 print >> fout, "#Model\tR-sq\tRCVE_Terms\tRCVE_Value"
100 all_combos = sorted(sscombs(s), key=len)
101 all_combos.reverse()
102 for j, cols in enumerate(all_combos):
103 #if len(cols) == len(s): #Same as the full model above
104 # continue
105 if len(cols) == 1:
106 x_vals1 = x_vals[int(cols)]
107 else:
108 x_v = []
109 for col in cols:
110 x_v.append(x_vals[int(col)])
111 x_vals1 = numpy.asarray(x_v).transpose()
112 dat = r.list(x=array(x_vals1), y=y_vals)
113 set_default_mode(NO_CONVERSION)
114 red = r.lm(r("y ~ x"), data= dat) #Reduced model
115 set_default_mode(BASIC_CONVERSION)
116 summary = r.summary(red)
117 redr2 = summary.get('r.squared','NA')
118 try:
119 rcve = (float(fullr2)-float(redr2))/float(fullr2)
120 except:
121 rcve = 'NA'
122 col_str = ""
123 for col in cols:
124 col_str = col_str + str(int(x_cols[int(col)]) + 1) + " "
125 col_str.strip()
126 rcve_col_str = ""
127 for col in s:
128 if col not in cols:
129 rcve_col_str = rcve_col_str + str(int(x_cols[int(col)]) + 1) + " "
130 rcve_col_str.strip()
131 if len(cols) == len(s): #full model
132 rcve_col_str = "-"
133 rcve = "-"
134 try:
135 redr2 = "%.4f" % (float(redr2))
136 except:
137 pass
138 try:
139 rcve = "%.4f" % (float(rcve))
140 except:
141 pass
142 print >> fout, "%s\t%s\t%s\t%s" % ( col_str, redr2, rcve_col_str, rcve )