rpy2 + 负二项式glm

Question

我正在尝试从rpy2中调用R的glm.nb函数：

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

测试：

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

返回一个错误：

    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

然而，当我运行简单的glm时，它可以正常运行。可能出现的问题是什么，如何进行调试？

- Dima Lituiev

Answer 1

这个问题涉及到R中矩阵和数组的数据结构。下面复制了一个在R中出现错误并进行修复的示例，以及在 rpy2 中复制修复过程的挑战和一个可行的解决方案：

R 错误和修复

library(MASS)

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

model2 <- glm.nb(y~x)

在 x[good, , drop = FALSE] * w 中存在非一致的数组

然而，有三种解决方法可用：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不能有效地将numpy矩阵传递到R矩阵中，因为在两个stat的简单glm()和MASS的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: 在 glm.fitter(x = X, y = Y, w = w, start = start, etastart = etastart) 中出现错误：未找到对象 'fit'。即使 numpy2ri.activate() 也无法转换 numpy 矩阵。

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

未实现错误：对象类型为'<class 'numpy.matrixlib.defmatrix.matrix'>'的转换“ri2py”未定义

---

解决方案

只需与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)) 

- Parfait

Answer 2

发生的情况是底层的 R 代码期望 "向量" 而不是数组，但 Python 对象是数组。

一个简单的解决方法是给你调用的 MASS 包中的 R 函数提供它想要/期望的东西。你测试中的以下行可以更改：

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

...to this:

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)))

- lgautier