rpy2:在python中调用R函数的一个实例

python与R是当前数据科学计算的两大支柱,就我个人的使用经验而言,R更直观、简单和上手一些。很多专业的统计分析Python并没有提供R中的对应体,而你想要使用Python做数据分析,这时候就需要使用rpy2包来构建这个桥梁了。

比如我最近遇到一个分析问题,Python无法计算非参数统计检验wilcoxon test的置信区间,如果你仔细查看Python提供的非参数检验,你会发现它使用的是正态逼近,这在样本量大的时候(根据中心极限定理服从正态分布)当然可以使用Python计算,当如果你是小样本,比如大多数生物医学数据处理与分析中普遍样本少的可怜。在Stack overflow上有人讨论过并检查统计检验的p值,结论是算检验,R更靠谱些,优先采用。

言归正传,在银行统计工作室rpy2使用示例一文中对rpy2包各方面的使用都有介绍,加rpy2官方文档基本上可以掌握rpy2的使用,这里我提供这两天实现的一个实例——从python中调用R的wilcox.test函数进行非参数检验,如果大家有这方面需求可以作为一个参考。

代码已经封装为一个函数,略写了一下文档。

# 推荐使用conda管理环境
# conda create --name test python=3.6
# source activate test
# conda install rpy2  # should add conda-forge channel
# reference link: <http://www.cnblogs.com/cloudtj/articles/6372200.html>, <https://rpy2.readthedocs.io/en/version_2.8.x/index.html>

def pyWilcox(x, y=None, alternative='two.sided', mu=0, paired=False, exact=None, correct=True, conf_interval=True, conf_level=0.95):
    ''' 
    Run wilcoxon test using wilcox.test in R stats package, default is 'two sided' test.
        return p value, statistical value and confidence interval in a dictory. 

    Arguments:
        
        x   
            numeric vector of data values. Non-finite (e.g., infinite or missing) values will be omitted.
        y   
            an optional numeric vector of data values: as with x non-finite values will be omitted.
        alternative 
            a character string specifying the alternative hypothesis, must be one of "two.sided" (default), "greater" or "less". You can specify just the initial letter.
        mu  
            a number specifying an optional parameter used to form the null hypothesis. See ‘Details’.
        paired  
            a logical indicating whether you want a paired test.
        exact   
            a logical indicating whether an exact p-value should be computed.
        correct 
            a logical indicating whether to apply continuity correction in the normal approximation for the p-value.
        conf_interval   
            a logical indicating whether a confidence interval should be computed.
        conf_level  
            confidence level of the interval.

    Example:
        
        x = [i for i in range(0, 10)]
        y = [i for i in range(10, 20)]
        # two sample test
        pyWilcox(x, y)
        # one sample test
        pyWilcox(x)

        # output:
        #
        # {'p_value': [1.082508822446903e-05], 'statistic': [0.0], 'conf_interval': [-13.0, -7.0]}
        # {'p_value': [0.009151688852650072], 'statistic': [45.0], 'conf_interval': [2.500027475911944, 7.499972524088056]}
        

    Note:
        
        More information please run help('wilcox.test.default') in R console
    '''

    # 载入r对象
    from rpy2 import robjects
    # 载入导入包函数
    from rpy2.robjects.packages import importr
    # 将stats包导入为模块
    stats = importr('stats')
# When one wants to create a vector from Python, either the class Vector or the convenience classes IntVector, FloatVector, BoolVector, StrVector can be used.

    # 将列表x转换为r可识别数据对象
    x = robjects.FloatVector(x)
    
    # 将参数中的.替换为_,解决不兼容问题, 来自rpy2文档函数部分
    def iamfeelinglucky(func):
        def f(*args, **kwargs):
            d = {}
            for k, v in kwargs.items():
                d[k.replace('_', '.')] = v
            return func(**d)
        return f

    # 矫正参数名
    wilcox = iamfeelinglucky(stats.wilcox_test_default)

    # None类型似乎没有相应的函数,只能用条件语句进行判断
    if y != None:
        y = robjects.FloatVector(y)
        if exact != None:
            pr = wilcox(x = x, y = y, alternative = alternative, mu = mu, paired = paired, exact = exact, correct = correct, conf_int = conf_interval, conf_level = conf_level)
        else:
            pr = wilcox(x = x, y = y, alternative = alternative, mu = mu, paired = paired, conf_int = conf_interval, conf_level = conf_level)       
    else:
        if exact != None:
            pr = wilcox(x = x,  alternative = alternative, mu = mu, exact = exact, correct = correct, conf_int = conf_interval, conf_level = conf_level)
        else:
            pr = wilcox(x = x,  alternative = alternative, mu = mu, correct = correct, conf_int = conf_interval, conf_level = conf_level) 
           
    print(pr)
    res = list(pr)
    # 返回结果中需要的值构建字典
    res = {"p_value":list(res[2]), "statistic":list(res[0]), "conf_interval":list(res[7])}   
    return(res)

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

推荐阅读更多精彩内容

  • 她知道自己是有些洁癖的。二十来岁的时候,她的洁癖比现在还严重。那时她不能容忍男朋友用她的杯子喝水,在床上说的最多的...
    冬惊阅读 2,256评论 5 16
  • 喜欢每天早上骑车去公司的路上,天空晴朗☀️,微风的温度随着季节在变换。 不确定自己会在一个地方待多久,但喜欢这时的...
    simpleour阅读 197评论 0 0
  • 不曾想自己有一种爱好,那就是手机摄影。说不上来专业,但却时常充满感动。好的照片不单是为了取悦众人,更是能感动自己。...
    Jack叔叔阅读 758评论 0 5
  • javascript观察者模式学习笔记 原文地址 首先是观察者模式代码 在观察着模式的运用当中。react组件间的...
    strong9527阅读 384评论 0 0