Monday, February 14, 2011

Correlation Test using R software

There is no reason why anyone interested in statistical computing will not choose R software, unless they are already wedded to existing software from their school days, or the choice was already made for them by their companies or their professors!. R boasts the most number of libraries and is open source and best of all FREE!
It takes time to master the syntax and we suggest that readers should try to download and test the R software in their machine if they had not done so.

In this post, we shall describe our online solver for testing correlation between two variables X and Y.
The solver is mainly a web interface to the cor.test function available in R. This function has the following spec for input arguments:


cor.test(x, y,
alternative = c("two.sided", "less", "greater"),
method = c("pearson", "kendall", "spearman"),
exact = NULL, conf.level = 0.95, continuity = FALSE, ...)


Here is how our solver for correlation test looks like at the moment.



When the user clicks on the submit button, the solver will call the R software, grabs its output, and print the results:

0001 > x <- c(44.4,45.9,41.9,53.3,44.7,44.1,50.7,45.2,60.1)

0002 > y <- c(2.6,3.1,2.5,5.0,3.6,4.0,5.2,2.8,3.8)

0003 >

0004 > cor.test(x,y, alternative="two.sided", method="pearson",conf.level=0.9,exact=TRUE, continuity=TRUE)

0005 

0006 Pearson's product-moment correlation

0007 

0008 data:  x and y

0009 t = 1.8411, df = 7, p-value = 0.1082

0010 alternative hypothesis: true correlation is not equal to 0

0011 90 percent confidence interval:

0012 -0.02223023  0.86697863

0013 sample estimates:

0014 cor

0015 0.5711816

0016 

0017 >

The printed p-value .1082 tells us not to reject the Null hypothesis that the correlation between
x and y is zero.

We will a scattergraph to the solver in the future.


Here is our interface code for the correlation test.

__version__ = "0.0.2"  
__catalog__ = "TEST-STAT-0036"
__date__    = "2011.02.14"
__author__  = "E.P. Adorio"


from   qp.fill.form     import Form, SingleSelectWidget, RadiobuttonsWidget, TextWidget, StringWidget, CheckboxWidget
from   qp.sites.extreme.lib.uicommon  import renderheader, renderfooter, processheader, processfooter
from   math             import sqrt
from   qp.sites.extreme.lib.webutils import vecRead
from   qp.sites.extreme.lib   import webutils
from   qp.fill.directory import Directory
from   qp.pub.common     import header, footer, page, redirect
from   qp.fill.css  import BASIC_FORM_CSS
from   qp.sites.extreme.lib.webutils import vecRead, matRead, matReadByColumn, runRcode
from   qp.sites.extreme.lib.qpyutils        import printRlines, showLogo


from   scipy          import stats

import time
from   scipy          import stats

import copy
import math


def getGrFile(grType):
    grfile   = GraphicsFile("%s" % grType)
    barefile = grfile.split(str("/"))[-1]
    return grfile, barefile

def getformDict(form):
    D={}
    for field in form.get_all_widgets():
        D[field.name] = form.get(field.name)
    return D    

def asColumns(S, ncolumns):
    """
    Args:
      ncolumns is number of columns.
      S is a string containing a space delimited columns of values.
    Return value
      a list consisting of the columns of S.
    """
    columns = [[] for i in range(ncolumns)]
    for i, v in enumerate(S.split()):
        columns[i % ncolumns].append( v)
    return columns

def asRvector(x):
    """
    Arg
      x - Python string array.
    Return value
      an R vector as a string c(v1,....vn)
    """
    S = "c(" 
    for v in x:
        S += str(v) + ","
    S += ")"
    return S.replace(",)", ")")

 


