rpy2 +负二项分布GLM

我试图从rpy2调用[R glm.nb：rpy2 +负二项分布GLM

from rpy2 import robjects 
from rpy2.robjects.packages import importr 

MASS = importr('MASS') 
stats = importr('stats') 

def glm_nb(x,y): 
    formula = robjects.Formula('y~x') 
    env = formula.environment 
    env["x"] = x 
    env["y"] = y 
    fitted = MASS.glm_nb(formula) 
#  fitted = stats.glm(formula) 
    return fitted

测试：当我运行简单

104   for k, v in kwargs.items(): 
    105    new_kwargs[k] = conversion.py2ri(v) 
--> 106   res = super(Function, self).__call__(*new_args, **new_kwargs) 
    107   res = conversion.ri2ro(res) 
    108   return res 

RRuntimeError: Error in x[good, , drop = FALSE] * w : non-conformable arrays

但是：

N = 100 
x = np.random.rand(N) 
y = x + np.random.poisson(10, N) 
fitted = glm_nb(x, np.round(y))

返回一个错误glm它运行正常。什么是问题，以及如何调试？

来源

2016-05-12 Dima Lituiev

的基本问题涉及在R.矩阵和阵列的数据结构下面再现你的错误中的R与修复，在rpy2复制修复程序的挑战，和工作溶液：

R错误和修复

library(MASS) 

# ARRAY 
x <- array(rnorm(100)) 
y <- as.integer(x) + array(rpois(100, 10)) 

model2 <- glm.nb(y~x)

Error in x[good, , drop = FALSE] * w : non-conformable arrays

然而，三个修复是可用的：1）使用矩阵（二维特殊类型的阵列）的; 2）等同定义的数组（指定dim参数）;和3）矩阵转换。请注意：重复限制的警告可能会出现取决于随机值，但仍会运行。

# MATRIX 
x <- matrix(rnorm(100)) 
y <- as.integer(x) + matrix(rpois(100, 10)) 

model1 <- glm.nb(y~x) 

# EQUIVALENT ARRAY 
x <- array(rnorm(100),c(100,1)) 
y <- as.integer(x) + matrix(rpois(100, 10),c(100,1)) 

model2 <- glm.nb(y~x) 

# EXPLICIT MATRIX CONVERSION (USED IN WORKING SOLUTION) 
x <- as.matrix(array(rnorm(100))) 
y <- as.integer(x) + as.matrix(array(rpois(100, 10))) 

model3 <- glm.nb(y~x)

挑战

Python的rpy2作为不同的错误都统计的简单glm()和大规模出现不能有效地从我的剧本的运作传递一个numpy的矩阵为R矩阵'glm.nb()：

import numpy as np 
from rpy2 import robjects 
from rpy2.robjects.packages import importr 
from rpy2.robjects.numpy2ri import numpy2ri 
MASS = importr('MASS') 

#rpy2 + negative binomial glm 
stats = importr('stats') 

def glm_nb(x,y): 
    formula = robjects.Formula('y~x') 
    env = formula.environment 
    env["x"] = x 
    env["y"] = y  
    fitted = MASS.glm_nb(formula) 
# fitted = stats.glm(formula) 
    return fitted 

N = 100 
x = np.random.rand(N) 
x = np.asmatrix(x)       # PYTHON CONVERSION TO MATRIX 
r_x = numpy2ri(x) 

# REPLACED NP.ROUND FOR AS.TYPE() TO COMPARE WITH R 
y = x.astype(int) + np.random.poisson(10, N) 
y = np.asmatrix(y)       # PYTHON CONVERSION TO MATRIX 
r_y = numpy2ri(y) 

fitted = glm_nb(r_x, r_y)

rpy2.rinterface.RRuntimeError: Error in glm.fitter(x = X, y = Y, w = w, start = start, etastart = etastart, : object 'fit' not found

即使numpy2ri.activate()未能将numpy的矩阵转换：

from rpy2.robjects import numpy2ri 
robjects.numpy2ri.activate() 
r_x = numpy2ri.ri2py(x) 
r_y = numpy2ri.ri2py(y)

NotImplementedError: Conversion 'ri2py' not defined for objects of type '<class 'numpy.matrixlib.defmatrix.matrix'>'

工作溶液

简单地与robjects.r()接口和具有R阵列对象转换为矩阵的工作。回想一下上面的第三个修复：

N = 100 
x = np.random.rand(N) 
r_x = numpy2ri(x) 

y = x.astype(int) + np.random.poisson(10, N) 
r_y = numpy2ri(y) 

from rpy2.robjects import r 
r.assign("y", r_y) 
r.assign("x", r_x) 
r("x <- as.matrix(x)") 
r("y <- as.matrix(y)") 
r("res <- glm.nb(y~x)") 

r_result = r("res[1:5]") 

# CONVERSION INTO PY DICTIONARY  
from rpy2.robjects import pandas2ri 
pandas2ri.activate() 
pyresult = pandas2ri.ri2py(r_result) 
print(pyresult)      # OUTPUTS COEFF, RESID, FITTED VALS, EFFECTS, R 

# OR OLDER DEPRECATED CONVERSION 
import pandas.rpy.common as com 
pyresult = com.convert_robj(r_result) 
print(pyresult)      # OUTPUTS COEFF, RESID, FITTED VALS, EFFECTS, R

命令行的解决方案

如果允许在你的应用程序，只需调用在Python将R模型脚本作为命令行子，绕过任何需要的rpy2，甚至通过根据需要提供参数：

from subprocess import Popen, PIPE 

command = 'Rscript.exe' 
path2Script = 'path/to/Script.R'  
args = ['arg1', 'arg2', 'arg3'] 

cmd = [command, path2Script] + args 

p = Popen(cmd,stdin= PIPE, stdout= PIPE, stderr= PIPE)    
output,error = p.communicate() 

if p.returncode == 0:    
    print('R OUTPUT:\n {0}'.format(output))    
else:     
    print('R ERROR:\n {0}'.format(error))

来源

2016-05-12 18:56:39 Parfait

发生的事情是底层R代码期望的是“矢量”而不是数组，但Python对象是数组。

一个简单的解决方法是给包MASS上的R函数调用它想要的/期望的。在测试下面的线是可以改变的：

fitted = glm_nb(x, np.round(y))

...这样的：

import array 
fitted = glm_nb(array.array('f', x), array.array('f', np.round(y)))

...或者这样：

from rpy2.robjects.vectors import FloatVector 
fitted = glm_nb(FloatVector(x), FloatVector(np.round(y)))

来源

2016-05-13 23:54:08 lgautier

rpy2 +负二项分布GLM

回答

相关问题