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

No estoy seguro de por qué estoy obteniendo 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 ipython aquí , que puede reproducir la inconsistencia.

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) 

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.