def SolveProblem(form):
    D    = getformDict(form)
    x, y = asColumns(D["pairedSample"], 2)
    alpha       = float(D["alpha"])
    method      = D["method"]
    alternative = D["alternative"]
    cont        = D["continuity_correction"]
    exact       = D["exact"]

    exact = "TRUE" if exact else "FALSE"
    cont  = "TRUE" if cont else  "FALSE"

    xstr = "x <- " + asRvector(x)+"\n"
    ystr = "y <- " + asRvector(y)+"\n"
    Rcode =xstr + ystr

    conflevel = 1-alpha

    Rcode += """
cor.test(x,y, alternative="%s", method="%s",conf.level=%s,exact=%s, continuity=%s)               
""" %(alternative, method,conflevel, exact, cont)
     
    # create the R command.
    Rcode = Rcode.replace(",)", ")")
    return runRcode(Rcode)


class CorrTestPage(Directory):
    def get_exports(self):
        yield ('',       'index',   'Correlation test for Paired Data', "R stat software corr.test for Correlation")

    def index[html](self):
        form = Form(enctype = "multipart/form-data")

        form.add(CheckboxWidget, name = "continuity_correction", \
                 value = True,
                )
        form.add(CheckboxWidget, name = "exact", \
                 value = True,
                )

        form.add(TextWidget,  name = "pairedSample",  \
                 rows = "10", cols = "30", value=""" 
          44.4 2.6
          45.9 3.1
          41.9 2.5
          53.3 5.0
          44.7 3.6
          44.1 4.0
          50.7 5.2
          45.2 2.8
          60.1 3.8
                """
                )

        form.add(StringWidget, name = "alpha", size = 5, value = "0.10")
        form.add(SingleSelectWidget, name="alternative", options=[
                 ("two.sided", "two.sided"),
                 ("less", "less"),
                 ("greater", "greater")])

        form.add(SingleSelectWidget, name="method", options=[
                 ("pearson", "Pearson"),
                 ("kendall", "Kendall"),
                 ("spearman", "Spearman")])


        form.add_hidden("time",   value = time.time())
        form.add_submit("submit", "submit")


        def render [html] ():
            renderheader("Correlation test using R-software")

            pairedSample = form.get_widget("pairedSample")
            alpha = form.get_widget("alpha")
            exact = form.get_widget("exact")
            method   = form.get_widget("method")
            alternative = form.get_widget("alternative")
            cont = form.get_widget("continuity_correction")
            """
            







MethodAlternativeAlphaContinuity CorrectionExact p value?
%s %s%s%s%s
X Y data
%s
""" %( method.render(), alternative.render(), alpha.render(),cont.render(), exact.render(), pairedSample.render())

renderfooter(form, __version__, __catalog__, __author__)

if not form.is_submitted():
return page('CorrTestPage', render(), style= BASIC_FORM_CSS)

def process [html] ():
processheader("Correlation test using R software corr.test() function")
calctime_start = time.time()

if True:
flag, output = SolveProblem(form)
printRlines(output)

if 0:
for widget in form.get_all_widgets():
"
%s %s
" % (widget.name, widget.value)
processfooter(form, calctime_start, "./", "./stest_t1samp")
process()

Tuesday, February 1, 2011

Creating a Unit Conversion Page

We had a unit conversion page before but I felt the initial code was too complicated to maintain. A rewrite was necessary and though the page is still at version 0.0.1, we feel that the rewrite is justified. Here is the code, using Python and QP. I need a web artist or graphic designer to add color and spice to the unit conversion page.

from qp.sites.extreme.lib.config import image_dir, homeurl



__version__ = "0.0.1 2011.01.28"
__date__    = "2011.01.28"
__catalog__ = "UNITS-CONV-001"
__author__  = "ernesto.adorio@gmail.com"
__url__     = homeurl + "solvers/units/"
__title__   = "Fundamental Units Conversion Page"


from qp.pub.common      import get_session
from qp.fill.css        import BASIC_FORM_CSS
from qp.fill.directory  import Directory
from quixote.errors     import PublishError
from quixote.util       import dump_request
from qp.fill.form       import Form
from qp.fill.widget    import CheckboxWidget,StringWidget,SingleSelectWidget, StringWidget, TextWidget
from qp.pub.common      import header, footer, page, redirect
from   qp.sites.extreme.lib.uicommon import renderheader, renderfooter, processheader, processfooter
import time
import urllib

