В коде основного класса OLS в Python Pandas , я ищу помощь, чтобы уточнить, какие соглашения используются для стандартной ошибки и t-статистики, сообщаемой при выполнении взвешенного OLS.
Вот мой пример набора данных с некоторым импортом для использования Pandas и прямого использования scikits.statsmodels WLS:
import pandas
import numpy as np
from statsmodels.regression.linear_model import WLS
# Make some random data.
np.random.seed(42)
df = pd.DataFrame(np.random.randn(10, 3), columns=['a', 'b', 'weights'])
# Add an intercept term for direct use in WLS
df['intercept'] = 1
# Add a number (I picked 10) to stabilize the weight proportions a little.
df['weights'] = df.weights + 10
# Fit the regression models.
pd_wls = pandas.ols(y=df.a, x=df.b, weights=df.weights)
sm_wls = WLS(df.a, df[['intercept','b']], weights=df.weights).fit()
Я использую %cpaste для выполнения этого в IPython, а затем печатаю сводки обеих регрессий:
In [226]: %cpaste
Pasting code; enter '--' alone on the line to stop or use Ctrl-D.
:import pandas
:import numpy as np
:from statsmodels.regression.linear_model import WLS
:
:# Make some random data.
np:np.random.seed(42)
:df = pd.DataFrame(np.random.randn(10, 3), columns=['a', 'b', 'weights'])
:
:# Add an intercept term for direct use in WLS
:df['intercept'] = 1
:
:# Add a number (I picked 10) to stabilize the weight proportions a little.
:df['weights'] = df.weights + 10
:
:# Fit the regression models.
:pd_wls = pandas.ols(y=df.a, x=df.b, weights=df.weights)
:sm_wls = WLS(df.a, df[['intercept','b']], weights=df.weights).fit()
:--
In [227]: pd_wls
Out[227]:
-------------------------Summary of Regression Analysis-------------------------
Formula: Y ~ <x> + <intercept>
Number of Observations: 10
Number of Degrees of Freedom: 2
R-squared: 0.2685
Adj R-squared: 0.1770
Rmse: 2.4125
F-stat (1, 8): 2.9361, p-value: 0.1250
Degrees of Freedom: model 1, resid 8
-----------------------Summary of Estimated Coefficients------------------------
Variable Coef Std Err t-stat p-value CI 2.5% CI 97.5%
--------------------------------------------------------------------------------
x 0.5768 1.0191 0.57 0.5869 -1.4206 2.5742
intercept 0.5227 0.9079 0.58 0.5806 -1.2567 2.3021
---------------------------------End of Summary---------------------------------
In [228]: sm_wls.summ
sm_wls.summary sm_wls.summary_old
In [228]: sm_wls.summary()
Out[228]:
<class 'statsmodels.iolib.summary.Summary'>
"""
WLS Regression Results
==============================================================================
Dep. Variable: a R-squared: 0.268
Model: WLS Adj. R-squared: 0.177
Method: Least Squares F-statistic: 2.936
Date: Wed, 17 Jul 2013 Prob (F-statistic): 0.125
Time: 15:14:04 Log-Likelihood: -10.560
No. Observations: 10 AIC: 25.12
Df Residuals: 8 BIC: 25.72
Df Model: 1
==============================================================================
coef std err t P>|t| [95.0% Conf. Int.]
------------------------------------------------------------------------------
intercept 0.5227 0.295 1.770 0.115 -0.158 1.204
b 0.5768 0.333 1.730 0.122 -0.192 1.346
==============================================================================
Omnibus: 0.967 Durbin-Watson: 1.082
Prob(Omnibus): 0.617 Jarque-Bera (JB): 0.622
Skew: 0.003 Prob(JB): 0.733
Kurtosis: 1.778 Cond. No. 1.90
==============================================================================
"""
Обратите внимание на несоответствие стандартных ошибок: Pandas утверждает, что стандартные ошибки равны [0.9079, 1.0191], а statsmodels говорит [0.295, 0.333]..
Вернемся к коду, ссылку на который я дал вверху сообщения Я пытался отследить, откуда берется несоответствие.
Во-первых, вы можете видеть, что стандартные ошибки являются отчетами функции:
def _std_err_raw(self):
"""Returns the raw standard err values."""
return np.sqrt(np.diag(self._var_beta_raw))
Итак, глядя на self._var_beta_raw, я нахожу:
def _var_beta_raw(self):
"""
Returns the raw covariance of beta.
"""
x = self._x.values
y = self._y.values
xx = np.dot(x.T, x)
if self._nw_lags is None:
return math.inv(xx) * (self._rmse_raw ** 2)
else:
resid = y - np.dot(x, self._beta_raw)
m = (x.T * resid).T
xeps = math.newey_west(m, self._nw_lags, self._nobs, self._df_raw,
self._nw_overlap)
xx_inv = math.inv(xx)
return np.dot(xx_inv, np.dot(xeps, xx_inv))
В моем случае использования self._nw_lags всегда будет None, поэтому первая часть вызывает недоумение. Поскольку xx — это просто стандартный продукт матрицы регрессора: x.T.dot(x), мне интересно, как на это влияют веса. Термин self._rmse_raw происходит непосредственно из регрессии статистических моделей, подогнанной в конструкторе OLS, так что он определенно включает веса.
Это вызывает такие вопросы:
- Почему стандартная ошибка сообщается с весами, примененными в части RMSE, но не к переменным регрессора.
- Является ли это стандартной практикой, если вам нужны «непреобразованные» переменные (не могли бы вы также захотеть непреобразованный RMSE ??) Есть ли способ заставить Pandas вернуть полностью взвешенную версию стандартной ошибки?
- Почему все неправильное направление? В конструкторе вычисляется полная регрессия, подобранная статистическими моделями. Почему бы абсолютно каждой сводной статистике не браться прямо оттуда? Почему это смешивается и сопоставляется так, что некоторые из них получены из выходных данных statsmodels, а некоторые — из домашних вычислений Pandas?
Похоже, я могу согласовать вывод Pandas, выполнив следующие действия:
In [238]: xs = df[['intercept', 'b']]
In [239]: trans_xs = xs.values * np.sqrt(df.weights.values[:,None])
In [240]: trans_xs
Out[240]:
array([[ 3.26307961, -0.45116742],
[ 3.12503809, -0.73173821],
[ 3.08715494, 2.36918991],
[ 3.08776136, -1.43092325],
[ 2.87664425, -5.50382662],
[ 3.21158019, -3.25278836],
[ 3.38609639, -4.78219647],
[ 2.92835309, 0.19774643],
[ 2.97472796, 0.32996453],
[ 3.1158155 , -1.87147934]])
In [241]: np.sqrt(np.diag(np.linalg.inv(trans_xs.T.dot(trans_xs)) * (pd_wls._rmse_raw ** 2)))
Out[241]: array([ 0.29525952, 0.33344823])
Меня просто очень смущают эти отношения. Это то, что распространено среди статистиков: включение весов с частью RMSE, но затем выбор, взвешивать переменные или нет при расчете стандартной ошибки коэффициента? Если это так, то почему сами коэффициенты также не будут отличаться между пандами и статистическими моделями, поскольку они аналогичным образом получаются из переменных, сначала преобразованных статистическими моделями?
Для справки, вот полный набор данных, использованный в моем игрушечном примере (на случай, если np.random.seed недостаточно для воспроизводимости):
In [242]: df
Out[242]:
a b weights intercept
0 0.496714 -0.138264 10.647689 1
1 1.523030 -0.234153 9.765863 1
2 1.579213 0.767435 9.530526 1
3 0.542560 -0.463418 9.534270 1
4 0.241962 -1.913280 8.275082 1
5 -0.562288 -1.012831 10.314247 1
6 -0.908024 -1.412304 11.465649 1
7 -0.225776 0.067528 8.575252 1
8 -0.544383 0.110923 8.849006 1
9 0.375698 -0.600639 9.708306 1
