Obtener errores estándar de las regresiones usando RPY2
Pregunta
Estoy usando RPY2 para regresiones. El objeto devuelto tiene una lista que incluye coeficientes, residuos, valores ajustados, rango del modelo ajustado, etc.)
Sin embargo, no puedo encontrar los errores estándar (ni el R^2) en el objeto FIT. Ejecutando el modelo LM directamente en R, los errores estándar se muestran con el comando Resumen, pero no puedo acceder a ellos directamente en el marco de datos del modelo.
¿Cómo puedo obtener esta información con RPY2?
El código de pitón de muestra es
from scipy import random
from numpy import hstack, array, matrix
from rpy2 import robjects
from rpy2.robjects.packages import importr
def test_regress():
stats=importr('stats')
x=random.uniform(0,1,100).reshape([100,1])
y=1+x+random.uniform(0,1,100).reshape([100,1])
x_in_r=create_r_matrix(x, x.shape[1])
y_in_r=create_r_matrix(y, y.shape[1])
formula=robjects.Formula('y~x')
env = formula.environment
env['x']=x_in_r
env['y']=y_in_r
fit=stats.lm(formula)
coeffs=array(fit[0])
resids=array(fit[1])
fitted_vals=array(fit[4])
return(coeffs, resids, fitted_vals)
def create_r_matrix(py_array, ncols):
if type(py_array)==type(matrix([1])) or type(py_array)==type(array([1])):
py_array=py_array.tolist()
r_vector=robjects.FloatVector(flatten_list(py_array))
r_matrix=robjects.r['matrix'](r_vector, ncol=ncols)
return r_matrix
def flatten_list(source):
return([item for sublist in source for item in sublist])
test_regress()
Solución
Entonces esto parece funcionar para mí:
def test_regress():
stats=importr('stats')
x=random.uniform(0,1,100).reshape([100,1])
y=1+x+random.uniform(0,1,100).reshape([100,1])
x_in_r=create_r_matrix(x, x.shape[1])
y_in_r=create_r_matrix(y, y.shape[1])
formula=robjects.Formula('y~x')
env = formula.environment
env['x']=x_in_r
env['y']=y_in_r
fit=stats.lm(formula)
coeffs=array(fit[0])
resids=array(fit[1])
fitted_vals=array(fit[4])
modsum = base.summary(fit)
rsquared = array(modsum[7])
se = array(modsum.rx2('coefficients')[2:4])
return(coeffs, resids, fitted_vals, rsquared, se)
Aunque, como dije, esta es literalmente mi primera incursión en RPY2, por lo que puede haber una mejor manera de hacerlo. Pero esta versión parece generar matrices que contienen el valor R cuadrado junto con los errores estándar.
Puedes usar print(modsum.names)
para ver los nombres de los componentes del objeto R (algo así como names(modsum)
en r) y luego .rx
y .rx2
son el equivalente de [
y [[
en R.
Otros consejos
@Joran: Muy bien. Yo diría que es más o menos la forma de hacerlo.
from rpy2 import robjects
from rpy2.robjects.packages import importr
base = importr('base')
stats = importr('stats') # import only once !
def test_regress():
x = base.matrix(stats.runif(100), nrow = 100)
y = (x.ro + base.matrix(stats.runif(100), nrow = 100)).ro + 1 # not so nice
formula = robjects.Formula('y~x')
env = formula.environment
env['x'] = x
env['y'] = y
fit = stats.lm(formula)
coefs = stats.coef(fit)
resids = stats.residuals(fit)
fitted_vals = stats.fitted(fit)
modsum = base.summary(fit)
rsquared = modsum.rx2('r.squared')
se = modsum.rx2('coefficients')[2:4]
return (coefs, resids, fitted_vals, rsquared, se)