from unitslib import Length, Mass, Time, convert


import urllib, urllib2, os

from  qp.sites.extreme.lib.webutils  import before, after
from  qp.sites.extreme.lib import config


def getstatusoutput(command):
     """
     Fredrick Lundh 
     http://mail.python.org/pipermail/python-list/2006-October/406444.html 
     """
     from subprocess import Popen, PIPE, STDOUT
     p = Popen(command, stdout=PIPE, stderr=STDOUT, shell=True)
     s = p.stdout.read().split("\n")
     return p.wait(), s


class UnitsDirectory(Directory):
    def get_exports(self):
        yield '', 'index', 'NetTools server', None



    def index (self):
        form  = Form(enctype = "multipart/form-data")  # enctype for file upload

 form.add(StringWidget, name = "qty_length", value = "1.0", size = 30)
        length_options = []
        for key in Length:
            length_options.append(key)
      
 form.add(SingleSelectWidget, name = "unit_length",
           value = "meter",  
           options=length_options
        )   


 form.add(StringWidget, name = "qty_mass", value = "1.0", size = 30)
        mass_options = []
        for key in Mass:
            mass_options.append(key)
      
 form.add(SingleSelectWidget, name = "unit_mass",
           value = "kilogram",  
           options=mass_options
        )   

 form.add(StringWidget, name = "qty_time", value = "1.0", size = 30)
        time_options = []
        for key in Time:
            time_options.append(key)
      
 form.add(SingleSelectWidget, name = "unit_time",
           value  = "second",  
           options=time_options
        )   

 form.add_hidden("time",   value = time.time()) 
        form.add_submit("submit", "submit")

 def render [html] ():
            renderheader(__title__)
            if 1:
              """
              

Length Unit

%s%s
""" % (form.get_widget("qty_length").render(), form.get_widget("unit_length").render()) """

Mass Unit

%s%s
""" % (form.get_widget("qty_mass").render(), form.get_widget("unit_mass").render()) """

Time Unit

%s%s
""" % (form.get_widget("qty_time").render(), form.get_widget("unit_time").render()) "The quantified unit will be converted in other equivalent units." renderfooter(form, __version__, __catalog__, __author__) if not form.is_submitted(): return page('UnitsDirectory', render(), style= BASIC_FORM_CSS) def process [html] (): processheader("Unit Utility") calctime_start = time.time() "" "" "" "" "
"
qty = float(form.get("qty_length"))
unit = form.get("unit_length")

"

Length Conversion

"


'' % (qty, unit) for u in Length.keys(): " " % (convert(Length, qty, unit, u),u) "
Input [%s %s]
%s%s
"
"
"
qty = float(form.get("qty_mass"))
unit = form.get("unit_mass")

"

Mass Conversion

"
'' % (qty, unit) for u in Mass.keys(): " " % (convert(Mass, qty, unit, u),u) "
Input [%s %s]
%s%s
"
"
"
qty = float(form.get("qty_time"))
unit = form.get("unit_time")
"

Time Conversion

"

'' % (qty, unit) for u in Time.keys(): " " % (convert(Time, qty, unit, u),u) "
Input [%s %s]
%s%s
"
"
" processfooter(form, calctime_start, homeurl, __url__) return process()

Yuck! the html output code is being translated by a browser! Use your browswerĊ› View Source facility to recover the original contents.

Codes are provided for educational use. Please make the proper attribution.


Only the fundamental units for Length, Mass and Time are considered here, we will add unit conversions for other quantities such as Area, Volume/Capacity, Work/Energy, Pressure and others.

The solver page may be accessed at http://extreme.adorio-research.org/solvers/units/.

The Python code units.py which is called by the above QP interface code above is at http://adorio-research.org/wordpress/?p=10377.

Our reference to conversion values is http://physics.nist.gov/Pubs/SP811/appenB9.html#LENGTH
for length units. The Mass and Time units are similarly obtained from NIST.