|
0
|
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 infile = sys.argv[1]
|
|
|
13 y_col = int(sys.argv[2])-1
|
|
|
14 x_cols = sys.argv[3].split(',')
|
|
|
15 outfile = sys.argv[4]
|
|
|
16 outfile2 = sys.argv[5]
|
|
|
17 print "Predictor columns: %s; Response column: %d" % ( x_cols, y_col+1 )
|
|
|
18 fout = open(outfile,'w')
|
|
|
19
|
|
|
20 for i, line in enumerate( file ( infile )):
|
|
|
21 line = line.rstrip('\r\n')
|
|
|
22 if len( line )>0 and not line.startswith( '#' ):
|
|
|
23 elems = line.split( '\t' )
|
|
|
24 break
|
|
|
25 if i == 30:
|
|
|
26 break # Hopefully we'll never get here...
|
|
|
27
|
|
|
28 if len( elems )<1:
|
|
|
29 stop_err( "The data in your input dataset is either missing or not formatted properly." )
|
|
|
30
|
|
|
31 y_vals = []
|
|
|
32 x_vals = []
|
|
|
33
|
|
|
34 for k, col in enumerate(x_cols):
|
|
|
35 x_cols[k] = int(col)-1
|
|
|
36 x_vals.append([])
|
|
|
37
|
|
|
38 NA = 'NA'
|
|
|
39 for ind, line in enumerate( file( infile ) ):
|
|
|
40 if line and not line.startswith( '#' ):
|
|
|
41 try:
|
|
|
42 fields = line.split("\t")
|
|
|
43 try:
|
|
|
44 yval = float(fields[y_col])
|
|
|
45 except Exception, ey:
|
|
|
46 yval = r('NA')
|
|
|
47 y_vals.append(yval)
|
|
|
48 for k, col in enumerate(x_cols):
|
|
|
49 try:
|
|
|
50 xval = float(fields[col])
|
|
|
51 except Exception, ex:
|
|
|
52 xval = r('NA')
|
|
|
53 x_vals[k].append(xval)
|
|
|
54 except:
|
|
|
55 pass
|
|
|
56
|
|
|
57 response_term = ""
|
|
|
58
|
|
|
59 x_vals1 = numpy.asarray(x_vals).transpose()
|
|
|
60
|
|
|
61 dat = r.list(x=array(x_vals1), y=y_vals)
|
|
|
62
|
|
|
63 r.library("leaps")
|
|
|
64
|
|
|
65 set_default_mode(NO_CONVERSION)
|
|
|
66 try:
|
|
|
67 leaps = r.regsubsets(r("y ~ x"), data= r.na_exclude(dat))
|
|
|
68 except RException, rex:
|
|
|
69 stop_err("Error performing linear regression on the input data.\nEither the response column or one of the predictor columns contain no numeric values.")
|
|
|
70 set_default_mode(BASIC_CONVERSION)
|
|
|
71
|
|
|
72 summary = r.summary(leaps)
|
|
|
73 tot = len(x_vals)
|
|
|
74 pattern = "["
|
|
|
75 for i in range(tot):
|
|
|
76 pattern = pattern + 'c' + str(int(x_cols[int(i)]) + 1) + ' '
|
|
|
77 pattern = pattern.strip() + ']'
|
|
|
78 print >> fout, "#Vars\t%s\tR-sq\tAdj. R-sq\tC-p\tbic" % (pattern)
|
|
|
79 for ind, item in enumerate(summary['outmat']):
|
|
|
80 print >> fout, "%s\t%s\t%s\t%s\t%s\t%s" % (str(item).count('*'), item, summary['rsq'][ind], summary['adjr2'][ind], summary['cp'][ind], summary['bic'][ind])
|
|
|
81
|
|
|
82
|
|
|
83 r.pdf( outfile2, 8, 8 )
|
|
|
84 r.plot(leaps, scale="Cp", main="Best subsets using Cp Criterion")
|
|
|
85 r.plot(leaps, scale="r2", main="Best subsets using R-sq Criterion")
|
|
|
86 r.plot(leaps, scale="adjr2", main="Best subsets using Adjusted R-sq Criterion")
|
|
|
87 r.plot(leaps, scale="bic", main="Best subsets using bic Criterion")
|
|
|
88
|
|
|
89 r.dev_off()
|
