# HG changeset patch
# User fubar
# Date 1446269667 14400
# Node ID 21b12c7c52e4a30ba905b684b64adc79b34246c1
# Parent 7d66bfa4fd56d8a6a748ca9d973e5ff913dced6e
Fixes to paths in git for deps
diff -r 7d66bfa4fd56 -r 21b12c7c52e4 .shed.yml
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/.shed.yml Sat Oct 31 01:34:27 2015 -0400
@@ -0,0 +1,9 @@
+categories:
+- Statistics
+description: glmnet and lars for galaxy
+long_description: |
+ Allows hybrid lasso-ridge regression using glmnet
+name: rglasso_1_9_8
+owner: fubar
+remote_repository_url: https://github.com/galaxyproject/tools-iuc/tree/master/tools/rglasso
+type: unrestricted
diff -r 7d66bfa4fd56 -r 21b12c7c52e4 readme.rst
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/readme.rst Sat Oct 31 01:34:27 2015 -0400
@@ -0,0 +1,39 @@
+glmnet wrappers
+===============
+
+This is a self installing Galaxy tool exposing the glmnet_ R package which has excellent documentation at
+glmnet_ Minimal details are provided in this wrapper - please RTM to get the best out of it.
+
+The tool exposes the entire range of penalised maximum likelihood
+GLM models ranging from pure lasso (set alpha to 1) to pure ridge-regression (set alpha to 0).
+
+These models can be k-fold internally cross validated to help select an "optimal" predictive or classification
+algorithm. Predictive coefficients for each included independent variable are output for each model.
+
+Predictors can be forced into models to adjust for known confounders or explanatory factors.
+
+The glmnet_ implementation of the coordinate descent algorithm is fast and efficient even on relatively large problems
+with tens of thousands of predictors and thousands of samples - such as normalised microarray intensities and anthropometry
+on a very large sample of obese patients.
+
+The user supplies a tabular file with rows as samples and columns containing observations, then chooses
+as many predictors as required. A separate model will be output for each of potentially multiple dependent
+variables. Models are reported as the coefficients for terms in an 'optimal' model.
+These optimal predictors are selected by repeatedly setting
+aside a random subsample, building a model in the remainder and estimating AUC or deviance
+using k (default 10) fold internal cross validation. For each of these steps, a random 1/k
+of the samples are set aside and used to estiamte performance of an optimal model estimated
+from the remaining samples. Plots are provided showing the range of these (eg 10) internal validation
+estimates and mean model AUC (binomial) or residual deviance plots at each penalty increment step.
+
+A full range of link functions are available including Gaussian, Poisson, Binomial and
+Cox proportional hazard time to failure for censored data in this wrapper.
+
+Note that multinomial and multiresponse gaussian models are NOT yet implemented since I have not yet
+had use for them - send code!
+
+.. _glmnet: http://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html
+
+Wrapper author: Ross Lazarus
+19 october 2014
+
diff -r 7d66bfa4fd56 -r 21b12c7c52e4 rgToolFactory.py
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/rgToolFactory.py Sat Oct 31 01:34:27 2015 -0400
@@ -0,0 +1,649 @@
+# rgToolFactory.py
+# see https://bitbucket.org/fubar/galaxytoolfactory/wiki/Home
+#
+# copyright ross lazarus (ross stop lazarus at gmail stop com) May 2012
+#
+# all rights reserved
+# Licensed under the LGPL
+# suggestions for improvement and bug fixes welcome at https://bitbucket.org/fubar/galaxytoolfactory/wiki/Home
+#
+# july 2014
+# added buffered read of sterror after run
+#
+# august 2013
+# found a problem with GS if $TMP or $TEMP missing - now inject /tmp and warn
+#
+# july 2013
+# added ability to combine images and individual log files into html output
+# just make sure there's a log file foo.log and it will be output
+# together with all images named like "foo_*.pdf
+# otherwise old format for html
+#
+# January 2013
+# problem pointed out by Carlos Borroto
+# added escaping for <>$ - thought I did that ages ago...
+#
+# August 11 2012
+# changed to use shell=False and cl as a sequence
+
+# This is a Galaxy tool factory for simple scripts in python, R or whatever ails ye.
+# It also serves as the wrapper for the new tool.
+#
+# you paste and run your script
+# Only works for simple scripts that read one input from the history.
+# Optionally can write one new history dataset,
+# and optionally collect any number of outputs into links on an autogenerated HTML page.
+
+# DO NOT install on a public or important site - please.
+
+# installed generated tools are fine if the script is safe.
+# They just run normally and their user cannot do anything unusually insecure
+# but please, practice safe toolshed.
+# Read the fucking code before you install any tool
+# especially this one
+
+# After you get the script working on some test data, you can
+# optionally generate a toolshed compatible gzip file
+# containing your script safely wrapped as an ordinary Galaxy script in your local toolshed for
+# safe and largely automated installation in a production Galaxy.
+
+# If you opt for an HTML output, you get all the script outputs arranged
+# as a single Html history item - all output files are linked, thumbnails for all the pdfs.
+# Ugly but really inexpensive.
+#
+# Patches appreciated please.
+#
+#
+# long route to June 2012 product
+# Behold the awesome power of Galaxy and the toolshed with the tool factory to bind them
+# derived from an integrated script model
+# called rgBaseScriptWrapper.py
+# Note to the unwary:
+# This tool allows arbitrary scripting on your Galaxy as the Galaxy user
+# There is nothing stopping a malicious user doing whatever they choose
+# Extremely dangerous!!
+# Totally insecure. So, trusted users only
+#
+# preferred model is a developer using their throw away workstation instance - ie a private site.
+# no real risk. The universe_wsgi.ini admin_users string is checked - only admin users are permitted to run this tool.
+#
+
+import sys
+import shutil
+import subprocess
+import os
+import time
+import tempfile
+import optparse
+import tarfile
+import re
+import shutil
+import math
+
+progname = os.path.split(sys.argv[0])[1]
+myversion = 'V000.2 June 2012'
+verbose = False
+debug = False
+toolFactoryURL = 'https://bitbucket.org/fubar/galaxytoolfactory'
+buffsize = 1048576
+
+
+def timenow():
+ """return current time as a string
+ """
+ return time.strftime('%d/%m/%Y %H:%M:%S', time.localtime(time.time()))
+
+html_escape_table = {
+ "&": "&",
+ ">": ">",
+ "<": "<",
+ "$": "\$"
+ }
+
+def html_escape(text):
+ """Produce entities within text."""
+ return "".join(html_escape_table.get(c,c) for c in text)
+
+def cmd_exists(cmd):
+ return subprocess.call("type " + cmd, shell=True,
+ stdout=subprocess.PIPE, stderr=subprocess.PIPE) == 0
+
+
+class ScriptRunner:
+ """class is a wrapper for an arbitrary script
+ """
+
+ def __init__(self,opts=None):
+ """
+ cleanup inputs, setup some outputs
+
+ """
+ self.useGM = cmd_exists('gm')
+ self.useIM = cmd_exists('convert')
+ self.useGS = cmd_exists('gs')
+ self.temp_warned = False # we want only one warning if $TMP not set
+ if opts.output_dir: # simplify for the tool tarball
+ os.chdir(opts.output_dir)
+ self.thumbformat = 'png'
+ self.opts = opts
+ self.toolname = re.sub('[^a-zA-Z0-9_]+', '', opts.tool_name) # a sanitizer now does this but..
+ self.toolid = self.toolname
+ self.myname = sys.argv[0] # get our name because we write ourselves out as a tool later
+ self.pyfile = self.myname # crude but efficient - the cruft won't hurt much
+ self.xmlfile = '%s.xml' % self.toolname
+ s = open(self.opts.script_path,'r').readlines()
+ s = [x.rstrip() for x in s] # remove pesky dos line endings if needed
+ self.script = '\n'.join(s)
+ fhandle,self.sfile = tempfile.mkstemp(prefix=self.toolname,suffix=".%s" % (opts.interpreter))
+ tscript = open(self.sfile,'w') # use self.sfile as script source for Popen
+ tscript.write(self.script)
+ tscript.close()
+ self.indentedScript = '\n'.join([' %s' % x for x in s]) # for restructured text in help
+ self.escapedScript = '\n'.join([html_escape(x) for x in s])
+ self.elog = os.path.join(self.opts.output_dir,"%s_error.log" % self.toolname)
+ self.tlog = os.path.join(self.opts.output_dir,"%s_runner.log" % self.toolname)
+ if opts.output_dir: # may not want these complexities
+ art = '%s.%s' % (self.toolname,opts.interpreter)
+ artpath = os.path.join(self.opts.output_dir,art) # need full path
+ artifact = open(artpath,'w') # use self.sfile as script source for Popen
+ artifact.write(self.script)
+ artifact.close()
+ self.html = []
+ self.cl = (opts.interpreter,self.sfile)
+ self.outFormats = 'tabular' # TODO make this an option at tool generation time
+ self.inputFormats = 'tabular' # TODO make this an option at tool generation time
+ self.test1Input = '%s_test1_input.xls' % self.toolname
+ self.test1Output = '%s_test1_output.xls' % self.toolname
+ self.test1HTML = '%s_test1_output.html' % self.toolname
+
+ def makeXML(self):
+ """
+ Create a Galaxy xml tool wrapper for the new script as a string to write out
+ fixme - use templating or something less fugly than this example of what we produce
+
+
+ a tabular file
+
+ reverse.py --script_path "$runMe" --interpreter "python"
+ --tool_name "reverse" --input_tab "$input1" --output_tab "$tab_file"
+
+
+
+
+
+
+
+
+
+
+
+**What it Does**
+
+Reverse the columns in a tabular file
+
+
+
+
+
+# reverse order of columns in a tabular file
+import sys
+inp = sys.argv[1]
+outp = sys.argv[2]
+i = open(inp,'r')
+o = open(outp,'w')
+for row in i:
+ rs = row.rstrip().split('\t')
+ rs.reverse()
+ o.write('\t'.join(rs))
+ o.write('\n')
+i.close()
+o.close()
+
+
+
+
+
+
+ """
+ newXML="""
+ %(tooldesc)s
+ %(command)s
+
+ %(inputs)s
+
+
+ %(outputs)s
+
+
+
+ %(script)s
+
+
+ %(tooltests)s
+
+ %(help)s
+
+ """ # needs a dict with toolname, toolid, interpreter, scriptname, command, inputs as a multi line string ready to write, outputs ditto, help ditto
+
+ newCommand="""
+ %(toolname)s.py --script_path "$runMe" --interpreter "%(interpreter)s"
+ --tool_name "%(toolname)s" %(command_inputs)s %(command_outputs)s
+ """ # may NOT be an input or htmlout
+ tooltestsTabOnly = """
+
+
+
+
+ """
+ tooltestsHTMLOnly = """
+
+
+
+
+ """
+ tooltestsBoth = """
+
+
+
+
+
+ """
+ xdict = {}
+ xdict['tool_version'] = self.opts.tool_version
+ xdict['test1Input'] = self.test1Input
+ xdict['test1HTML'] = self.test1HTML
+ xdict['test1Output'] = self.test1Output
+ if self.opts.make_HTML and self.opts.output_tab <> 'None':
+ xdict['tooltests'] = tooltestsBoth % xdict
+ elif self.opts.make_HTML:
+ xdict['tooltests'] = tooltestsHTMLOnly % xdict
+ else:
+ xdict['tooltests'] = tooltestsTabOnly % xdict
+ xdict['script'] = self.escapedScript
+ # configfile is least painful way to embed script to avoid external dependencies
+ # but requires escaping of <, > and $ to avoid Mako parsing
+ if self.opts.help_text:
+ xdict['help'] = open(self.opts.help_text,'r').read()
+ else:
+ xdict['help'] = 'Please ask the tool author for help as none was supplied at tool generation'
+ coda = ['**Script**','Pressing execute will run the following code over your input file and generate some outputs in your history::']
+ coda.append(self.indentedScript)
+ coda.append('**Attribution** This Galaxy tool was created by %s at %s\nusing the Galaxy Tool Factory.' % (self.opts.user_email,timenow()))
+ coda.append('See %s for details of that project' % (toolFactoryURL))
+ coda.append('Please cite: Creating re-usable tools from scripts: The Galaxy Tool Factory. Ross Lazarus; Antony Kaspi; Mark Ziemann; The Galaxy Team. ')
+ coda.append('Bioinformatics 2012; doi: 10.1093/bioinformatics/bts573')
+ xdict['help'] = '%s\n%s' % (xdict['help'],'\n'.join(coda))
+ if self.opts.tool_desc:
+ xdict['tooldesc'] = '%s' % self.opts.tool_desc
+ else:
+ xdict['tooldesc'] = ''
+ xdict['command_outputs'] = ''
+ xdict['outputs'] = ''
+ if self.opts.input_tab <> 'None':
+ xdict['command_inputs'] = '--input_tab "$input1" ' # the space may matter a lot if we append something
+ xdict['inputs'] = ' \n' % self.inputFormats
+ else:
+ xdict['command_inputs'] = '' # assume no input - eg a random data generator
+ xdict['inputs'] = ''
+ xdict['inputs'] += ' \n' % self.toolname
+ xdict['toolname'] = self.toolname
+ xdict['toolid'] = self.toolid
+ xdict['interpreter'] = self.opts.interpreter
+ xdict['scriptname'] = self.sfile
+ if self.opts.make_HTML:
+ xdict['command_outputs'] += ' --output_dir "$html_file.files_path" --output_html "$html_file" --make_HTML "yes" '
+ xdict['outputs'] += ' \n'
+ if self.opts.output_tab <> 'None':
+ xdict['command_outputs'] += ' --output_tab "$tab_file"'
+ xdict['outputs'] += ' \n' % self.outFormats
+ xdict['command'] = newCommand % xdict
+ xmls = newXML % xdict
+ xf = open(self.xmlfile,'w')
+ xf.write(xmls)
+ xf.write('\n')
+ xf.close()
+ # ready for the tarball
+
+
+ def makeTooltar(self):
+ """
+ a tool is a gz tarball with eg
+ /toolname/tool.xml /toolname/tool.py /toolname/test-data/test1_in.foo ...
+ """
+ retval = self.run()
+ if retval:
+ print >> sys.stderr,'## Run failed. Cannot build yet. Please fix and retry'
+ sys.exit(1)
+ self.makeXML()
+ tdir = self.toolname
+ os.mkdir(tdir)
+ if self.opts.input_tab <> 'None': # no reproducible test otherwise? TODO: maybe..
+ testdir = os.path.join(tdir,'test-data')
+ os.mkdir(testdir) # make tests directory
+ shutil.copyfile(self.opts.input_tab,os.path.join(testdir,self.test1Input))
+ if self.opts.output_tab <> 'None':
+ shutil.copyfile(self.opts.output_tab,os.path.join(testdir,self.test1Output))
+ if self.opts.make_HTML:
+ shutil.copyfile(self.opts.output_html,os.path.join(testdir,self.test1HTML))
+ if self.opts.output_dir:
+ shutil.copyfile(self.tlog,os.path.join(testdir,'test1_out.log'))
+ op = '%s.py' % self.toolname # new name
+ outpiname = os.path.join(tdir,op) # path for the tool tarball
+ pyin = os.path.basename(self.pyfile) # our name - we rewrite ourselves (TM)
+ notes = ['# %s - a self annotated version of %s generated by running %s\n' % (op,pyin,pyin),]
+ notes.append('# to make a new Galaxy tool called %s\n' % self.toolname)
+ notes.append('# User %s at %s\n' % (self.opts.user_email,timenow()))
+ pi = open(self.pyfile,'r').readlines() # our code becomes new tool wrapper (!) - first Galaxy worm
+ notes += pi
+ outpi = open(outpiname,'w')
+ outpi.write(''.join(notes))
+ outpi.write('\n')
+ outpi.close()
+ stname = os.path.join(tdir,self.sfile)
+ if not os.path.exists(stname):
+ shutil.copyfile(self.sfile, stname)
+ xtname = os.path.join(tdir,self.xmlfile)
+ if not os.path.exists(xtname):
+ shutil.copyfile(self.xmlfile,xtname)
+ tarpath = "%s.gz" % self.toolname
+ tar = tarfile.open(tarpath, "w:gz")
+ tar.add(tdir,arcname=self.toolname)
+ tar.close()
+ shutil.copyfile(tarpath,self.opts.new_tool)
+ shutil.rmtree(tdir)
+ ## TODO: replace with optional direct upload to local toolshed?
+ return retval
+
+
+ def compressPDF(self,inpdf=None,thumbformat='png'):
+ """need absolute path to pdf
+ note that GS gets confoozled if no $TMP or $TEMP
+ so we set it
+ """
+ assert os.path.isfile(inpdf), "## Input %s supplied to %s compressPDF not found" % (inpdf,self.myName)
+ our_env = os.environ.copy()
+ if not (our_env.get('TMP',None) or our_env.get('TEMP',None)):
+ our_env['TMP'] = '/tmp'
+ if not self.temp_warned:
+ print >> sys.stdout,'## WARNING - no $TMP or $TEMP!!! Please fix - using /tmp temporarily'
+ self.temp_warned = True
+ hlog = os.path.join(self.opts.output_dir,"compress_%s.txt" % os.path.basename(inpdf))
+ sto = open(hlog,'w')
+ outpdf = '%s_compressed' % inpdf
+ cl = ["gs", "-sDEVICE=pdfwrite", "-dNOPAUSE", "-dUseCIEColor", "-dBATCH","-dPDFSETTINGS=/printer", "-sOutputFile=%s" % outpdf,inpdf]
+ x = subprocess.Popen(cl,stdout=sto,stderr=sto,cwd=self.opts.output_dir,env=our_env)
+ retval1 = x.wait()
+ sto.close()
+ if retval1 == 0:
+ os.unlink(inpdf)
+ shutil.move(outpdf,inpdf)
+ os.unlink(hlog)
+ else:
+ x = open(hlog,'r').readlines()
+ print >> sys.stdout,x
+ hlog = os.path.join(self.opts.output_dir,"thumbnail_%s.txt" % os.path.basename(inpdf))
+ sto = open(hlog,'w')
+ outpng = '%s.%s' % (os.path.splitext(inpdf)[0],thumbformat)
+ if self.useGM:
+ cl2 = ['gm', 'convert', inpdf, outpng]
+ else: # assume imagemagick
+ cl2 = ['convert', inpdf, outpng]
+ x = subprocess.Popen(cl2,stdout=sto,stderr=sto,cwd=self.opts.output_dir,env=our_env)
+ retval2 = x.wait()
+ sto.close()
+ if retval2 <> 0:
+ x = open(hlog,'r').readlines()
+ print >> sys.stdout,x
+ else:
+ os.unlink(hlog)
+ retval = retval1 or retval2
+ return retval
+
+
+ def getfSize(self,fpath,outpath):
+ """
+ format a nice file size string
+ """
+ size = ''
+ fp = os.path.join(outpath,fpath)
+ if os.path.isfile(fp):
+ size = '0 B'
+ n = float(os.path.getsize(fp))
+ if n > 2**20:
+ size = '%1.1f MB' % (n/2**20)
+ elif n > 2**10:
+ size = '%1.1f KB' % (n/2**10)
+ elif n > 0:
+ size = '%d B' % (int(n))
+ return size
+
+ def makeHtml(self):
+ """ Create an HTML file content to list all the artifacts found in the output_dir
+ """
+
+ galhtmlprefix = """
+
+
\n"""
+
+ flist = os.listdir(self.opts.output_dir)
+ flist = [x for x in flist if x <> 'Rplots.pdf']
+ flist.sort()
+ html = []
+ html.append(galhtmlprefix % progname)
+ html.append('
Galaxy Tool "%s" run at %s
' % (self.toolname,timenow()))
+ fhtml = []
+ if len(flist) > 0:
+ logfiles = [x for x in flist if x.lower().endswith('.log')] # log file names determine sections
+ logfiles.sort()
+ logfiles = [x for x in logfiles if os.path.abspath(x) <> os.path.abspath(self.tlog)]
+ logfiles.append(os.path.abspath(self.tlog)) # make it the last one
+ pdflist = []
+ npdf = len([x for x in flist if os.path.splitext(x)[-1].lower() == '.pdf'])
+ for rownum,fname in enumerate(flist):
+ dname,e = os.path.splitext(fname)
+ sfsize = self.getfSize(fname,self.opts.output_dir)
+ if e.lower() == '.pdf' : # compress and make a thumbnail
+ thumb = '%s.%s' % (dname,self.thumbformat)
+ pdff = os.path.join(self.opts.output_dir,fname)
+ retval = self.compressPDF(inpdf=pdff,thumbformat=self.thumbformat)
+ if retval == 0:
+ pdflist.append((fname,thumb))
+ else:
+ pdflist.append((fname,fname))
+ if (rownum+1) % 2 == 0:
+ fhtml.append('
' % (fname,fname,sfsize))
+ for logfname in logfiles: # expect at least tlog - if more
+ if os.path.abspath(logfname) == os.path.abspath(self.tlog): # handled later
+ sectionname = 'All tool run'
+ if (len(logfiles) > 1):
+ sectionname = 'Other'
+ ourpdfs = pdflist
+ else:
+ realname = os.path.basename(logfname)
+ sectionname = os.path.splitext(realname)[0].split('_')[0] # break in case _ added to log
+ ourpdfs = [x for x in pdflist if os.path.basename(x[0]).split('_')[0] == sectionname]
+ pdflist = [x for x in pdflist if os.path.basename(x[0]).split('_')[0] <> sectionname] # remove
+ nacross = 1
+ npdf = len(ourpdfs)
+
+ if npdf > 0:
+ nacross = math.sqrt(npdf) ## int(round(math.log(npdf,2)))
+ if int(nacross)**2 != npdf:
+ nacross += 1
+ nacross = int(nacross)
+ width = min(400,int(1200/nacross))
+ html.append('
%s images and outputs
' % sectionname)
+ html.append('(Click on a thumbnail image to download the corresponding original PDF image) ')
+ ntogo = nacross # counter for table row padding with empty cells
+ html.append('
\n
')
+ for i,paths in enumerate(ourpdfs):
+ fname,thumb = paths
+ s= """
\n""" % (fname,thumb,fname,width,fname)
+ if ((i+1) % nacross == 0):
+ s += '
\n'
+ ntogo = 0
+ if i < (npdf - 1): # more to come
+ s += '
\n')
+ else:
+ if ntogo > 0: # pad
+ html.append('
'*ntogo)
+ html.append('\n')
+ logt = open(logfname,'r').readlines()
+ logtext = [x for x in logt if x.strip() > '']
+ html.append('
%s log output
' % sectionname)
+ if len(logtext) > 1:
+ html.append('\n
\n')
+ html += logtext
+ html.append('\n
\n')
+ else:
+ html.append('%s is empty ' % logfname)
+ if len(fhtml) > 0:
+ fhtml.insert(0,'
Output File Name (click to view)
Size
\n')
+ fhtml.append('
')
+ html.append('
All output files available for downloading
\n')
+ html += fhtml # add all non-pdf files to the end of the display
+ else:
+ html.append('
### Error - %s returned no files - please confirm that parameters are sane
' % self.opts.interpreter)
+ html.append(galhtmlpostfix)
+ htmlf = file(self.opts.output_html,'w')
+ htmlf.write('\n'.join(html))
+ htmlf.write('\n')
+ htmlf.close()
+ self.html = html
+
+ def run(self):
+ """
+ scripts must be small enough not to fill the pipe!
+ """
+ my_env = os.environ.copy()
+ if self.opts.output_dir:
+ ste = open(self.elog,'w')
+ sto = open(self.tlog,'w')
+ sto.write('## Toolfactory running %s as %s script\n' % (self.toolname,self.opts.interpreter))
+ sto.flush()
+ p = subprocess.Popen(self.cl,shell=False,stdout=sto,stderr=ste,cwd=self.opts.output_dir,env=my_env)
+ retval = p.wait()
+ sto.close()
+ ste.close()
+ # get stderr, allowing for case where it's very large
+ tmp_stderr = open( self.elog, 'rb' )
+ stderr = ''
+ try:
+ while True:
+ stderr += tmp_stderr.read( buffsize )
+ if not stderr or len( stderr ) % buffsize != 0:
+ break
+ except OverflowError:
+ pass
+ tmp_stderr.close()
+ else:
+ p = subprocess.Popen(self.cl,shell=False,env=my_env)
+ retval = p.wait()
+ if self.opts.make_HTML:
+ self.makeHtml()
+ return retval
+
+
+ def remove_me_runBash(self):
+ """
+ cannot use - for bash so use self.sfile
+ """
+ if self.opts.output_dir:
+ s = '## Toolfactory generated command line = %s\n' % ' '.join(self.cl)
+ ste = open(self.elog,'w')
+ sto = open(self.tlog,'w')
+ sto.write(s)
+ sto.flush()
+ p = subprocess.Popen(self.cl,shell=False,stdout=sto,stderr=ste,cwd=self.opts.output_dir)
+ else:
+ p = subprocess.Popen(self.cl,shell=False)
+ retval = p.wait()
+ if self.opts.output_dir:
+ sto.close()
+ ste.close()
+ # get stderr, allowing for case where it's very large
+ tmp_stderr = open(self.elog, 'rb' )
+ stderr = ''
+ try:
+ while True:
+ stderr += tmp_stderr.read( buffsize )
+ if not stderr or len( stderr ) % buffsize != 0:
+ break
+ except OverflowError:
+ pass
+ tmp_stderr.close()
+ return retval
+
+
+def main():
+ u = """
+ This is a Galaxy wrapper. It expects to be called by a special purpose tool.xml as (eg):
+ rgToolFactory.py --script_path "$scriptPath" --tool_name "foo" --interpreter "Rscript"
+
+ The tool writes a script to a scriptPath using a configfile.
+ Some things in the script are templates.
+ The idea here is that depending on how this code is called, it uses the specified interpreter
+ to run a (hopefully correct) script/template. Optionally generates a clone of itself
+ which will run that same script/template as a toolshed repository tarball for uploading to a toolshed.
+ There's now an updated version which allows arbitrary parameters.
+ And so it goes.
+ """
+ op = optparse.OptionParser()
+ a = op.add_option
+ a('--script_path',default=None)
+ a('--tool_name',default=None)
+ a('--interpreter',default=None)
+ a('--output_dir',default=None)
+ a('--output_html',default=None)
+ a('--input_tab',default="None")
+ a('--output_tab',default="None")
+ a('--user_email',default='Unknown')
+ a('--bad_user',default=None)
+ a('--make_Tool',default=None)
+ a('--make_HTML',default=None)
+ a('--help_text',default=None)
+ a('--tool_desc',default=None)
+ a('--new_tool',default=None)
+ a('--tool_version',default=None)
+ opts, args = op.parse_args()
+ assert not opts.bad_user,'UNAUTHORISED: %s is NOT authorized to use this tool until Galaxy admin adds %s to admin_users in universe_wsgi.ini' % (opts.bad_user,opts.bad_user)
+ assert opts.tool_name,'## Tool Factory expects a tool name - eg --tool_name=DESeq'
+ assert opts.interpreter,'## Tool Factory wrapper expects an interpreter - eg --interpreter=Rscript'
+ assert os.path.isfile(opts.script_path),'## Tool Factory wrapper expects a script path - eg --script_path=foo.R'
+ if opts.output_dir:
+ try:
+ os.makedirs(opts.output_dir)
+ except:
+ pass
+ r = ScriptRunner(opts)
+ if opts.make_Tool:
+ retcode = r.makeTooltar()
+ else:
+ retcode = r.run()
+ os.unlink(r.sfile)
+ if retcode:
+ sys.exit(retcode) # indicate failure to job runner
+
+
+if __name__ == "__main__":
+ main()
+
+
diff -r 7d66bfa4fd56 -r 21b12c7c52e4 rg_nri.xml
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/rg_nri.xml Sat Oct 31 01:34:27 2015 -0400
@@ -0,0 +1,633 @@
+
+ and other model improvement measures
+
+ R_3_2_2
+ graphicsmagick
+ ghostscript
+ glmnet_lars_3_2
+
+
+ rgToolFactory.py --script_path "$runme" --interpreter "Rscript" --tool_name "rg_NRI"
+ --output_dir "$html_file.files_path" --output_html "$html_file" --make_HTML "yes"
+
+
+
+
+ 1) {
+ print.noquote("INPUT ERROR: x1 and x2 must be in [0,1]")
+ sink()
+ quit(save="no",status=3)
+ }
+ incidence=sum(y)/n
+ if (missing(threshvec)) {
+ threshvec=c(0, incidence,1)
+ print(paste('threshvec missing. using',paste(threshvec,collapse=',')))
+ }
+ a = (y == 1)
+ b = (y == 0)
+ na = sum(a)
+ nb = sum(b)
+ d = x2 - x1
+ ### NRI
+ n.thresh=length(threshvec)-1
+ risk.class.x1.ev=cut2(x1[a],threshvec)
+ risk.class.x2.ev=cut2(x2[a],threshvec)
+ thresh=c()
+ lt = length(threshvec)
+ for (i in 1:(lt-1)) {
+ thresh[i] = paste("[",toString(threshvec[i]),",",toString(threshvec[i+1]),"]")
+ }
+ levels(risk.class.x1.ev)=thresh
+ levels(risk.class.x2.ev)=thresh
+ cM.ev=confusionMatrix(risk.class.x2.ev,risk.class.x1.ev)
+ pup.ev=0
+ pdown.ev=0
+ for (i in 1:(n.thresh-1)) { pup.ev = pup.ev + sum(cM.ev\$table[(i+1):n.thresh,i])}
+ for (i in 2:n.thresh) { pdown.ev = pdown.ev + sum(cM.ev\$table[1:(i-1),i])}
+ pup.ev=pup.ev/na
+ pdown.ev=pdown.ev/na
+ risk.class.x1.ne=cut2(x1[b],threshvec)
+ risk.class.x2.ne=cut2(x2[b],threshvec)
+ levels(risk.class.x1.ne)=thresh
+ levels(risk.class.x2.ne)=thresh
+ cM.ne=confusionMatrix(risk.class.x2.ne,risk.class.x1.ne)
+ pup.ne=0
+ pdown.ne=0
+ for (i in 1:(n.thresh-1)){pup.ne=pup.ev+sum(cM.ne\$table[(i+1):n.thresh,i])}
+ for (i in 2:n.thresh){pdown.ne=pdown.ne+sum(cM.ne\$table[1:(i-1),i])}
+ pdown.ne=pdown.ne/nb
+ pup.ne=pup.ne/nb
+ nri = pup.ev - pdown.ev - (pup.ne - pdown.ne)
+ se.nri = sqrt((pup.ev + pdown.ev)/na + (pup.ne + pdown.ne)/nb)
+ z.nri = nri/se.nri
+ nri.ev = pup.ev - pdown.ev
+ se.nri.ev = sqrt((pup.ev + pdown.ev)/na)
+ z.nri.ev = nri.ev/se.nri.ev
+ nri.ne = pdown.ne - pup.ne
+ se.nri.ne = sqrt((pdown.ne + pup.ne)/nb)
+ z.nri.ne = nri.ne/se.nri.ne
+ ### Category Free NRI calculations
+ cfpup.ev = mean(d[a] > 0)
+ cfpup.ne = mean(d[b] > 0)
+ cfpdown.ev = mean(d[a] < 0)
+ cfpdown.ne = mean(d[b] < 0)
+ cfnri = cfpup.ev - cfpdown.ev - (cfpup.ne - cfpdown.ne)
+ se.cfnri = sqrt((cfpup.ev + cfpdown.ev)/na + (cfpup.ne + cfpdown.ne)/nb)
+ z.cfnri = cfnri/se.cfnri
+ cfnri.ev = cfpup.ev - cfpdown.ev
+ se.cfnri.ev = sqrt((cfpup.ev + cfpdown.ev)/na)
+ z.cfnri.ev = cfnri.ev/se.cfnri.ev
+ cfnri.ne = cfpdown.ne - cfpup.ne
+ se.cfnri.ne = sqrt((cfpdown.ne + cfpup.ne)/nb)
+ z.cfnri.ne = cfnri.ne/se.cfnri.ne
+ ### IDI calculations
+ improveSens = sum(d[a])/na
+ improveSpec = -sum(d[b])/nb
+ idi.ev = mean(improveSens)
+ idi.ne = mean(improveSpec)
+ idi = idi.ev - idi.ne
+ var.ev = var(d[a])/na
+ se.idi.ev = sqrt(var.ev)
+ z.idi.ev = idi.ev/se.idi.ev
+ var.ne = var(d[b])/nb
+ se.idi.ne = sqrt(var.ne)
+ z.idi.ne = idi.ne/se.idi.ne
+ se.idi = sqrt(var.ev + var.ne)
+ z.idi = idi/se.idi
+ ### AUC calculations
+ roc.x1 = roc(y, x1)
+ auc.x1 = auc(roc.x1)
+ ci.auc.x1 = ci.auc(roc.x1)
+ se.auc.x1 = (ci.auc.x1[3] - auc.x1)/qnorm(0.975)
+ roc.x2 = roc(y, x2)
+ auc.x2 = auc(roc.x2)
+ ci.auc.x2 = ci.auc(roc.x2)
+ se.auc.x2 = (ci.auc.x2[3] - auc.x2)/qnorm(0.975)
+ roc.test.x1.x2 = roc.test(roc.x1, roc.x2) ###Uses the default Delong method
+ sens.x1 = roc.x1\$sensitivities
+ spec.x1 = 1 - roc.x1\$specificities
+ n.x1 = length(sens.x1)
+ x1 = roc.x1\$thresholds
+ x1 = x1[c(-1,-n.x1)]
+ sens.x1 = sens.x1[c(-1,-n.x1)]
+ spec.x1 = spec.x1[c(-1,-n.x1)]
+ sens.x2 = roc.x2\$sensitivities
+ spec.x2 = 1 - roc.x2\$specificities
+ n.x2 = length(sens.x2)
+ x2 = roc.x2\$thresholds
+ x2 = x2[c(-1,-n.x2)]
+ sens.x2 = sens.x2[c(-1,-n.x2)]
+ spec.x2 = spec.x2[c(-1,-n.x2)]
+ ### Integrated sensitivity and 1-specificity calculations
+ is.x1 = trapz(x = x1, y = sens.x1) ### area under curves (relates to integrated sens, 1-spec)
+ is.x2 = trapz(x = x2, y = sens.x2)
+ ip.x1 = trapz(x = x1, y = spec.x1)
+ ip.x2 = trapz(x = x2, y = spec.x2)
+
+ ### Output
+ output = c(n, na, nb, pup.ev, pup.ne, pdown.ev, pdown.ne, nri, se.nri, z.nri,
+ nri.ev, se.nri.ev, z.nri.ev, nri.ne, se.nri.ne, z.nri.ne,
+ cfpup.ev, cfpup.ne, cfpdown.ev, cfpdown.ne, cfnri, se.cfnri, z.cfnri,
+ cfnri.ev, se.cfnri.ev, z.cfnri.ev, cfnri.ne, se.cfnri.ne, z.cfnri.ne,
+ improveSens, improveSpec, idi.ev, se.idi.ev, z.idi.ev, idi.ne,
+ se.idi.ne, z.idi.ne, idi, se.idi, z.idi, is.x1, NA, is.x2, NA,
+ ip.x1, NA, ip.x2, NA, auc.x1, se.auc.x1, auc.x2, se.auc.x2,
+ roc.test.x1.x2\$p.value,incidence)
+ names(output) = c("n", "na", "nb", "pup.ev", "pup.ne", "pdown.ev", "pdown.ne",
+ "nri", "se.nri", "z.nri", "nri.ev", "se.nri.ev", "z.nri.ev",
+ "nri.ne", "se.nri.ne", "z.nri.ne",
+ "cfpup.ev", "cfpup.ne", "cfpdown.ev", "cfpdown.ne",
+ "cfnri", "se.cfnri", "z.cfnri", "cfnri.ev", "se.cfnri.ev", "z.cfnri.ev",
+ "cfnri.ne", "se.cfnri.ne", "z.cfnri.ne", "improveSens", "improveSpec",
+ "idi.ev", "se.idi.ev", "z.idi.ev", "idi.ne", "se.idi.ne",
+ "z.idi.ne", "idi", "se.idi", "z.idi", "is.x1", "se.is.x1",
+ "is.x2", "se.is.x2", "ip.x1", "se.ip.x1", "ip.x2", "se.ip.x2",
+ "auc.x1", "se.auc.x1", "auc.x2", "se.auc.x2",
+ "roc.test.x1.x2.pvalue","incidence")
+ resdf = data.frame(N=n, Na=na, Nb=nb, pup.ev=pup.ev, pup.ne=pup.ne, pdown.ev=pdown.ev, pdown.ne=pdown.ne, NRI=nri, NRI.se=se.nri, NRI.z=z.nri,
+ NRI.ev=nri.ev, NRI.ev.se=se.nri.ev, NRI.ev.z=z.nri.ev, NRI.ne=nri.ne, NRI.ne.se=se.nri.ne, NRI.ne.z=z.nri.ne,
+ cfpup.ev=cfpup.ev, cfpup.ne=cfpup.ne, cfpdown.ev=cfpdown.ev, cfpdown.ne=cfpdown.ne, CFNRI=cfnri, CFNRI.se=se.cfnri, CFNRI.z=z.cfnri,
+ CFNRI.ev=cfnri.ev, CFNRI.ev.se=se.cfnri.ev, CFNRI.ev.z=z.cfnri.ev, CFNRI.ne=cfnri.ne, CFNRI.ne.se=se.cfnri.ne, CFNRI.ne.z=z.cfnri.ne,
+ improvSens=improveSens, improvSpec=improveSpec, IDI.ev=idi.ev, IDI.ev.se=se.idi.ev, IDI.ev.z=z.idi.ev, IDI.ne=idi.ne,
+ IDI.ne.se=se.idi.ne, IDI.ne.z=z.idi.ne, IDI=idi, IDI.se=se.idi, IDI.z=z.idi, isx1=is.x1, isx2=is.x2,
+ ipxi=ip.x1, ipx2=ip.x2, AUC.x1=auc.x1, AUC.x1.se=se.auc.x1, AUC.x2=auc.x2, AUC.x2.se=se.auc.x2,
+ roctestpval=roc.test.x1.x2\$p.value,incidence=incidence)
+ tr = t(resdf)
+ tresdf = data.frame(measure=colnames(resdf),value=tr[,1])
+ return(list(resdf=tresdf,output=output))
+}
+
+
+### More comprehensive summary statistics from a raplot
+### Choice of confidence intervals determined through asymptotics or bootstrapping (n.boot = ### of bootstrap resamples)
+### dp is number of decimal places for results table
+
+summary.raplot = function(x1, x2, y, threshvec, cis = "boot", conf.level = 0.95, n.boot = 2000, dp = 4, stat_ra=NA)
+{
+ results = stat_ra
+ if (cis == "boot") {
+ print.noquote("Bootstrap estimates for SE")
+ results.boot = matrix(NA, n.boot, length(names(results)))
+
+ colnames(results.boot) = names(results)
+
+ for (i in 1:n.boot) {
+ ###boot.index = sample(length(cc.status), replace = TRUE)
+ ###risk.model1.boot = risk.model1[boot.index]
+ ###risk.model2.boot = risk.model2[boot.index]
+ ###cc.status.boot = cc.status[boot.index]
+ boot.index = sample(length(y), replace = TRUE)
+ risk.model1.boot = x1[boot.index]
+ risk.model2.boot = x2[boot.index]
+ cc.status.boot = y[boot.index]
+ r = statistics.raplot(x1 = risk.model1.boot, x2 = risk.model2.boot, y = cc.status.boot)
+ results.boot[i, ] = r\$output
+ }
+
+ results.se.boot = apply(results.boot, 2, sd)
+ print(paste(results.se.boot,collapse=','))
+
+
+ results[grep("se", names(results))] = results.se.boot[grep("se", names(results)) - 1]
+
+ }
+
+
+
+ ### calculate cis and return
+
+ z = abs(qnorm((1 - conf.level)/2))
+
+ results.matrix = matrix(NA, 24, 2)
+
+ results.matrix[1, ] = c("Total (n)", results["n"])
+ results.matrix[2, ] = c("Events (n)", results["na"])
+ results.matrix[3, ] = c("Non-events (n)", results["nb"])
+ results.matrix[4, ] = c("Category free NRI and summary statistics","-------------------------")
+ results.matrix[5, ] = c("cfNRI events (%)",
+ paste(round(100*results["cfnri.ev"], dp-2), " (",
+ round(100*results["cfnri.ev"] - z * 100*results["se.cfnri.ev"], dp-2),
+ " to ", round(100*results["cfnri.ev"] +
+ z * 100*results["se.cfnri.ev"], dp-2), ")", sep = ""))
+ results.matrix[6, ] = c("cfNRI non-events (%)",
+ paste(round(100*results["cfnri.ne"], dp-2), " (",
+ round(100*results["cfnri.ne"] - z * 100*results["se.cfnri.ne"], dp)-2,
+ " to ", round(100*results["cfnri.ne"] + z * 100*results["se.cfnri.ne"],
+ dp-2), ")", sep = ""))
+ results.matrix[7, ] = c("cfNRI (%)",
+ paste(round(100*results["cfnri"], dp-2), " (",
+ round(100*results["cfnri"] - z * 100*results["se.cfnri"], dp-2),
+ " to ", round(100*results["cfnri"] + z * 100*results["se.cfnri"],
+ dp-2), ")", sep = ""))
+ results.matrix[8, ] = c("NRI and summary statistics","-------------------------")
+ results.matrix[9, ] = c("NRI events (%)",
+ paste(round(100*results["nri.ev"], dp-2), " (",
+ round(100*results["nri.ev"] - z * 100*results["se.nri.ev"], dp-2),
+ " to ", round(100*results["nri.ev"] +
+ z * 100*results["se.nri.ev"], dp-2), ")", sep = ""))
+ results.matrix[10, ] = c("NRI non-events (%)",
+ paste(round(100*results["nri.ne"], dp-2), " (",
+ round(100*results["nri.ne"] - z * 100*results["se.nri.ne"], dp-2),
+ " to ", round(100*results["nri.ne"] + z * 100*results["se.nri.ne"],
+ dp-2), ")", sep = ""))
+ results.matrix[11, ] = c("NRI (%)",
+ paste(round(100*results["nri"], dp-2), " (",
+ round(100*results["nri"] - z * 100*results["se.nri"], dp-2),
+ " to ", round(100*results["nri"] + z * 100*results["se.nri"],
+ dp-2), ")", sep = ""))
+ results.matrix[12, ] = c("IDI and summary statistics","-------------------------")
+ results.matrix[13, ] = c("IDI events",
+ paste(round(results["idi.ev"], dp), " (",
+ round(results["idi.ev"] - z * results["se.idi.ev"], dp),
+ " to ", round(results["idi.ev"] + z * results["se.idi.ev"],
+ dp), ")", sep = ""))
+ results.matrix[14, ] = c("IDI non-events",
+ paste(round(results["idi.ne"], dp), " (",
+ round(results["idi.ne"] - z * results["se.idi.ne"], dp),
+ " to ", round(results["idi.ne"] + z * results["se.idi.ne"],
+ dp), ")", sep = ""))
+ results.matrix[15, ] = c("IDI",
+ paste(round(results["idi"], dp), " (",
+ round(results["idi"] - z * results["se.idi"], dp),
+ " to ", round(results["idi"] + z * results["se.idi"],
+ dp), ")", sep = ""))
+ results.matrix[16, ] = c("IS (null model)",
+ paste(round(results["is.x1"], dp), " (",
+ round(results["is.x1"] - z * results["se.is.x1"], dp),
+ " to ", round(results["is.x1"] + z * results["se.is.x1"],
+ dp), ")", sep = ""))
+ results.matrix[17, ] = c("IS (alt model)",
+ paste(round(results["is.x2"], dp), " (",
+ round(results["is.x2"] - z * results["se.is.x2"], dp),
+ " to ", round(results["is.x2"] + z * results["se.is.x2"],
+ dp), ")", sep = ""))
+ results.matrix[18, ] = c("IP (null model)",
+ paste(round(results["ip.x1"], dp), " (",
+ round(results["ip.x1"] - z * results["se.ip.x1"], dp),
+ " to ", round(results["ip.x1"] + z * results["se.ip.x1"],
+ dp), ")", sep = ""))
+ results.matrix[19, ] = c("IP (alt model)",
+ paste(round(results["ip.x2"], dp), " (",
+ round(results["ip.x2"] - z * results["se.ip.x2"], dp),
+ " to ", round(results["ip.x2"] + z * results["se.ip.x2"],
+ dp), ")", sep = ""))
+ results.matrix[20, ] = c("AUC","-------------------------")
+ results.matrix[21, ] = c("AUC (null model)",
+ paste(round(results["auc.x1"], dp), " (",
+ round(results["auc.x1"] - z * results["se.auc.x1"], dp),
+ " to ", round(results["auc.x1"] + z * results["se.auc.x1"],
+ dp), ")", sep = ""))
+ results.matrix[22, ] = c("AUC (alt model)",
+ paste(round(results["auc.x2"], dp), " (",
+ round(results["auc.x2"] - z * results["se.auc.x2"], dp),
+ " to ", round(results["auc.x2"] + z * results["se.auc.x2"],
+ dp), ")", sep = ""))
+ results.matrix[23, ] = c("difference (P)", round(results["roc.test.x1.x2.pvalue"], dp))
+ results.matrix[24, ] = c("Incidence", round(results["incidence"], dp))
+
+ return(results.matrix)
+}
+
+
+
+]]>
+
+options(width=120)
+options(digits=5)
+logf = file("rgNRI.log", open = "a")
+sink(logf,type = c("output", "message"))
+Out_Dir = "$html_file.files_path"
+Input1 = "$input1"
+Input2 = "$input2"
+myTitle = "$title"
+outtab = "$nri_file"
+input1_obs = $input1_observed
+input1_pred = $input1_predicted
+input1_id = $input1_id
+input2_obs = $input2_observed
+input2_pred = $input2_predicted
+input2_id = $input2_id
+in1 = read.table(Input1,head=T,sep='\t')
+in2 = read.table(Input2,head=T,sep='\t')
+id1 = in1[,input1_id]
+id2 = in2[,input2_id]
+useme1 = in1[which(id1 %in% id2),]
+useme2 = in2[which(id2 %in% id1),]
+id1 = useme1[,input1_id]
+id2 = useme2[,input2_id]
+useme1 = useme1[order(id1),]
+useme2 = useme2[order(id2),]
+x1 = useme1[,input1_pred]
+x2 = useme2[,input2_pred]
+y1 = useme1[,input1_obs]
+y2 = useme2[,input2_obs]
+n.boot = $CImeth.nboot
+conf.level = 0.95
+cis = "$CImeth.cis"
+digits = 4
+nydiff = sum(y1 != y2)
+if (nydiff > 0) {
+ print.noquote(paste('Input error: observed status column has',nydiff,'differences - cannot reliably proceed'))
+ quit(save="no",status=1)
+ }
+y = y2
+outplot = 'rgNRI_EventRisk.pdf'
+res = raplot(x1=x1, x2=x2, y=y, outplot=outplot,title=myTitle)
+
+stats = statistics.raplot(x1=x1, x2=x2, y=y)
+res1 = stats\$resdf
+out1 = stats\$output
+print.noquote('Results:')
+print.noquote(res1,digits=4)
+res2 = summary.raplot(x1=x1, x2=x2, y=y, cis = cis, conf.level = conf.level, n.boot = n.boot, dp = digits, stat_ra=out1)
+print.noquote('Summary:')
+print.noquote(res2,digits=4)
+write.table(format(res1,digits=4),outtab,quote=F, col.names=F,sep="\t",row.names=F)
+print.noquote('SessionInfo for this R session:')
+sessionInfo()
+print.noquote('warnings for this R session:')
+warnings()
+sink()
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+**Before you start**
+
+This is a simple tool to calculate various measures of improvement in prediction between two models described in pickering_paper_
+It is based on an R script pickering_code_ written by Dr John W Pickering and Dr David Cairns from sunny Otago University which
+has been debugged and slightly adjusted to fit a Galaxy tool wrapper.
+
+
+**What it does**
+
+Copied from the documentation in pickering_code_ ::
+
+
+ Functions to create risk assessment plots and associated summary statistics
+
+
+ (c) 2012 Dr John W Pickering, john.pickering@otago.ac.nz, and Dr David Cairns
+ Last modified August 2014
+
+ Redistribution and use in source and binary forms, with or without
+ modification, are permitted provided that the following conditions are met:
+ * Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+ * Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer in
+ the documentation and/or other materials provided with the distribution
+
+ FUNCTIONS
+ raplot
+ Produces a Risk Assessment Plot and outputs the coordinates of the four curves
+ Based on: Pickering, J. W. and Endre, Z. H. (2012). New Metrics for Assessing Diagnostic Potential of
+ Candidate Biomarkers. Clinical Journal of the American Society of Nephrology, 7, 1355–1364. doi:10.2215/CJN.09590911
+
+ statistics.raplot
+ Produces the NRIs, IDIs, IS, IP, AUCs.
+ Based on: Pencina, M. J., D'Agostino, R. B. and Steyerberg, E. W. (2011). Extensions of net reclassification improvement calculations to
+ measure usefulness of new biomarkers. Statistics in Medicine, 30(1), 11–21. doi:10.1002/sim.4085
+ Pencina, M. J., D'Agostino, R. B. and Vasan, R. S. (2008). Evaluating the added predictive ability of a new marker: From area under the
+ ROC curve to reclassification and beyond.
+ Statistics in Medicine, 27(2), 157–172. doi:10.1002/sim.2929
+ DeLong, E., DeLong, D. and Clarke-Pearson, D. (1988). Comparing the areas under 2 or more correlated receiver operating characteristic curves - a nonparametric approach.
+ Biometrics, 44(3), 837–845.
+
+ summary.raplot
+ Produces the NRIs, IDIs, IS, IP, AUCs with confidence intervals using a bootstrap or asymptotic procedure. (I prefer bootstrap which is chosed by cis=c("boot"))
+
+
+ Required arguments for all functions:
+ x1 is calculated risk (eg from a glm) for the null model, i.e. predict(,type="response") on a glm object
+ x2 is calculated risk (eg from a glm) for the alternative model
+ y is the case-control indicator (0 for controls, 1 for cases)
+ Optional argument
+ t are the boundaries of the risks for each group (ie 0, 1 and the thresholds beteween. eg c(0,0,3,0,7,1)). If missing, defaults to c(0, the incidence, 1)
+
+
+**Input**
+
+The observed and predicted outcomes from two models to be compared.
+
+**Output**
+
+Lots'o'measures (TM) see pickering_paper_ for details
+
+**Attributions**
+
+pickering_paper_ is the paper the caclulations performed by this tool is based on
+
+pickering_code_ is the R function from John Pickering exposed by this Galaxy tool with minor modifications and hacks by Ross Lazarus.
+
+Galaxy_ (that's what you are using right now!) for gluing everything together
+
+Otherwise, all code and documentation comprising this tool was written by Ross Lazarus and is
+licensed to you under the LGPL_ like other rgenetics artefacts
+
+.. _LGPL: http://www.gnu.org/copyleft/lesser.html
+.. _pickering_code: http://www.researchgate.net/publication/264672640_R_function_for_Risk_Assessment_Plot__reclassification_metrics_NRI_IDI_cfNRI
+.. _pickering_paper: http://cjasn.asnjournals.org/content/early/2012/05/24/CJN.09590911.full
+.. _Galaxy: http://getgalaxy.org
+
+
+
+
+
+ doi: 10.2215/CJN.09590911
+
+
+
+
+
diff -r 7d66bfa4fd56 -r 21b12c7c52e4 rglasso_cox.xml
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/rglasso_cox.xml Sat Oct 31 01:34:27 2015 -0400
@@ -0,0 +1,910 @@
+
+ and cox regression using elastic net
+
+ R_3_2_2
+ graphicsmagick
+ ghostscript
+ glmnet_lars_3_2
+
+
+ rgToolFactory.py --script_path "$runme" --interpreter "Rscript" --tool_name "rglasso"
+ --output_dir "$html_file.files_path" --output_html "$html_file" --make_HTML "yes"
+
+
+
+ 0 ) {
+ predcols = inmodel
+ xmat = as.matrix(x[,predcols])
+ colnames(xmat) = preds
+ bestpred = getpredp_cox(x=xmat,time=y[,'time'],status=y[,'status'],id=cox_id, predict_at=predict_at)
+ pred = data.frame(responsep=bestpred, best_lambda=best_lambda,lamchoice=optLambda,alpha=alpha)
+ write.table(pred,out_pred_file,quote=FALSE, sep="\t",row.names=F)
+ } else { print.noquote('WARNING: No coefficients in selected model to predict with - no predictions made') }
+ }
+ if (debugOn) {
+ print.noquote(paste('best_lambda=',best_lambda,'saving cox respreds=',paste(names(coefs),collapse=','),'as predictors of survival. Coefs=',paste(coefs,collapse=',')))
+ }
+ p = try({res = data.frame(regulator=names(coefs),partial_likelihood=coefs,forced_in=pen,glmnet_model='cox',best_lambda=best_lambda,
+ lambdaChoice=optLambda,alpha=alpha)},T)
+ if (class(p) == "try-error") {
+ message(paste('@@@ unable to return a dataframe',p))
+ sink()
+ return(NA)
+ }
+ print.noquote('@@@ Results preview:')
+ print.noquote(res,digits=5)
+ sink()
+ return(res)
+
+}
+
+
+do_lasso = function(x=NA,y=NA,do_standard=T,debugOn=T,defaultFam="gaussian",optLambda='minLambda',descr='description', indx=1,target='target',sane=F,
+ alpha=0.9,nfold=10,penalty=c(),out_pred=F,out_pred_file='outpred',yvarname='yvar',id=c(),mdsplots=F)
+{
+ logf = file(paste(target,"rglasso.log",sep='_'), open = "a")
+ sink(logf,type = c("output", "message"))
+ res = NA
+ phe_is_bin = (length(unique(y)) == 2)
+ forcedin = paste(colnames(x)[which(penalty == 0)],collapse=',')
+ fam = "gaussian"
+ if (defaultFam %in% c("poisson","binomial","gaussian","multinomial")) fam=defaultFam
+ if (phe_is_bin == T) {
+ fam = "binomial"
+ }
+ print.noquote(paste('target=',target,'is binary=',phe_is_bin,'dim(x)=',paste(dim(x),collapse=','),'length(y)=',length(y),'force=',forcedin,'fam=',fam))
+ standardize = do_standard
+ p = try({larsres = glmnet(x,y,family=fam,standardize=standardize,maxit=10000,alpha=alpha,penalty.factor=penalty) },T)
+ if (class(p) == "try-error")
+ {
+ print(paste('ERROR: unable to run glmnet for target',target,'error=',p))
+ sink()
+ return(NA)
+ }
+
+ mt = paste('Glmnet fraction deviance for',target)
+ outpdf = paste(target,'glmnetPath.pdf',sep='_')
+ pdf(outpdf)
+ plot(larsres,main=mt,label=T)
+ grid()
+ dev.off()
+
+ outpdf = paste(target,'glmnetDeviance.pdf',sep='_')
+
+ mt2 = paste('Glmnet lambda for',target)
+
+ pdf(outpdf)
+ plot(larsres,xvar="lambda",main=mt2,label=T)
+ grid()
+ dev.off()
+
+ larscv = NA
+ if (fam=="binomial") {
+ tmain = paste(target,'AUC')
+ outpdf = paste(target,'glmnetCV_AUC.pdf',sep='_')
+ p = try({larscv = ross.cv.glmnet(x=x,y=y,family=fam,type.measure='auc')},T)
+ } else {
+ tmain = paste(target,'CV MSE')
+ outpdf = paste(target,'glmnetCV_MSE.pdf',sep='_')
+ p = try({larscv = ross.cv.glmnet(x=x,y=y,family=fam,type.measure='mse')},T)
+ }
+ if (class(p) == "try-error")
+ {
+ print(paste('ERROR: unable to run cross validation for target',target,'error=',p))
+ sink()
+ return(NA)
+ }
+
+ pdf(outpdf)
+ plot(larscv,main=tmain)
+ grid()
+ dev.off()
+
+ lse = larscv\$cvhits.1se
+ lmin = larscv\$cvhits.min
+ tot = lse + lmin
+ allhits = data.frame(pred=colnames(x),hits_lambda_1se = lse,hits_lambda_min = lmin)
+ nzhits = allhits[which(tot != 0),]
+ message('Total hit count for each predictor over all CV models (excluding zero count predictors)')
+ print.noquote(nzhits)
+ out_nz_file = paste(target,'cross_validation_model_counts.xls',sep='_')
+ write.table(nzhits,out_nz_file,quote=FALSE, sep="\t",row.names=F)
+
+ ipenalty = c(0,penalty)
+ if (optLambda == 'lambda.min') {
+ best_lambda = larscv\$lambda.min
+ bestpred = as.matrix(coef(larscv, s = "lambda.min"))
+ } else {
+ best_lambda = larscv\$lambda.1se
+ bestpred = as.matrix(coef(larscv, s = "lambda.1se"))
+ }
+ inmodel = which(bestpred != 0)
+ coefs = bestpred[inmodel,1]
+ preds = rownames(bestpred)[inmodel]
+ iforced = ipenalty[inmodel]
+ forced = ! as.logical(iforced)
+ names(coefs) = preds
+ ncoef = length(coefs) - 1
+ if (out_pred==T && fam=="binomial")
+ {
+ print.noquote(paste('Predicting',target,'probabilities from binomial glmnet at alpha',alpha,'and lambda',best_lambda))
+ bestpred = predict.glmnet(larsres,s=best_lambda,newx=x,type="response")
+ bestpred = exp(bestpred)/(1+exp(bestpred))
+ pred = data.frame(id=id,y=y,predp=as.vector(bestpred), best_lambda=best_lambda)
+ write.table(pred,out_pred_file,quote=FALSE, sep="\t",row.names=F)
+ }
+ if (debugOn) {cat(indx,'best_lambda=',best_lambda,'saving',fam,'respreds=',names(coefs),'as predictors of',target,'coefs=',coefs,'\n')}
+ res = try(data.frame(i=indx,pred=target,regulator=names(coefs),coef=coefs,forced_in=forced,glmnet_model=fam,ncoef=ncoef,
+ best_lambda=best_lambda,lambdaChoice=optLambda,alpha=alpha),T)
+ if (class(res) == "try-error") {
+ sink()
+ return(NA) }
+ print.noquote(res)
+ sink()
+ return(res)
+}
+
+
+dolasso_generic = function(predvars=NA,depvars=NA,debugOn=T,maxsteps=100, alpha=0.9,nfold=10,xcolnames=c(),ycolnames=c(),optLambda='minLambda', out_pred_file=NA,
+ descr="describe me",do_standard=F,defaultFam="gaussian",penalty=c(),out_pred=F,cox_id=c(),mdsplots=F,xfilt=0.95)
+{
+ logf = file("rglasso.log", open = "a")
+ sink(logf,type = c("output", "message"))
+ xdat = predvars
+ xm = data.matrix(xdat)
+ res = NULL
+ id = cox_id
+ depnames = ycolnames
+ ndep = length(depnames)
+ if (mdsplots==T) {
+ outpdfname = 'rglasso_x_in_sample_space_MDS.pdf'
+ p = try({ mdsPlot(xm,'measurements in sample space',groups=NA,outpdfname=outpdfname,transpose=T) },T)
+ if (class(p) == "try-error")
+ {
+ print.noquote(paste('Unable to produce predictors in sample space mds plot. Error:',p))
+ }
+ outpdfname = 'rglasso_samples_in_x_space_MDS.pdf'
+ p = try({mdsPlot(xm,'samples in measurement space',groups=NA,outpdfname=outpdfname,transpose=F) },T)
+ if (class(p) == "try-error")
+ {
+ print.noquote(paste('Unable to produce samples in measurement space mds plot. Error:',p))
+ }
+ }
+ ndat = nrow(xm)
+ cfracs = colSums(! is.na(xm))/ndat
+ keepme = (cfracs >= xfilt)
+ print.noquote(paste('Removing', sum(! keepme), 'xvars with more than',xfilt,'fraction missing'))
+ vars = apply(xm,2,var,na.rm=T)
+ xm = xm[,keepme]
+ for (i in c(1:max(1,ndep))) {
+ target = depnames[i]
+ if (length(target) < 1) { target='y' }
+ if (i %% 100 == 0) { cat(i,target,'\n') }
+ if (ndep <= 1) {
+ y=depvars
+ } else {
+ y = depvars[,i]
+ }
+ if (fam == "binomial") {y = as.factor(y)}
+ x = xm
+ id = cox_id
+ if (fam != "cox") {
+ ok = complete.cases(x,y)
+ if (sum(! ok) > 0) {
+ message(paste('@@@ Removing',sum(! ok),'cases with missing y of',length(y),'@@@'))
+ y = y[(ok)]
+ x = x[(ok),]
+ id = id[(ok)]
+ }
+ }
+ ok = complete.cases(y)
+ if (sum(ok) == 0 ) {
+ print(paste("No complete cases found for",target,"in input x dim =",paste(dim(xm),collapse=','),"length y=",length(y)))
+ } else {
+ if (i == 1) { outpred = out_pred_file
+ } else {
+ outpred = paste(target,'predicted_output.xls')
+ }
+ regres = do_lasso(x=x,y=y,do_standard=do_standard,debugOn=debugOn,defaultFam=defaultFam,optLambda=optLambda,out_pred_file=outpred,
+ descr=descr,indx=i,target=target,alpha=alpha,nfold=nfold,penalty=penalty,out_pred=out_pred,yvarname=target,id=id,mdsplots=mdsplots)
+ if (! is.na(regres)) { res = rbind(res,regres) }
+ }
+ }
+ print.noquote('@@@ Results preview:')
+ print.noquote(res,digits=5)
+ sink()
+ return(res)
+}
+
+
+corPlot=function(xdat=c(),main='main title',is_raw=T)
+{
+ library(pheatmap)
+ library(gplots)
+ if (is_raw) {
+ cxdat = cor(xdat,method="spearman",use="pairwise.complete.obs")
+ } else {
+ cxdat=xdat
+ }
+ xro = nrow(cxdat)
+ if (xro > 1000) stop("Too many rows for heatmap, who can read?!")
+ fontsize_col = 5.0
+ pheatmap(cxdat, main=main, show_colnames = F, width=30, height=30,
+ fontsize_row=fontsize_col, border_color=NA)
+}
+
+
+runTest = function(n=10)
+{
+ set.seed (NULL)
+ Y = data.frame(y1=runif (n),y2=runif(n))
+ Xv <- runif(n*n)
+ X <- matrix(Xv, nrow = n, ncol = n)
+
+ mydf <- data.frame(Y, X)
+
+ regres_out = dolasso_generic(predvars=X,depvars=Y,debugOn=T,p.cutoff = 0.05,maxsteps=10000,nfold=10,
+ descr='randomdata',do_standard=do_standard,defaultFam="gaussian",alpha=0.05)
+ return(regres_out)
+}
+]]>
+options(width=512)
+options(digits=5)
+alpha = $alpha
+nfold = $nfold
+optLambda = "$optLambda"
+Out_Dir = "$html_file.files_path"
+Input = "$input1"
+indat = read.table(Input,head=T,sep='\t')
+datcols = colnames(indat)
+myTitle = "$title"
+outtab = "$model_file"
+do_standard = as.logical("$do_standard")
+mdsplots = as.logical("$mdsplots")
+fam = "$model.fam"
+xvar_cols_in = "$xvar_cols"
+force_xvar_cols_in = "$force_xvar_cols"
+xvar_cols = as.numeric(strsplit(xvar_cols_in,",")[[1]])
+force_xvar_cols = c()
+penalties = rep(1,length(datcols))
+forced_in = NA
+
+logxform = "$logxform_cols"
+if (logxform != "None") {
+ logxform_cols = as.numeric(strsplit(logxform,",")[[1]])
+ if (length(logxform_cols) > 0) {
+ small = 1e-10
+ sset = indat[,logxform_cols]
+ zeros = which(sset==0,arr.ind=T)
+ nz = nrow(zeros)
+ if (nz > 0) {
+ message(paste('Log transforming encountered',nz,'zeros - added 1e-10'))
+ sset[zeros] = sset[zeros] + small
+ lset = log(sset)
+ indat[,logxform_cols] = lset
+ }
+ }
+}
+if (force_xvar_cols_in != "None")
+{
+ force_xvar_cols = as.numeric(strsplit(force_xvar_cols_in,",")[[1]])
+ allx = c(xvar_cols,force_xvar_cols)
+ xvar_cols = unique(allx)
+ xvar_cols = xvar_cols[order(xvar_cols)]
+ penalties[force_xvar_cols] = 0
+}
+penalty = penalties[xvar_cols]
+forcedin = paste(datcols[which(penalties == 0)],collapse=',')
+cox_id_col = NA
+cox_id = NA
+
+message(paste('@@@ Using alpha =',alpha,'for all models'))
+x = indat[,xvar_cols]
+nx = nrow(x)
+cx = ncol(x)
+message(paste('@@@@ Input has',nx,'samples and',cx,'predictors'))
+if (cx > nx) {
+message('@@@ WARNING: Models will have more variables than cases so glmnet will likely return one of many possible solutions! Please DO NOT expect reliable results - glmnet is clever but not magical @@@')
+}
+
+xcolnames = datcols[xvar_cols]
+
+if (file.exists(Out_Dir) == F) dir.create(Out_Dir)
+out_full = F
+out_full_file = NA
+out_pred_file = ""
+out_pred = as.logical("$model.output_pred")
+
+#if $model.fam == "binomial" or $model.fam == "cox":
+ cox_id_col = $model.cox_id
+ cox_id = indat[,cox_id_col]
+ if (out_pred == T) {
+ out_pred_file="$output_pred_file"
+ rownames(x) = cox_id
+ }
+#end if
+#if $model.fam == "cox":
+ cox_time = $model.cox_time
+ cox_status = $model.cox_status
+ out_full = as.logical("$model.output_full")
+ if (out_full == T) { out_full_file="$output_full_file" }
+ yvar_cols = c(cox_time,cox_status)
+ ycolnames = c('time','status')
+ istat = as.double(indat[,cox_status])
+ itime = as.double(indat[,cox_time])
+ predict_at = quantile(itime)
+ if ("$model.predict_at" > "")
+ {
+ pa = "$model.predict_at"
+ predict_at = as.numeric(strsplit(pa,",")[[1]])
+ }
+ y = data.frame(time = itime, status = istat)
+ ustat = unique(istat)
+ if ((length(ustat) != 2) | (! 1 %in% ustat ) | (! 0 %in% ustat))
+ {
+ print.noquote(paste('INPUT ERROR: status must have 0 (censored) or 1 (event) but found',paste(ustat,collapse=',') ))
+ quit(save='no',status=1)
+ }
+ y = as.matrix(y)
+ x = as.matrix(x)
+ print.noquote(paste('@@@ Cox model will predict yvar=',datcols[cox_status],'using cols=',paste(xcolnames,collapse=','),'n preds=',length(xcolnames),
+ 'forced in=',forcedin))
+ regres_out = dolasso_cox(x=x,y=y,debugOn=F,maxsteps=10000,nfold=nfold,xcolnames=xcolnames,ycolnames=ycolnames,optLambda=optLambda,out_full=out_full,out_full_file=out_full_file,
+ out_pred=out_pred,out_pred_file=out_pred_file,cox_id=cox_id,descr=myTitle,do_standard=do_standard,alpha=alpha,penalty=penalty,predict_at=predict_at,mdsplots=mdsplots)
+#else:
+ yvar_cols = "$model.yvar_cols"
+ yvar_cols = as.numeric(strsplit(yvar_cols,",")[[1]])
+ ycolnames = datcols[yvar_cols]
+ print.noquote(paste('@@@',fam,'model will predict yvar=',paste(ycolnames,collapse=','),'using cols=',paste(xcolnames,collapse=','),'n preds=',length(xcolnames),
+ 'forced in=',forcedin))
+ y = data.matrix(indat[,yvar_cols])
+ print.noquote(paste('Model will use',fam,'link function to predict yvar=',paste(ycolnames,collapse=','),'n preds=',length(xcolnames),'forced in=',forcedin))
+ regres_out = dolasso_generic(predvars=x,depvars=y,debugOn=F, maxsteps=10000,nfold=nfold,xcolnames=xcolnames,ycolnames=ycolnames,optLambda=optLambda,out_pred_file=out_pred_file,
+ descr=myTitle,do_standard=do_standard,defaultFam=fam,alpha=alpha,penalty=penalty,out_pred=out_pred,cox_id=cox_id,mdsplots=mdsplots)
+#end if
+
+write.table(format(regres_out,digits=5),outtab,quote=FALSE, sep="\t",row.names=F)
+print.noquote('@@@ SessionInfo for this R session:')
+sessionInfo()
+warnings()
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ l
+
+
+
+ l
+
+
+
+
+
+
+
+
+
+ model['output_full'] == 'T'
+
+
+ model['output_pred'] == 'T'
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+**Before you start**
+
+Please read the glmnet documentation @ glmnet_
+
+This Galaxy wrapper merely exposes that code and the glmnet_ documentation is essential reading
+before getting useful results here.
+
+**What it does**
+
+From documentation at glmnet_ ::
+
+ Glmnet is a package that fits a generalized linear model via penalized maximum likelihood.
+ The regularization path is computed for the lasso or elasticnet penalty at a grid of values for the regularization parameter lambda.
+ The algorithm is extremely fast, and can exploit sparsity in the input matrix x.
+ It fits linear, logistic and multinomial, poisson, and Cox regression models.
+ A variety of predictions can be made from the fitted models.
+
+Internal cross validation is used to optimise the choice of lambda based on CV AUC for logistic (binomial outcome) models, or CV mse for gaussian.
+
+**Warning about the tyrany of dimensionality**
+
+Yes, this package will select 'optimal' models even when you (optimistically) supply more predictors than you have cases.
+The model returned is unlikely to represent the only informative regularisation path through your data - if you run repeatedly with
+exactly the same settings, you will probably see many different models being selected.
+This is not a software bug - the real problem is that you just don't have enough information in your data.
+
+Sufficiently big jobs will take a while (eg each lasso regression with 20k features on 1k samples takes about 2-3 minutes on our aged cluster)
+
+**Input**
+
+Assuming you have more measurements than samples, you supply data as a tabular text file where each row is a sample and columns
+are variables. You specify which columns are dependent (predictors) and which are observations for each sample. Each of multiple
+dependent variable columns will be run and reported independently. Predictors can be forced in to the model.
+
+**Output**
+
+For each selected dependent regression variable, a brief report of the model coefficients predicted at the
+'optimal' nfold CV value of lambda.
+
+**Predicted event probabilities for Cox and Logistic models**
+
+If you want to compare (eg) two competing clinical predictions, there's a companion generic NRI tool
+for predicted event probabilities. Estimates dozens of measures of improvement in prediction. Currently only works for identical id subjects
+but can probably be extended to independent sample predictions.
+
+Given a model, we can generate a predicted p (for status 1) in binomial or cox frameworks so models can be evaluated in terms of NRI.
+Of course, estimates are likely substantially inflated over 'real world' performance by being estimated from the same sample - but you probably
+already knew that since you were smart enough to reach this far down into the on screen help. The author salutes you, intrepid reader!
+
+It may seem an odd thing to do, but we can predict p for an event for each subject from our original data, given a parsimonious model. Doing
+this for two separate models (eg, forcing in an additional known explanatory measurement to the new model) allows comparison of the two models
+predicted status for each subject, or the same model in independent populations to see how badly it does
+
+**Attributions**
+
+glmnet_ is the R package exposed by this Galaxy tool.
+
+Galaxy_ (that's what you are using right now!) for gluing everything together
+
+Otherwise, all code and documentation comprising this tool was written by Ross Lazarus and is
+licensed to you under the LGPL_ like other rgenetics artefacts
+
+.. _LGPL: http://www.gnu.org/copyleft/lesser.html
+.. _glmnet: http://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html
+.. _Galaxy: http://getgalaxy.org
+
+
+
+
+@Article{Friedman2010, title = {Regularization Paths for Generalized Linear Models via Coordinate Descent},
+ author = {Jerome Friedman and Trevor Hastie and Robert Tibshirani},
+ journal = {Journal of Statistical Software},
+ year = {2010},
+ volume = {33},
+ number = {1},
+ pages = {1--22},
+ url = {http://www.jstatsoft.org/v33/i01/}
+ }
+
+
+10.1093/bioinformatics/bts573
+
+
+
diff -r 7d66bfa4fd56 -r 21b12c7c52e4 test-data/cox_coxlassotest_glmnet_cvdeviance.pdf
Binary file test-data/cox_coxlassotest_glmnet_cvdeviance.pdf has changed
diff -r 7d66bfa4fd56 -r 21b12c7c52e4 test-data/cox_coxlassotest_glmnetdev.pdf
Binary file test-data/cox_coxlassotest_glmnetdev.pdf has changed
diff -r 7d66bfa4fd56 -r 21b12c7c52e4 test-data/cox_test.xls
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/cox_test.xls Sat Oct 31 01:34:27 2015 -0400
@@ -0,0 +1,502 @@
+time status rhubarb vegemite apple
+575.966708 0 452.405468 30.339584 32286.089057
+1539.245319 1 329.689929 30.603839 15735.863202
+2072.798422 0 534.379263 32.474983 22639.136685
+1638.450154 0 362.522161 30.925996 31108.364370
+1771.625630 0 417.167751 31.444652 31317.491931
+4413.484706 1 459.947526 30.069033 11232.037882
+1519.376431 1 576.868414 30.102709 16539.948172
+4805.620290 0 554.756164 30.685550 9979.125726
+657.497124 0 439.481951 32.305768 22301.057122
+181.498205 1 581.498091 32.549519 10499.737992
+1630.552909 0 580.670671 30.825057 22321.821037
+2285.043061 1 588.034529 32.799742 8067.710675
+3697.405336 0 400.195630 31.826099 20975.730587
+3821.226480 0 561.466202 32.333230 17887.547850
+4810.837196 0 447.004938 32.547717 25383.925946
+3195.877681 1 482.169788 31.420740 27659.409774
+4556.750833 1 415.144072 31.960725 18401.138992
+3797.192061 1 484.492359 30.789882 7431.214452
+3734.484446 1 587.339066 32.052766 29443.923798
+567.665985 1 377.599777 31.203617 30148.076752
+2551.367122 0 447.623689 31.795708 15922.972458
+3.724234 1 534.080187 31.175131 8259.786167
+3528.906130 0 471.293653 31.225006 18307.543198
+3649.459162 1 443.804109 32.390564 9720.974163
+1490.640766 1 456.837379 32.099378 15292.782705
+4023.629633 1 443.869307 30.022985 23145.513801
+1389.288557 0 497.850711 30.361120 28757.227200
+535.427950 1 437.634304 32.307966 18289.733136
+3115.819953 0 587.593182 31.317353 21533.733397
+1561.895566 1 472.250365 31.550546 8724.067134
+2579.785003 0 331.955830 32.992176 12105.988070
+1592.819387 0 361.242026 31.011367 12823.339723
+4473.715814 0 334.539193 30.448049 9621.600487
+4328.198133 0 537.782739 32.683122 19090.292579
+3761.244687 1 392.916387 32.665889 22119.375302
+4638.730675 0 514.168146 30.548414 20651.266975
+1531.289211 0 451.404672 32.437627 21748.461852
+3844.140593 0 427.404840 32.894070 25824.198908
+4542.766339 0 366.992470 30.814372 6090.765585
+3296.236419 0 575.910564 32.509193 4050.864668
+710.915095 0 395.062092 32.940952 9086.464980
+3211.769396 0 428.608771 30.918481 29996.819084
+4448.442217 0 300.153203 31.160941 8735.458385
+4215.329607 0 504.064526 32.040819 29036.000397
+3249.725224 0 489.538009 32.863547 27520.919952
+3388.033079 0 485.938597 30.920035 20395.348731
+1465.239279 0 435.641977 30.193093 20469.748911
+2022.519649 0 362.011214 32.756787 18367.094394
+1474.662644 0 448.727985 31.194883 18456.157990
+2680.178585 0 553.521680 32.339108 17386.094437
+254.953131 1 467.930054 32.705913 28823.817908
+2080.749522 0 387.588403 30.520279 3389.418182
+4609.398123 1 471.568096 31.101077 27574.777543
+2506.254687 1 466.233702 32.928613 24693.847389
+1418.686322 0 519.402535 32.685931 9055.300824
+1716.501599 0 305.671024 31.453871 25382.505028
+2080.668376 0 589.514396 31.143248 24546.670311
+4051.441999 1 431.028785 30.854820 32291.736106
+221.739141 1 435.123612 30.465214 19077.175842
+3945.564918 0 559.566964 32.871624 18320.616591
+3543.190633 1 548.960006 30.269711 6259.113819
+4141.263169 1 418.013926 30.159280 12806.181141
+4613.643735 1 467.593153 30.584036 28998.377681
+2526.696757 1 442.602416 31.075723 24392.589804
+661.056610 0 306.186743 30.071900 4931.492293
+420.059232 0 557.043309 30.889234 30410.871000
+3777.621646 1 564.397202 30.810439 25359.485010
+4210.010862 0 541.943504 31.109403 18240.759045
+2251.752128 1 442.108369 31.308836 14136.141189
+4037.391365 1 467.892619 32.089289 14540.468274
+2478.012769 0 305.600832 30.653583 7117.580065
+2435.693015 0 551.918277 30.355254 29969.743529
+3816.690227 0 314.516061 31.511033 9880.067119
+531.400047 0 565.141705 30.722177 11252.478408
+4367.861436 0 421.500739 31.696812 16879.486148
+634.857316 1 476.587731 30.741441 31351.695315
+2757.953233 1 484.171370 31.510266 7767.961089
+676.813083 0 412.948756 30.174301 11123.425270
+3672.369695 0 374.771725 30.531392 32349.940422
+2572.842389 0 554.137199 32.777709 12727.962812
+4960.699040 0 560.164589 30.206776 3453.492367
+4316.475328 1 484.168551 32.069106 31984.273406
+4015.881070 1 400.274241 30.658886 11783.352436
+3052.451190 1 398.216589 30.358316 27928.803920
+1398.823206 0 446.521178 30.185715 31193.500034
+3922.909840 1 575.609185 32.724842 10233.363838
+182.803149 0 370.875903 30.856525 23954.116536
+2632.176915 1 426.948756 30.641145 30860.455084
+1060.390433 0 300.347626 31.143117 23615.898639
+3178.728283 1 514.261332 31.356224 26944.147061
+763.218442 1 463.367562 31.274306 7509.665226
+953.708474 0 384.301310 31.065789 15102.331493
+4434.527995 1 431.967331 32.369140 7177.785265
+623.914857 0 495.381873 32.915097 22739.985454
+3157.386533 1 515.293513 32.210421 31073.191185
+521.293142 0 593.449301 30.668076 29606.157493
+2845.587324 1 367.180408 32.424788 7052.350189
+4832.119270 1 465.679571 30.019029 5249.534096
+2856.445302 1 591.090965 31.133381 32498.377584
+2572.527257 1 409.021421 31.885768 13743.016360
+1413.194580 1 330.844512 32.740409 5441.606816
+2918.271731 0 304.169439 31.256918 22384.990257
+1674.024061 1 589.100300 32.763039 13305.351548
+2808.138084 0 491.357574 31.150814 17202.472869
+1424.663181 0 359.990029 31.343898 28564.024266
+2718.446177 1 580.778432 31.842978 29231.311831
+3735.730841 0 573.549637 30.079704 31952.172593
+4641.585855 0 365.172295 32.894223 17289.384320
+4071.637795 0 505.251583 32.472570 9633.518042
+471.788564 1 322.257982 32.546866 30076.426907
+1266.346159 1 335.749836 31.646670 32477.736291
+2695.370600 1 509.475837 31.947731 22351.940960
+3631.475811 0 584.575140 30.366829 13700.704821
+488.174184 0 522.383785 30.637536 29035.941989
+2645.805895 1 390.848623 31.639290 3004.701971
+4037.703358 1 425.444485 32.775358 28720.752294
+2196.610999 1 547.619331 32.640598 29108.134173
+1081.264261 1 350.044271 32.021783 19844.288986
+1074.305164 1 497.881588 32.048148 16325.089356
+1643.851420 0 372.156064 30.815571 12924.845649
+4919.287039 0 470.345027 31.043515 6953.198088
+3961.891671 1 353.750322 32.386195 10074.185995
+1643.427847 0 521.758887 31.784944 10303.308946
+2714.177818 0 363.519932 31.954710 4089.970927
+1958.737870 0 372.752456 32.778447 25735.088968
+72.430537 1 374.015311 30.406160 3519.910494
+666.824775 0 443.814245 31.382885 32222.429084
+179.203304 0 339.335011 30.911475 29294.393947
+3193.179459 1 329.239965 32.270432 7539.157039
+1405.955780 1 575.062041 31.769778 26056.208652
+4378.622335 1 583.779711 30.280724 12958.062814
+4933.985803 0 474.998284 30.835755 8298.116085
+479.970673 0 311.469610 31.110600 16498.771798
+194.917302 1 533.735116 32.682195 17363.293404
+3126.230367 1 568.043474 32.209623 22525.398240
+3605.775657 0 519.068951 30.750083 21588.128190
+68.236883 1 316.956971 32.735304 5303.699685
+2956.632516 1 359.835070 31.735527 26701.676218
+1273.924283 1 346.804469 31.185254 29977.854210
+1530.707021 1 568.442488 31.510433 20716.251168
+3847.696838 1 487.304670 32.242322 5276.023316
+3121.510941 0 518.612631 30.448810 10702.449428
+2693.401144 1 533.008711 32.964097 31390.744850
+1774.742268 1 588.145493 30.487710 21800.531169
+2216.272824 0 467.093275 31.533770 7991.454914
+1526.314833 0 433.703579 32.809884 11209.353005
+2089.304275 0 321.084472 31.021678 27856.425967
+3046.990363 1 460.307711 32.631130 7182.074497
+1317.038946 0 593.108863 31.275244 26560.416891
+1352.101840 0 306.529955 30.979817 21571.506540
+3163.866962 0 304.262526 30.303302 24717.381911
+3417.487105 1 367.150160 31.345893 13987.483107
+4582.950578 0 478.091198 32.280854 23336.546777
+1183.214096 1 435.761798 31.182880 32631.545653
+549.532217 0 436.107384 32.534391 23409.562965
+4018.079504 1 514.150908 32.582211 11007.862101
+4537.862661 0 474.720396 31.380999 6160.620917
+2023.882772 1 576.659734 31.873967 24564.400593
+301.339657 0 525.192309 30.612934 19844.142010
+989.521636 0 323.610331 30.471812 6456.387005
+3647.414661 1 353.160773 30.687803 20893.860619
+184.390450 0 367.556798 32.894773 7004.911101
+2851.709060 1 420.820133 32.238140 4715.183972
+3988.318894 0 352.860472 31.705620 18426.123706
+2244.796525 0 332.927674 31.617697 20768.384056
+726.779082 0 370.641547 32.668547 5629.086601
+2379.683675 1 566.970531 30.287987 26943.017249
+2039.968848 0 550.042019 31.203815 6686.226792
+949.545757 1 302.786734 30.845201 7393.079915
+2509.519553 1 378.264518 31.980310 16762.542739
+2949.998835 0 472.701800 32.896473 12708.457939
+3958.575946 1 402.750276 30.318068 10667.350272
+550.167142 0 394.589505 30.431805 32262.087079
+1462.278890 1 534.727448 31.514861 27759.667272
+2696.454522 1 535.656182 32.096399 32431.581354
+3756.029506 1 332.505291 32.486460 26887.712404
+1948.023625 1 470.099353 31.100519 31437.137850
+2008.758594 1 489.846790 31.656842 18513.611970
+889.052994 0 378.045215 30.360086 30999.351259
+956.860624 0 498.874188 31.238091 32986.760129
+1596.180150 1 468.477647 32.876509 10708.862352
+2044.623310 0 341.364451 30.831209 4377.873160
+1670.382802 0 344.655922 31.528413 25861.700033
+2751.419144 0 525.374048 31.427506 23893.998314
+1139.794772 0 510.822662 32.192052 27853.949410
+1951.275598 1 355.760660 30.362940 26312.111688
+3652.432308 0 330.985338 32.585428 25167.952966
+500.427055 1 453.806271 30.774837 16661.345691
+2884.409660 0 529.815414 31.683479 3317.970560
+2365.322179 0 540.382433 31.999929 29498.364669
+3853.672492 1 367.171163 32.671321 7512.370365
+938.719173 0 569.038125 31.387703 23444.121607
+3058.327259 1 411.675057 31.700811 28738.835254
+4215.442937 0 425.942827 31.617841 19505.555179
+1639.143130 1 310.802062 31.996376 24378.789768
+3015.816015 1 397.865213 30.774914 32208.055822
+1361.137094 1 314.093792 30.698064 21667.130306
+1900.651819 1 582.449408 30.643572 26755.743826
+2038.834435 1 359.535511 31.758170 26389.582514
+4217.395801 1 532.389233 32.367523 28128.901119
+4048.824843 1 415.408543 30.661216 6695.688224
+3278.940683 0 590.282168 31.148497 17706.150830
+2590.370816 1 467.475790 32.894899 31984.209283
+2961.953200 1 409.350826 30.412768 26800.587858
+757.305077 0 531.036597 32.692111 23056.009036
+4823.679979 0 547.031467 31.700254 4679.396986
+2159.887871 1 461.722894 31.554845 23118.215328
+1990.158592 0 436.583574 32.276433 22516.752446
+4010.528172 0 315.169873 31.546055 4521.472077
+527.527936 0 429.812828 32.161079 26933.467576
+2517.395024 1 442.643281 32.903345 11566.976782
+3973.530406 0 444.298671 30.860264 30738.467879
+909.115735 1 548.278268 31.969666 3061.067051
+1003.138129 0 599.492409 32.731826 5268.191640
+496.214281 1 525.539783 32.288968 18314.877509
+3321.884314 1 518.880235 31.593901 27585.403619
+531.060305 0 311.564795 30.660864 26630.390327
+1985.865186 0 468.954849 30.022861 7126.223629
+1075.080978 0 566.983564 32.548614 4189.576013
+2899.914447 0 501.981481 31.522821 3699.352956
+1209.686792 1 397.096875 30.589148 26885.709293
+1061.484419 0 399.802412 31.376845 21631.836258
+2341.148192 0 599.041273 32.705976 7273.381892
+2350.754927 1 578.118121 32.505017 16776.353449
+3715.290307 1 374.328239 31.642496 24018.734911
+353.266473 1 306.976330 30.935877 17228.695180
+739.123546 1 529.022431 30.659845 18290.776722
+465.537407 0 308.998940 31.824504 3898.463033
+3524.840397 1 385.275512 30.348328 6065.321632
+751.616003 0 322.281940 30.576017 20828.524508
+4989.405682 1 381.074571 30.437314 16844.697669
+488.563077 0 318.125000 31.138994 28854.306278
+4558.130890 0 479.058472 31.076337 3808.209751
+3644.754920 0 437.984557 30.038643 24302.584238
+4595.447449 0 323.287824 31.904699 24911.538572
+1064.168831 1 541.625023 30.168168 12380.733046
+4197.068049 0 418.764938 31.133507 14481.902057
+3593.303853 0 545.685419 30.904012 27860.611155
+3842.381434 0 596.665096 31.650039 32071.065698
+781.763222 0 544.786508 32.895055 4905.803177
+3941.386557 1 362.501808 32.489031 8386.976863
+4285.917626 1 373.105515 31.909823 18733.487831
+4152.435642 1 440.036083 30.439913 4669.690889
+2326.656610 0 409.061810 31.219523 17268.378337
+887.243489 1 458.991888 32.904742 32560.558550
+2383.154926 0 402.990252 31.501777 32165.195285
+1575.020117 1 567.576879 31.223480 11851.574021
+2681.793474 0 382.275765 30.672297 31505.220294
+3745.755933 1 524.831384 31.755858 20277.860327
+2310.372497 0 381.887619 32.549019 23884.208743
+1934.629896 0 430.992059 30.090286 16064.349955
+228.815046 1 510.344577 30.272399 13314.791074
+106.494754 0 421.895617 30.249703 6022.257615
+107.758807 1 352.048161 31.015999 27010.353099
+2366.497161 1 589.257754 32.599191 29060.870687
+3652.314844 1 325.392214 30.378902 9562.454775
+4628.382437 0 389.680904 30.488690 21333.408550
+3992.522719 1 460.945057 30.259246 21474.935828
+3116.879118 0 470.551054 30.384387 23511.934422
+3669.891947 1 468.454570 30.792238 6696.682439
+2455.615625 1 536.995869 30.130133 27385.537567
+334.254171 0 510.509005 31.195189 15166.910603
+2392.466005 1 503.464732 31.609503 28757.028622
+563.741633 0 535.281334 31.242527 11455.193241
+2683.204132 0 530.400221 30.122120 21995.611523
+116.891969 0 458.723324 30.381523 5852.247596
+2817.127744 1 547.927894 32.877328 22278.789600
+323.704077 0 432.630710 31.643437 8793.259513
+2223.194666 1 584.593612 31.798636 4232.455807
+4836.806252 0 328.278741 31.551587 28407.563886
+2302.391657 1 304.186017 30.683099 21760.214902
+1597.778913 1 593.734794 32.275770 26311.401980
+3231.228313 0 469.955959 31.155902 21789.637085
+4329.078171 1 326.299894 30.629070 28598.928501
+3004.531335 0 345.517768 31.874809 20075.187439
+3426.658667 0 407.558711 32.882949 32961.004750
+2942.293151 1 474.436479 30.681717 10297.918976
+2790.908728 1 372.634696 31.883438 18617.426989
+3645.736003 0 338.768900 30.271779 13340.959790
+4954.230936 0 475.204218 32.305344 31368.726270
+1374.169666 0 322.794825 32.054615 29507.251607
+952.311143 0 466.066676 31.283949 31527.207980
+1582.935333 1 430.497522 32.636945 6365.677867
+1574.488078 1 512.547475 32.000036 3214.800133
+3675.314941 0 420.357568 30.365188 24034.010725
+3737.900636 0 310.463467 32.944674 19212.098849
+453.453905 1 307.998144 31.734006 6696.312537
+180.430188 0 500.121289 31.647589 25439.918801
+2223.386978 0 342.168201 31.413275 10321.858985
+3347.152707 0 450.255234 32.419409 21979.139997
+2032.988552 1 359.268392 32.444386 4151.429949
+1855.874006 1 486.902129 32.194183 28734.066535
+576.472090 0 491.131111 30.250436 21183.624322
+2629.739140 1 490.880496 30.706030 20741.966223
+3234.449151 0 591.702238 31.565719 9467.501732
+4568.779861 0 391.821980 31.778473 12171.711242
+1369.135728 1 352.424832 32.501343 24914.054783
+1203.818730 1 453.408627 30.021814 9344.375802
+1328.193901 1 326.714975 31.957412 22506.010391
+3143.392191 0 488.722719 30.906428 20404.329854
+1516.300551 1 524.551137 30.783381 26784.582992
+4401.773609 1 490.129039 31.474209 27673.989380
+471.342147 0 345.653763 30.714091 27614.068357
+2846.556659 0 387.773914 32.101213 25758.358083
+3113.332960 0 340.017571 31.250862 8517.797718
+569.863929 1 507.285752 32.723247 7350.729977
+1351.433366 1 355.255156 31.629009 10035.728199
+2057.479962 0 504.044203 31.725231 14699.764489
+2029.207221 1 326.501058 31.883646 25956.067323
+4987.734533 0 391.140321 30.918262 15723.785267
+2728.957211 1 498.936107 32.017269 21531.521853
+2885.926679 0 458.513681 31.221876 23111.352960
+1214.796671 0 531.593343 32.409879 26260.051485
+4073.228494 1 563.045005 31.629048 5742.795393
+632.393708 1 428.087852 32.360714 20236.813182
+4548.671698 0 387.645357 30.048297 28095.140888
+735.589017 0 510.876342 30.202972 28887.002260
+2683.928468 1 352.751328 31.385745 3968.105767
+3761.692766 1 470.650685 31.713693 21645.943688
+3949.678238 0 572.913454 32.230970 18919.434641
+1695.346973 0 387.076598 32.045632 8265.498458
+4805.534106 0 322.040787 30.256378 20922.671974
+3834.130476 0 548.536590 30.581470 12530.815196
+3700.095266 0 434.154537 30.131024 28914.087338
+2083.132166 1 435.468758 31.961937 24155.150686
+2566.436697 1 455.368285 30.113525 23903.526223
+3372.416477 1 439.390127 32.380970 8312.312507
+3341.243804 1 496.573292 32.055384 9484.908468
+1289.638171 1 346.709851 31.116968 13853.905474
+441.798476 1 395.115019 32.820316 27178.509065
+3876.357092 1 501.568592 30.100469 3872.161975
+3082.204963 0 445.140351 31.403349 21141.060076
+1222.633709 1 510.019589 32.832816 12522.872061
+3257.092925 1 567.760441 30.935411 17780.698060
+4978.318869 0 488.434181 31.349966 20399.317688
+268.617967 0 563.410236 30.677835 13353.114481
+3668.370880 1 432.604476 31.336602 18298.862821
+1616.920156 1 528.362539 30.664065 18290.365790
+1582.470140 1 566.928554 31.829775 31526.762426
+1894.952274 0 576.611659 32.930078 5451.111476
+130.599064 1 357.672426 30.192175 5127.046275
+1658.869611 0 518.418202 30.222626 21787.150970
+1930.327624 0 539.744396 31.409771 11998.451350
+1360.478477 1 357.741711 32.358872 4894.722814
+2762.405397 1 498.202228 32.110196 20926.903902
+66.089100 0 525.721147 30.742496 28773.910408
+4617.260796 0 388.807122 31.052491 19010.672351
+4046.674710 1 363.050835 32.600458 8208.205191
+103.101768 0 401.874765 30.885852 29276.307064
+2870.170146 1 479.961929 31.997551 19841.555525
+1069.648823 1 598.655109 31.936390 10655.006702
+1988.650561 0 483.709224 31.937744 11014.379454
+2983.772162 0 328.380916 32.023817 8974.076899
+4534.002381 0 311.260565 30.136494 25947.551507
+3064.855784 1 369.080218 31.164070 29129.542680
+1746.878916 0 313.512424 30.992758 20731.078669
+4018.774695 1 497.866125 31.834090 3971.821700
+3742.871097 1 353.541863 30.115952 20815.885581
+2828.315773 1 454.655490 32.424599 26891.927439
+637.033003 1 424.553960 31.428681 20585.183759
+3087.856465 1 536.195443 31.597533 12661.823700
+604.873967 1 378.515493 30.925004 22439.100233
+3748.531033 1 515.115463 32.625710 6408.313067
+1558.535013 1 433.937532 32.380082 8383.657584
+2365.992186 0 409.216710 31.115409 32511.614098
+1183.022513 0 362.544664 30.959448 5065.906749
+3789.375438 0 595.724663 30.914010 20244.524890
+850.273991 1 387.136972 31.618840 5094.567581
+1231.544018 1 322.683917 31.457545 20767.599384
+1581.181089 1 305.296620 31.856051 17539.248499
+4880.147078 1 316.598007 31.317619 11086.097072
+4013.027211 0 485.648876 31.953258 30138.577997
+31.538821 0 486.966298 30.424936 3940.149283
+1735.735564 1 453.551097 30.972695 29273.785203
+4682.372724 1 426.064369 30.915967 29848.714453
+1419.673303 1 410.849400 31.957769 20423.937410
+4735.693777 0 331.994083 30.314568 21510.813276
+796.562868 0 502.136561 31.584336 14765.778268
+2637.389588 1 550.253910 31.759164 29206.835372
+1596.336472 0 367.679269 32.352327 5586.166324
+4673.788569 0 380.575280 32.897426 12723.911545
+2339.522676 0 405.444884 30.501805 20468.692308
+2935.072522 0 458.797029 30.585909 9955.177397
+2865.796979 0 497.765242 31.556959 27770.325805
+1032.502528 1 390.204545 31.050812 26330.543896
+1571.491134 1 454.897576 32.799730 4600.369582
+453.244362 0 513.834562 31.693168 18448.030313
+2566.490430 1 370.231825 31.371622 3530.831525
+2470.800433 1 506.755104 32.707049 15365.573438
+1674.578908 1 484.850818 31.657126 31034.118967
+1569.707807 0 380.339060 31.625216 26765.694375
+2048.318765 1 502.814703 32.061583 26794.094870
+4878.355942 1 308.718223 31.785469 9019.516491
+3394.948104 0 354.869800 31.553766 18461.635011
+3784.598905 0 506.145921 32.828899 22372.102070
+3648.947834 0 503.931222 32.468882 27102.757907
+1722.690965 1 389.086536 31.876640 5309.945334
+1266.660624 0 452.533007 30.468222 28572.903114
+794.701939 1 459.408945 32.075307 21311.040821
+1461.841041 0 370.527899 32.775555 10211.935095
+1976.852418 0 492.898069 30.817372 17854.725954
+859.962977 1 304.267257 30.348074 3439.856318
+2900.907710 1 396.919281 31.272958 4858.314485
+3693.615130 1 597.052502 30.903533 25716.531228
+1281.197465 1 470.854828 30.164368 9315.105349
+126.181622 1 381.237663 30.660067 28411.425041
+2517.244075 0 559.007071 30.373940 11135.538898
+2069.720704 1 544.708646 31.610539 16441.322232
+578.677623 1 591.867328 32.146388 21419.935158
+1918.809273 1 308.991680 31.457030 31410.129144
+1074.094412 0 336.860228 30.308287 17215.627286
+3034.994121 0 560.991229 31.862778 17867.418001
+4794.119298 1 480.292221 30.685158 11676.367480
+3307.911419 0 389.291016 30.558138 31108.361766
+457.891802 1 588.926068 32.330020 21119.054157
+1258.403059 0 332.016717 31.309245 16741.712001
+3312.429446 0 568.695418 31.010073 14155.107714
+3839.897879 0 370.810835 31.722406 16960.037012
+2068.214455 0 392.107828 30.058746 32853.457910
+2526.811696 0 510.287902 30.913056 22333.505644
+914.090787 1 323.087744 32.697177 14094.308393
+1890.199376 0 597.011219 31.259415 19326.502842
+1234.225123 0 536.105198 30.519191 9869.430841
+350.596239 1 320.162343 30.373749 5065.307636
+2846.841176 0 403.168597 31.624641 25354.766351
+1853.530932 1 564.965357 30.095729 11551.597590
+1183.876207 1 572.411649 32.013383 17167.971369
+3953.600477 0 404.253211 30.364392 17515.400202
+3289.696058 1 576.771217 31.015263 14654.467190
+2975.113909 1 587.011738 31.223106 7437.225919
+1404.141852 0 372.793625 30.325109 13500.525950
+1634.182114 0 569.705673 32.338410 27233.776345
+462.646125 1 467.892571 30.013450 27556.195432
+4138.443866 0 318.632678 30.348871 21795.067964
+4139.367443 1 304.173328 32.473154 29126.297266
+969.741373 1 599.886523 32.409973 16299.878199
+1494.463287 1 429.995117 30.484126 17461.996880
+4336.770459 1 416.241963 30.152207 32032.419174
+3064.782012 1 428.313584 31.688383 31731.608033
+2581.013362 0 490.480997 32.990983 17803.475008
+413.230864 1 430.932119 30.358681 18338.089371
+4469.602591 0 544.949864 30.599950 17579.676286
+2299.774314 1 532.263729 32.087371 25491.390739
+3570.957635 0 429.164688 30.798996 21783.084247
+2730.725224 0 434.512962 32.353024 8584.173078
+1859.613968 0 465.494240 30.853468 24199.875583
+3239.212454 1 585.752947 31.392592 17429.754381
+1626.243010 0 548.642457 32.586064 28320.002213
+4552.584886 0 411.908494 32.541147 18556.711124
+719.687465 1 437.155107 32.439905 18592.734838
+4528.629292 1 598.029535 30.261328 28073.717443
+1445.996134 1 407.243137 30.545202 28549.910302
+2143.373187 1 581.134422 30.728743 26291.822745
+4543.013656 0 406.670256 30.716994 24689.785142
+3479.794580 0 585.305633 31.528087 18146.824396
+2206.562888 0 526.292818 32.250928 6586.281804
+1796.678722 0 441.391846 32.716749 16294.189483
+2705.347154 0 551.894682 31.578422 11390.827042
+172.103511 1 320.919515 31.863888 29876.119334
+315.943273 1 512.197902 32.302445 17856.311999
+2750.695088 1 565.301177 31.514433 19939.221077
+4398.626158 0 418.977196 32.337330 4024.240756
+1522.050792 0 554.491924 32.249187 30588.435646
+3385.739992 1 496.030916 31.480018 17815.393181
+955.825566 1 393.353810 31.434487 26026.721071
+3971.619403 1 552.518065 30.319670 5068.002389
+2010.848206 1 569.817785 31.925723 26117.233780
+2181.752401 1 535.925058 30.952605 3509.328188
+905.194864 1 588.512406 30.229715 15643.423750
+3274.240917 0 590.044688 32.642379 13352.913630
+2061.196669 0 457.883033 30.745459 19655.487967
+4056.062516 0 354.337157 30.418176 24440.653038
+611.999600 0 432.598681 30.201091 19387.208118
+45.061679 1 437.997979 31.939322 21882.265038
+2091.842258 1 445.122275 32.963802 4496.459612
+2785.821696 1 535.953261 32.295237 5166.117813
+1272.363761 0 498.289820 31.101707 13536.602174
+2161.870929 0 487.379529 30.927239 22791.304718
+1230.121506 0 326.397921 31.953971 21880.979367
+2065.316190 0 405.291555 30.660933 28025.816633
+3230.173916 1 365.768097 30.055539 23327.051800
+4572.451029 1 323.399890 31.386897 21097.082477
+4916.172714 0 384.428209 32.754127 10981.900793
+4882.524497 0 399.434058 30.795923 18117.477126
+1457.322605 1 412.012101 30.834341 20230.063544
+1988.228627 0 566.906278 31.371606 15413.227181
+4337.262787 1 316.056609 32.556132 11299.235579
+2762.610907 1 454.177332 32.779747 14436.150973
+4756.268480 1 562.493863 32.727774 26209.569931
+3031.941096 0 313.505196 31.386484 9197.983181
+4586.200757 0 329.983728 30.115395 27039.250998
+300.892434 0 520.668472 32.713615 14579.471742
+392.962829 0 513.293270 30.568538 10311.182381
+2116.838847 1 355.318371 30.883049 18887.246656
+779.846218 1 407.952315 30.145534 22063.619725
+4236.655978 0 326.344600 32.797447 16004.751246
+4008.051049 1 319.382840 30.262423 20092.936132
+1550.891721 1 570.189663 32.238898 11686.528747
+2616.214276 1 448.726264 31.633295 26368.520969
+2775.847491 1 532.436984 30.938190 21751.684582
+
diff -r 7d66bfa4fd56 -r 21b12c7c52e4 test-data/coxlassotest.html
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/coxlassotest.html Sat Oct 31 01:34:27 2015 -0400
@@ -0,0 +1,167 @@
+
+
+
+
+
+
+
+
+
+
+
Galaxy Tool "rglasso" run at 18/02/2015 22:06:08
+
cox images and outputs
+(Click on a thumbnail image to download the corresponding original PDF image)
+
+
+
+
+
+
+
+
+
+
cox log output
+
+
+
+[1] @@@ Cox model will be predicted at times = 3.724234,1325.40516225,2373.090418,3650.1730825,4989.405682
+
+[1] Times each predictor was selected in CV models (excluding zero count predictors)
+
+ hits_lambda_1se hits_lambda_min
+
+rhubarb 10 10
+
+vegemite 6 6
+
+[1] @@@ Results preview:
+
+ regulator partial_likelihood forced_in glmnet_model best_lambda lambdaChoice alpha
+
+rhubarb rhubarb 0.0012215 TRUE cox 0.04197 lambda.1se 0.95
+
+
+
+
+
rglasso log output
+
+
+
+Loading required package: Matrix
+
+Loading required package: methods
+
+Loaded glmnet 1.9-8
+
+Loaded lars 1.2
+
+Loading required package: splines
+
+Warning messages:
+
+1: In if (is.na(predict_at)) { :
+
+ the condition has length > 1 and only the first element will be used
+
+2: In if (class(p) == "try-error") { :
+
+ the condition has length > 1 and only the first element will be used
+
+3: In plot.window(...) : "label" is not a graphical parameter
+
+4: In plot.xy(xy, type, ...) : "label" is not a graphical parameter
+
+5: In axis(side = side, at = at, labels = labels, ...) :
+
+ "label" is not a graphical parameter
+
+6: In axis(side = side, at = at, labels = labels, ...) :
+
+ "label" is not a graphical parameter
+
+7: In box(...) : "label" is not a graphical parameter
+
+8: In title(...) : "label" is not a graphical parameter
+
+
+
+
+
Other log output
+
+
+
+## Toolfactory running rglasso as Rscript script
+
+[1] @@@ Using alpha = 0.95 for all models
+
+[1] @@@@ Input has 500 samples and 3 predictors
+
+[1] @@@ Cox model will predict yvar= status using cols= rhubarb,vegemite,apple n preds= 3 forced in= rhubarb
+
+[1] @@@ SessionInfo for this R session:
+
+R version 3.1.0 (2014-04-10)
+
+Platform: x86_64-unknown-linux-gnu (64-bit)
+
+locale:
+
+ [1] LC_CTYPE=en_AU.UTF-8 LC_NUMERIC=C LC_TIME=en_AU.UTF-8 LC_COLLATE=en_AU.UTF-8 LC_MONETARY=en_AU.UTF-8 LC_MESSAGES=en_AU.UTF-8 LC_PAPER=en_AU.UTF-8 LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=en_AU.UTF-8 LC_IDENTIFICATION=C
+
+attached base packages:
+
+[1] splines methods stats graphics grDevices utils datasets base
+
+other attached packages:
+
+[1] pec_2.4.4 survival_2.37-7 lars_1.2 glmnet_1.9-8 Matrix_1.1-5
+
+loaded via a namespace (and not attached):
+
+[1] codetools_0.2-10 foreach_1.4.2 grid_3.1.0 iterators_1.0.7 lattice_0.20-29 lava_1.3 prodlim_1.5.1
+
+Warning messages:
+
+1: In if (is.na(predict_at)) { ... :
+
+ the condition has length > 1 and only the first element will be used
+
+2: In if (class(p) == "try-error") { ... :
+
+ the condition has length > 1 and only the first element will be used
+
+3: In plot.window(...) : "label" is not a graphical parameter
+
+4: In plot.xy(xy, type, ...) : "label" is not a graphical parameter
+
+5: In axis(side = side, at = at, labels = labels, ...) :
+
+ "label" is not a graphical parameter
+
+6: In axis(side = side, at = at, labels = labels, ...) :
+
+ "label" is not a graphical parameter
+
+7: In box(...) : "label" is not a graphical parameter
+
+8: In title(...) : "label" is not a graphical parameter
+
+
+