Mercurial > repos > devteam > best_regression_subsets

comparison best_regression_subsets.py @ 0:2bb8843930ac

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