Diferencia en los modelos de estadísticas de Python OLS y R''s lm

ols in python statsmodels

Question

No estoy seguro de por qué obtengo resultados ligeramente diferentes para un OLS simple, dependiendo de si utilizo la interfaz rpy experimental de panda para realizar la regresión en R o si uso statsmodels en Python.

import pandas
from rpy2.robjects import r

from functools import partial

loadcsv = partial(pandas.DataFrame.from_csv,
                  index_col="seqn", parse_dates=False)

demoq = loadcsv("csv/DEMO.csv")
rxq = loadcsv("csv/quest/RXQ_RX.csv")

num_rx = {}
for seqn, num in rxq.rxd295.iteritems():
    try:
        val = int(num)
    except ValueError:
        val = 0
    num_rx[seqn] = val

series = pandas.Series(num_rx, name="num_rx")
demoq = demoq.join(series)

import pandas.rpy.common as com
df = com.convert_to_r_dataframe(demoq)
r.assign("demoq", df)
r(''lmout <- lm(demoq$num_rx ~ demoq$ridageyr)'')  # run the regression
r(''print(summary(lmout))'')  # print from R

De R , obtengo el siguiente resumen:

Call:
lm(formula = demoq$num_rx ~ demoq$ridageyr)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.9086 -0.6908 -0.2940  0.1358 15.7003 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)    -0.1358216  0.0241399  -5.626 1.89e-08 ***
demoq$ridageyr  0.0358161  0.0006232  57.469  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 1.545 on 9963 degrees of freedom
Multiple R-squared: 0.249,  Adjusted R-squared: 0.2489 
F-statistic:  3303 on 1 and 9963 DF,  p-value: < 2.2e-16

Usando statsmodels.api para hacer el OLS:

import statsmodels.api as sm
results = sm.OLS(demoq.num_rx, demoq.ridageyr).fit()
results.summary()

Los resultados son similares a la salida de R pero no son los mismos:

OLS Regression Results
Adj. R-squared:  0.247
Log-Likelihood:  -18488.
No. Observations:    9965    AIC:   3.698e+04
Df Residuals:    9964    BIC:   3.698e+04
             coef   std err  t     P>|t|    [95.0% Conf. Int.]
ridageyr     0.0331  0.000   82.787    0.000        0.032 0.034

El proceso de instalación es un poco engorroso. Pero, hay un cuaderno de ipython here , que puede reproducir la inconsistencia.

Answer 1

Parece que Python no agrega una intercepción de forma predeterminada a tu expresión, mientras que R lo hace cuando usas la interfaz de fórmula ...

Esto significa que encajaste en dos modelos diferentes. Tratar

lm( y ~ x - 1, data)

en R para excluir la intercepción, o en su caso y con una notación algo más estándar

lm(num_rx ~ ridageyr - 1, data=demoq)

Answer 2

Tenga en cuenta que todavía puede usar statsmodels.formula.api desde statsmodels.formula.api :

from statsmodels.formula.api import ols

results = ols(''num_rx ~ ridageyr'', demoq).fit()
results.summary()

Creo que usa patsy en el backend para traducir la expresión de fórmula, y la intercepción se agrega automáticamente.