mirror of
https://github.com/galaxyproject/galaxy.git
synced 2026-09-24 16:30:27 +08:00
Added 3 tools: 1) to compute weighted average values based on interval intersections (under Regional variation), 2) to perform logistic regression (under Multiple regression), 3) to compute partial R-squared (under Multiple regression).
This commit is contained in:
Guruprasad Anada
@@ -175,6 +175,7 @@
|
||||
<section name="Regional Variation" id="regVar">
|
||||
<tool file="regVariation/windowSplitter.xml" />
|
||||
<tool file="regVariation/featureCounter.xml" />
|
||||
<tool file="regVariation/WeightedAverage.xml" />
|
||||
<tool file="regVariation/quality_filter.xml" />
|
||||
<tool file="regVariation/maf_cpg_filter.xml" />
|
||||
<tool file="regVariation/getIndels_2way.xml" />
|
||||
@@ -194,8 +195,10 @@
|
||||
</section>
|
||||
<section name="Multiple regression" id="multReg">
|
||||
<tool file="regVariation/linear_regression.xml" />
|
||||
<tool file="regVariation/logistic_regression_vif.xml" />
|
||||
<tool file="regVariation/best_regression_subsets.xml" />
|
||||
<tool file="regVariation/rcve.xml" />
|
||||
<tool file="regVariation/partialR_square.xml" />
|
||||
</section>
|
||||
<section name="Multivariate Analysis" id="multVar">
|
||||
<tool file="multivariate_stats/pca.xml" />
|
||||
|
||||
Executable
+103
@@ -0,0 +1,103 @@
|
||||
#!/usr/bin/env python
|
||||
"""
|
||||
|
||||
usage: %prog bed_file_1 bed_file_2 out_file
|
||||
-1, --cols1=N,N,N,N: Columns for chr, start, end, strand in first file
|
||||
-2, --cols2=N,N,N,N,N: Columns for chr, start, end, strand, name/value in second file
|
||||
"""
|
||||
from galaxy import eggs
|
||||
|
||||
import collections
|
||||
import sys, string
|
||||
#import numpy
|
||||
from galaxy import eggs
|
||||
import pkg_resources
|
||||
pkg_resources.require( "bx-python" )
|
||||
import sys, traceback, fileinput
|
||||
from warnings import warn
|
||||
from galaxy.tools.util.galaxyops import *
|
||||
from bx.cookbook import doc_optparse
|
||||
|
||||
|
||||
#export PYTHONPATH=~/galaxy/lib/
|
||||
#running command python WeightedAverage.py interval_interpolate.bed value_interpolate.bed interpolate_result.bed
|
||||
|
||||
def stop_err(msg):
|
||||
sys.stderr.write(msg)
|
||||
sys.exit()
|
||||
|
||||
|
||||
def FindRate(chromosome,start_stop,dictType):
|
||||
OverlapList=[]
|
||||
for tempO in dictType[chromosome]:
|
||||
DatabaseInterval=[tempO[0],tempO[1]]
|
||||
Overlap=GetOverlap(start_stop,DatabaseInterval)
|
||||
if Overlap>0:
|
||||
OverlapList.append([Overlap,tempO[2]])
|
||||
|
||||
if len(OverlapList)>0:
|
||||
|
||||
SumRecomb=0
|
||||
SumOverlap=0
|
||||
for member in OverlapList:
|
||||
SumRecomb+=member[0]*member[1]
|
||||
SumOverlap+=member[0]
|
||||
averageRate=SumRecomb/SumOverlap
|
||||
|
||||
return averageRate
|
||||
|
||||
else:
|
||||
return 'NA'
|
||||
|
||||
|
||||
|
||||
def GetOverlap(a,b):
|
||||
return min(a[1],b[1])-max(a[0],b[0])
|
||||
|
||||
options, args = doc_optparse.parse( __doc__ )
|
||||
|
||||
try:
|
||||
chr_col_1, start_col_1, end_col_1, strand_col1 = parse_cols_arg( options.cols1 )
|
||||
chr_col_2, start_col_2, end_col_2, strand_col2, name_col_2 = parse_cols_arg( options.cols2 )
|
||||
input1, input2, input3 = args
|
||||
except Exception, eee:
|
||||
print eee
|
||||
stop_err( "Data issue: click the pencil icon in the history item to correct the metadata attributes." )
|
||||
|
||||
|
||||
|
||||
fd2=open(input2)
|
||||
lines2=fd2.readlines()
|
||||
RecombChrDict=collections.defaultdict(list)
|
||||
|
||||
skipped=0
|
||||
for line in lines2:
|
||||
temp=line.strip().split()
|
||||
try:
|
||||
assert float(temp[int(name_col_2)])
|
||||
except:
|
||||
skipped+=1
|
||||
continue
|
||||
tempIndex=[int(temp[int(start_col_2)]),int(temp[int(end_col_2)]),float(temp[int(name_col_2)])]
|
||||
RecombChrDict[temp[int(chr_col_2)]].append(tempIndex)
|
||||
|
||||
print "Skipped %d features with invalid values" %(skipped)
|
||||
|
||||
fd1=open(input1)
|
||||
lines=fd1.readlines()
|
||||
finalProduct=''
|
||||
for line in lines:
|
||||
temp=line.strip().split('\t')
|
||||
chromosome=temp[int(chr_col_1)]
|
||||
start=int(temp[int(start_col_1)])
|
||||
stop=int(temp[int(end_col_1)])
|
||||
start_stop=[start,stop]
|
||||
RecombRate=FindRate(chromosome,start_stop,RecombChrDict)
|
||||
try:
|
||||
RecombRate="%.4f" %(float(RecombRate))
|
||||
except:
|
||||
RecombRate=RecombRate
|
||||
finalProduct+=line.strip()+'\t'+str(RecombRate)+'\n'
|
||||
fdd=open(input3,'w')
|
||||
fdd.writelines(finalProduct)
|
||||
fdd.close()
|
||||
Executable
+71
@@ -0,0 +1,71 @@
|
||||
<tool id="wtavg" name="Assign weighted-average" version="1.0.0">
|
||||
<description> of the values of features overlapping an interval </description>
|
||||
<command interpreter="python">WeightedAverage.py $genomic_interval $genomic_feature $out_file1 -1 ${genomic_interval.metadata.chromCol},${genomic_interval.metadata.startCol},${genomic_interval.metadata.endCol},${genomic_interval.metadata.strandCol} -2 ${genomic_feature.metadata.chromCol},${genomic_feature.metadata.startCol},${genomic_feature.metadata.endCol},${genomic_feature.metadata.strandCol},${genomic_feature.metadata.nameCol}</command>
|
||||
|
||||
<inputs>
|
||||
<param format="interval" name="genomic_interval" type="data" label="Genomic intervals (first dataset)" help="Dataset missing? See Note below."/>
|
||||
<param format="interval" name="genomic_feature" label="Genomic features (second dataset)" type="data" help="Make sure the value column is specified. See Note below." />
|
||||
</inputs>
|
||||
<outputs>
|
||||
<data format="input" name="out_file1" metadata_source="genomic_interval" />
|
||||
</outputs>
|
||||
<tests>
|
||||
<!-- Test data with valid values -->
|
||||
<test>
|
||||
<param name="genomic_interval" value="interval_interpolate.bed"/>
|
||||
<param name="genomic_feature" value="value_interpolate.bed"/>
|
||||
<output name="out_file1" file="interpolate_result.bed"/>
|
||||
</test>
|
||||
|
||||
</tests>
|
||||
<help>
|
||||
|
||||
|
||||
.. class:: infomark
|
||||
|
||||
**What it does**
|
||||
|
||||
For each interval in your first dataset, this tool calculates the weighted average value of the overlapping features in your second dataset.
|
||||
|
||||
- When a genomic interval partially or totally overlaps a single genomic feature, the value of that genomic feature is assigned to the genomic interval.
|
||||
- When a genomic interval partially or totally overlaps with more than one genomic features, the average of the values of the overlapping genomic features weighted by the corresponding number of overlapping bases is assigned to the genomic interval.
|
||||
- When a genomic interval does not overlap with any genomic feature, 'NA' will be assigned as it's value.
|
||||
|
||||
-----
|
||||
|
||||
.. class:: warningmark
|
||||
|
||||
**Note**
|
||||
|
||||
The input datasets should be in **bed** or **interval** format. Please use "edit attributes"/pencil icon to specify the column containing the values for the features in the second dataset as **name/identifier** column.
|
||||
|
||||
The output will contain all the columns in the first input plus a new column containing the assigned value for each interval.
|
||||
|
||||
-----
|
||||
|
||||
**Example**
|
||||
|
||||
- Suppose our first dataset contains the following **genomic intervals**::
|
||||
|
||||
chr start stop
|
||||
chr1 1000 2000
|
||||
chr1 3000 5000
|
||||
chr1 8000 9000
|
||||
|
||||
- and our second dataset contains the following **genomic features** each having an associated value (in fourth column) ::
|
||||
|
||||
chr start stop name
|
||||
chr1 900 1200 0.5
|
||||
chr1 2900 3100 0.2
|
||||
chr1 4800 5100 0.8
|
||||
|
||||
- For each **genomic interval** in our first dataset, this tool calculates the weighted average value of the overlapping **genomic features** in our second dataset ::
|
||||
|
||||
chr1 1000 2000 0.5
|
||||
chr1 3000 5000 0.6
|
||||
chr1 8000 9000 NA
|
||||
|
||||
|
||||
|
||||
</help>
|
||||
</tool>
|
||||
+185
@@ -0,0 +1,185 @@
|
||||
#!/usr/bin/env python
|
||||
|
||||
from galaxy import eggs
|
||||
import sys, string
|
||||
from rpy import *
|
||||
import numpy
|
||||
|
||||
#export PYTHONPATH=~/galaxy/lib/
|
||||
|
||||
def stop_err(msg):
|
||||
sys.stderr.write(msg)
|
||||
sys.exit()
|
||||
#infile = 'logreg_inp.tab'
|
||||
#y_col=3
|
||||
#x_cols=[1,2,3]
|
||||
#outfile='logreg_out.txt'
|
||||
#python logistic_regression_vif.py logreg_inp.tab 4 1,2,3 logreg_out2.tabular # running test
|
||||
infile = sys.argv[1]
|
||||
y_col = int(sys.argv[2])-1
|
||||
x_cols = sys.argv[3].split(',')
|
||||
outfile = sys.argv[4]
|
||||
|
||||
|
||||
print "Predictor columns: %s; Response column: %d" %(x_cols,y_col+1)
|
||||
fout = open(outfile,'w')
|
||||
elems = []
|
||||
for i, line in enumerate( file ( infile )):
|
||||
line = line.rstrip('\r\n')
|
||||
if len( line )>0 and not line.startswith( '#' ):
|
||||
elems = line.split( '\t' )
|
||||
break
|
||||
if i == 30:
|
||||
break # Hopefully we'll never get here...
|
||||
|
||||
if len( elems )<1:
|
||||
stop_err( "The data in your input dataset is either missing or not formatted properly." )
|
||||
|
||||
y_vals = []
|
||||
x_vals = []
|
||||
|
||||
for k,col in enumerate(x_cols):
|
||||
x_cols[k] = int(col)-1
|
||||
x_vals.append([])
|
||||
|
||||
NA = 'NA'
|
||||
for ind,line in enumerate( file( infile )):
|
||||
if line and not line.startswith( '#' ):
|
||||
try:
|
||||
fields = line.split("\t")
|
||||
try:
|
||||
yval = float(fields[y_col])
|
||||
except:
|
||||
yval = r('NA')
|
||||
y_vals.append(yval)
|
||||
for k,col in enumerate(x_cols):
|
||||
try:
|
||||
xval = float(fields[col])
|
||||
except:
|
||||
xval = r('NA')
|
||||
x_vals[k].append(xval)
|
||||
except:
|
||||
pass
|
||||
|
||||
x_vals1 = numpy.asarray(x_vals).transpose()
|
||||
|
||||
check1=0
|
||||
check0=0
|
||||
for i in y_vals:
|
||||
if i == 1:
|
||||
check1=1
|
||||
if i == 0:
|
||||
check0=1
|
||||
if check1==0 or check0==0:
|
||||
sys.exit("Warning: logistic regression must have at least two classes")
|
||||
|
||||
for i in y_vals:
|
||||
if i not in [1,0,r('NA')]:
|
||||
print >>fout, str(i)
|
||||
sys.exit("Warning: the current version of this tool can run only with two classes and need to be labeled as 0 and 1.")
|
||||
|
||||
|
||||
dat= r.list(x=array(x_vals1), y=y_vals)
|
||||
novif=0
|
||||
set_default_mode(NO_CONVERSION)
|
||||
try:
|
||||
linear_model = r.glm(r("y ~ x"), data = r.na_exclude(dat),family="binomial")
|
||||
#r('library(car)')
|
||||
#r.assign('dat',dat)
|
||||
#r.assign('ncols',len(x_cols))
|
||||
#r.vif(r('glm(dat$y ~ ., data = na.exclude(data.frame(as.matrix(dat$x,ncol=ncols))->datx),family="binomial")')).as_py()
|
||||
|
||||
except RException, rex:
|
||||
stop_err("Error performing logistic regression on the input data.\nEither the response column or one of the predictor columns contain only non-numeric or invalid values.")
|
||||
if len(x_cols)>1:
|
||||
try:
|
||||
|
||||
r('library(car)')
|
||||
r.assign('dat',dat)
|
||||
r.assign('ncols',len(x_cols))
|
||||
vif=r.vif(r('glm(dat$y ~ ., data = na.exclude(data.frame(as.matrix(dat$x,ncol=ncols))->datx),family="binomial")'))
|
||||
except RException, rex:
|
||||
print rex
|
||||
else:
|
||||
novif=1
|
||||
|
||||
set_default_mode(BASIC_CONVERSION)
|
||||
|
||||
coeffs=linear_model.as_py()['coefficients']
|
||||
null_deviance=linear_model.as_py()['null.deviance']
|
||||
residual_deviance=linear_model.as_py()['deviance']
|
||||
yintercept= coeffs['(Intercept)']
|
||||
summary = r.summary(linear_model)
|
||||
co = summary.get('coefficients', 'NA')
|
||||
"""
|
||||
if len(co) != len(x_vals)+1:
|
||||
stop_err("Stopped performing logistic regression on the input data, since one of the predictor columns contains only non-numeric or invalid values.")
|
||||
"""
|
||||
|
||||
try:
|
||||
yintercept = r.round(float(yintercept), digits=10)
|
||||
pvaly = r.round(float(co[0][3]), digits=10)
|
||||
except:
|
||||
pass
|
||||
print >>fout, "response column\tc%d" %(y_col+1)
|
||||
tempP=[]
|
||||
for i in x_cols:
|
||||
tempP.append('c'+str(i+1))
|
||||
tempP=','.join(tempP)
|
||||
print >>fout, "predictor column(s)\t%s" %(tempP)
|
||||
print >>fout, "Y-intercept\t%s" %(yintercept)
|
||||
print >>fout, "p-value (Y-intercept)\t%s" %(pvaly)
|
||||
|
||||
if len(x_vals) == 1: #Simple linear regression case with 1 predictor variable
|
||||
try:
|
||||
slope = r.round(float(coeffs['x']), digits=10)
|
||||
except:
|
||||
slope = 'NA'
|
||||
try:
|
||||
pval = r.round(float(co[1][3]), digits=10)
|
||||
except:
|
||||
pval = 'NA'
|
||||
print >>fout, "Slope (c%d)\t%s" %(x_cols[0]+1,slope)
|
||||
print >>fout, "p-value (c%d)\t%s" %(x_cols[0]+1,pval)
|
||||
else: #Multiple regression case with >1 predictors
|
||||
ind=1
|
||||
while ind < len(coeffs.keys()):
|
||||
try:
|
||||
slope = r.round(float(coeffs['x'+str(ind)]), digits=10)
|
||||
except:
|
||||
slope = 'NA'
|
||||
print >>fout, "Slope (c%d)\t%s" %(x_cols[ind-1]+1,slope)
|
||||
try:
|
||||
pval = r.round(float(co[ind][3]), digits=10)
|
||||
except:
|
||||
pval = 'NA'
|
||||
print >>fout, "p-value (c%d)\t%s" %(x_cols[ind-1]+1,pval)
|
||||
ind+=1
|
||||
|
||||
rsq = summary.get('r.squared','NA')
|
||||
|
||||
|
||||
try:
|
||||
rsq= r.round(float((null_deviance-residual_deviance)/null_deviance), digits=5)
|
||||
null_deviance= r.round(float(null_deviance), digits=5)
|
||||
residual_deviance= r.round(float(residual_deviance), digits=5)
|
||||
|
||||
#rsq = r.round(float(rsq), digits=5)
|
||||
|
||||
except:
|
||||
pass
|
||||
|
||||
print >>fout, "Null deviance\t%s" %(null_deviance)
|
||||
print >>fout, "Residual deviance\t%s" %(residual_deviance)
|
||||
print >>fout, "pseudo R-squared\t%s" %(rsq)
|
||||
print >>fout, "\n"
|
||||
print >>fout, 'vif'
|
||||
|
||||
if novif==0:
|
||||
py_vif=vif.as_py()
|
||||
count=0
|
||||
for i in sorted(py_vif.keys()):
|
||||
print >>fout,'c'+str(x_cols[count]+1) ,str(py_vif[i])
|
||||
count+=1
|
||||
elif novif==1:
|
||||
print >>fout, "vif can calculate only when model have more than 1 predictor"
|
||||
+74
@@ -0,0 +1,74 @@
|
||||
<tool id="LogisticRegression" name="Perform Logistic Regression with vif" version="1.0.1">
|
||||
<description> </description>
|
||||
<command interpreter="python">
|
||||
logistic_regression_vif.py
|
||||
$input1
|
||||
$response_col
|
||||
$predictor_cols
|
||||
$out_file1
|
||||
1>/dev/null
|
||||
</command>
|
||||
<inputs>
|
||||
<param format="tabular" name="input1" type="data" label="Select data" help="Dataset missing? See TIP below."/>
|
||||
<param name="response_col" label="Response column (Y)" type="data_column" data_ref="input1" numerical="True"/>
|
||||
<param name="predictor_cols" label="Predictor columns (X)" type="data_column" data_ref="input1" numerical="True" multiple="true" >
|
||||
<validator type="no_options" message="Please select at least one column."/>
|
||||
</param>
|
||||
</inputs>
|
||||
<outputs>
|
||||
<data format="input" name="out_file1" metadata_source="input1" />
|
||||
|
||||
</outputs>
|
||||
<requirements>
|
||||
<requirement type="python-module">rpy</requirement>
|
||||
</requirements>
|
||||
<tests>
|
||||
<test>
|
||||
<param name="input1" value="logreg_inp.tabular"/>
|
||||
<param name="response_col" value="4"/>
|
||||
<param name="predictor_cols" value="1,2,3"/>
|
||||
<output name="out_file1" file="logreg_out2.tabular"/>
|
||||
|
||||
</test>
|
||||
</tests>
|
||||
<help>
|
||||
|
||||
|
||||
.. class:: infomark
|
||||
|
||||
**TIP:** If your data is not TAB delimited, use *Edit Datasets->Convert characters*
|
||||
|
||||
-----
|
||||
|
||||
.. class:: infomark
|
||||
|
||||
**What it does**
|
||||
|
||||
This tool uses the **'glm'** function from R statistical package to perform logistic regression on the input data. It outputs one file containing the summary statistics of the performed regression. Also, it calculates VIF(Variance Inflation Factor) with **'vif'** function from library (car) in R.
|
||||
|
||||
|
||||
*R Development Core Team (2010). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL http://www.R-project.org.*
|
||||
|
||||
-----
|
||||
|
||||
.. class:: warningmark
|
||||
|
||||
**Note**
|
||||
|
||||
- This tool currently treats all predictor variables as continuous numeric variables and response variable as categorical variable. Currently, the response variable can have only two classes, namely 0 and 1. The program will take 0 as base class.
|
||||
|
||||
- Rows containing non-numeric (or missing) data in any of the chosen columns will be skipped from the analysis.
|
||||
|
||||
- The summary statistics in the output are described below:
|
||||
|
||||
- Pseudo R-squared: the proportion of model improvement from null model
|
||||
- p-value: p-value for the z-test of the null hypothesis that the corresponding slope is equal to zero against the two-sided alternative.
|
||||
- Coefficient indicates log ratio of (probability to be class 1 / probability to be class 0)
|
||||
|
||||
- This tool also provides **Variance Inflation Factor or VIF** which quantifies the level of multicollinearity. The tool will automatic generate VIF if the model has more than one predictor. The higher the VIF, the higher is the multicollinearity. Multicollinearity will inflate standard error and reduce level of significance of the predictor. In the worst case, it can reverse direction of slope for highly correlated predictors if one of them is significant. A general thumb-rule is to use those predictors having VIF lower than 10 or 5.
|
||||
- **vif** is calculated by
|
||||
- First, regressing each predictor over all other predictors, and recording R-squared for each regression.
|
||||
- Second, computing vif as 1/(1- R_squared)
|
||||
|
||||
</help>
|
||||
</tool>
|
||||
Executable
+146
@@ -0,0 +1,146 @@
|
||||
#!/usr/bin/env python
|
||||
|
||||
from galaxy import eggs
|
||||
|
||||
import sys, string
|
||||
from rpy import *
|
||||
import numpy
|
||||
|
||||
#export PYTHONPATH=~/galaxy/lib/
|
||||
#running command python partialR_square.py reg_inp.tab 4 1,2,3 partialR_result.tabular
|
||||
|
||||
def stop_err(msg):
|
||||
sys.stderr.write(msg)
|
||||
sys.exit()
|
||||
|
||||
def sscombs(s):
|
||||
if len(s) == 1:
|
||||
return [s]
|
||||
else:
|
||||
ssc = sscombs(s[1:])
|
||||
return [s[0]] + [s[0]+comb for comb in ssc] + ssc
|
||||
|
||||
|
||||
infile = sys.argv[1]
|
||||
y_col = int(sys.argv[2])-1
|
||||
x_cols = sys.argv[3].split(',')
|
||||
outfile = sys.argv[4]
|
||||
|
||||
print "Predictor columns: %s; Response column: %d" %(x_cols,y_col+1)
|
||||
fout = open(outfile,'w')
|
||||
|
||||
for i, line in enumerate( file ( infile )):
|
||||
line = line.rstrip('\r\n')
|
||||
if len( line )>0 and not line.startswith( '#' ):
|
||||
elems = line.split( '\t' )
|
||||
break
|
||||
if i == 30:
|
||||
break # Hopefully we'll never get here...
|
||||
|
||||
if len( elems )<1:
|
||||
stop_err( "The data in your input dataset is either missing or not formatted properly." )
|
||||
|
||||
y_vals = []
|
||||
x_vals = []
|
||||
|
||||
for k,col in enumerate(x_cols):
|
||||
x_cols[k] = int(col)-1
|
||||
x_vals.append([])
|
||||
"""
|
||||
try:
|
||||
float( elems[x_cols[k]] )
|
||||
except:
|
||||
try:
|
||||
msg = "This operation cannot be performed on non-numeric column %d containing value '%s'." %( col, elems[x_cols[k]] )
|
||||
except:
|
||||
msg = "This operation cannot be performed on non-numeric data."
|
||||
stop_err( msg )
|
||||
"""
|
||||
NA = 'NA'
|
||||
for ind,line in enumerate( file( infile )):
|
||||
if line and not line.startswith( '#' ):
|
||||
try:
|
||||
fields = line.split("\t")
|
||||
try:
|
||||
yval = float(fields[y_col])
|
||||
except Exception, ey:
|
||||
yval = r('NA')
|
||||
#print >>sys.stderr, "ey = %s" %ey
|
||||
y_vals.append(yval)
|
||||
for k,col in enumerate(x_cols):
|
||||
try:
|
||||
xval = float(fields[col])
|
||||
except Exception, ex:
|
||||
xval = r('NA')
|
||||
#print >>sys.stderr, "ex = %s" %ex
|
||||
x_vals[k].append(xval)
|
||||
except:
|
||||
pass
|
||||
|
||||
x_vals1 = numpy.asarray(x_vals).transpose()
|
||||
dat= r.list(x=array(x_vals1), y=y_vals)
|
||||
|
||||
set_default_mode(NO_CONVERSION)
|
||||
try:
|
||||
full = r.lm(r("y ~ x"), data= r.na_exclude(dat)) #full model includes all the predictor variables specified by the user
|
||||
except RException, rex:
|
||||
stop_err("Error performing linear regression on the input data.\nEither the response column or one of the predictor columns contain no numeric values.")
|
||||
set_default_mode(BASIC_CONVERSION)
|
||||
|
||||
summary = r.summary(full)
|
||||
fullr2 = summary.get('r.squared','NA')
|
||||
|
||||
if fullr2 == 'NA':
|
||||
stop_error("Error in linear regression")
|
||||
|
||||
if len(x_vals) < 10:
|
||||
s = ""
|
||||
for ch in range(len(x_vals)):
|
||||
s += str(ch)
|
||||
else:
|
||||
stop_err("This tool only works with less than 10 predictors.")
|
||||
|
||||
print >>fout, "#Model\tR-sq\tpartial_R_Terms\tpartial_R_Value"
|
||||
all_combos = sorted(sscombs(s), key=len)
|
||||
all_combos.reverse()
|
||||
for j,cols in enumerate(all_combos):
|
||||
#if len(cols) == len(s): #Same as the full model above
|
||||
# continue
|
||||
if len(cols) == 1:
|
||||
x_vals1 = x_vals[int(cols)]
|
||||
else:
|
||||
x_v = []
|
||||
for col in cols:
|
||||
x_v.append(x_vals[int(col)])
|
||||
x_vals1 = numpy.asarray(x_v).transpose()
|
||||
dat= r.list(x=array(x_vals1), y=y_vals)
|
||||
set_default_mode(NO_CONVERSION)
|
||||
red = r.lm(r("y ~ x"), data= dat) #Reduced model
|
||||
set_default_mode(BASIC_CONVERSION)
|
||||
summary = r.summary(red)
|
||||
redr2 = summary.get('r.squared','NA')
|
||||
try:
|
||||
partial_R = (float(fullr2)-float(redr2))/(1-float(redr2))
|
||||
except:
|
||||
partial_R = 'NA'
|
||||
col_str = ""
|
||||
for col in cols:
|
||||
col_str = col_str + str(int(x_cols[int(col)]) + 1) + " "
|
||||
col_str.strip()
|
||||
partial_R_col_str = ""
|
||||
for col in s:
|
||||
if col not in cols:
|
||||
partial_R_col_str = partial_R_col_str + str(int(x_cols[int(col)]) + 1) + " "
|
||||
partial_R_col_str.strip()
|
||||
if len(cols) == len(s): #full model
|
||||
partial_R_col_str = "-"
|
||||
partial_R = "-"
|
||||
try:
|
||||
redr2 = "%.4f" %(float(redr2))
|
||||
except:
|
||||
pass
|
||||
try:
|
||||
partial_R = "%.4f" %(float(partial_R))
|
||||
except:
|
||||
pass
|
||||
print >>fout, "%s\t%s\t%s\t%s" %(col_str,redr2,partial_R_col_str,partial_R)
|
||||
Executable
+68
@@ -0,0 +1,68 @@
|
||||
<tool id="partialRsq" name="Compute partial R square" version="1.0.0">
|
||||
<description> </description>
|
||||
<command interpreter="python">
|
||||
partialR_square.py
|
||||
$input1
|
||||
$response_col
|
||||
$predictor_cols
|
||||
$out_file1
|
||||
1>/dev/null
|
||||
</command>
|
||||
<inputs>
|
||||
<param format="tabular" name="input1" type="data" label="Select data" help="Dataset missing? See TIP below."/>
|
||||
<param name="response_col" label="Response column (Y)" type="data_column" data_ref="input1" />
|
||||
<param name="predictor_cols" label="Predictor columns (X)" type="data_column" data_ref="input1" multiple="true">
|
||||
<validator type="no_options" message="Please select at least one column."/>
|
||||
</param>
|
||||
</inputs>
|
||||
<outputs>
|
||||
<data format="input" name="out_file1" metadata_source="input1" />
|
||||
</outputs>
|
||||
<requirements>
|
||||
<requirement type="python-module">rpy</requirement>
|
||||
</requirements>
|
||||
<tests>
|
||||
<!-- Test data with vlid values -->
|
||||
<test>
|
||||
<param name="input1" value="regr_inp.tabular"/>
|
||||
<param name="response_col" value="3"/>
|
||||
<param name="predictor_cols" value="1,2"/>
|
||||
<output name="out_file1" file="partialR_result.tabular"/>
|
||||
</test>
|
||||
|
||||
</tests>
|
||||
<help>
|
||||
|
||||
.. class:: infomark
|
||||
|
||||
**TIP:** If your data is not TAB delimited, use *Edit Datasets->Convert characters*
|
||||
|
||||
-----
|
||||
|
||||
.. class:: infomark
|
||||
|
||||
**What it does**
|
||||
|
||||
This tool computes the Partial R squared for all possible variable subsets using the following formula:
|
||||
|
||||
**Partial R squared = [SSE(without i: 1,2,...,p-1) - SSE (full: 1,2,..,i..,p-1) / SSE(without i: 1,2,...,p-1)]**, which denotes the case where the 'i'th predictor is dropped.
|
||||
|
||||
|
||||
|
||||
In general, **Partial R squared = [SSE(without i: 1,2,...,p-1) - SSE (full: 1,2,..,i..,p-1) / SSE(without i: 1,2,...,p-1)]**, where,
|
||||
|
||||
- SSE (full: 1,2,..,i..,p-1) = Sum of Squares left out by the full set of predictors SSE(X1, X2 … Xp)
|
||||
- SSE (full: 1,2,..,i..,p-1) = Sum of Squares left out by the set of predictors excluding; for example, if we omit the first predictor, it will be SSE(X2 … Xp).
|
||||
|
||||
|
||||
The 4 columns in the output are described below:
|
||||
|
||||
- Column 1 (Model): denotes the variables present in the model
|
||||
- Column 2 (R-sq): denotes the R-squared value corresponding to the model in Column 1
|
||||
- Column 3 (Partial R squared_Terms): denotes the variable/s for which Partial R squared is computed. These are the variables that are absent in the reduced model in Column 1. A '-' in this column indicates that the model in Column 1 is the Full model.
|
||||
- Column 4 (Partial R squared): denotes the Partial R squared value corresponding to the variable/s in Column 3. A '-' in this column indicates that the model in Column 1 is the Full model.
|
||||
|
||||
*R Development Core Team (2010). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL http://www.R-project.org.*
|
||||
|
||||
</help>
|
||||
</tool>
|
||||
Reference in New Issue
Block a user