我的书签
添加书签
移除书签3.1.3.1 用“公式” 来在Python中指定统计模型
3.1.3.1.1 简单线性回归
给定两组观察值，x和y, 我们想要检验假设y是x的线性函数，换句话说:
$y = x * coef + intercept + e$
其中$e$是观察噪音。我们将使用statmodels module:
- 拟合一个线性模型。我们将使用简单的策略，普通最小二乘 (OLS)。
- 检验系数是否是非0。
首先，我们生成模型的虚拟数据:
In [9]:
import numpy as npx = np.linspace(-5, 5, 20)np.random.seed(1)# normal distributed noisey = -5 + 3*x + 4 * np.random.normal(size=x.shape)# Create a data frame containing all the relevant variablesdata = pandas.DataFrame({'x': x, 'y': y})
Python中的统计公式
然后我们指定一个OLS模型并且拟合它:
In [10]:
from statsmodels.formula.api import olsmodel = ols("y ~ x", data).fit()
我们可以检查fit产生的各种统计量:
In [26]:
print(model.summary())
OLS Regression Results==============================================================================Dep. Variable: y R-squared: 0.804Model: OLS Adj. R-squared: 0.794Method: Least Squares F-statistic: 74.03Date: Wed, 18 Nov 2015 Prob (F-statistic): 8.56e-08Time: 17:10:03 Log-Likelihood: -57.988No. Observations: 20 AIC: 120.0Df Residuals: 18 BIC: 122.0Df Model: 1Covariance Type: nonrobust==============================================================================coef std err t P>|t| [95.0% Conf. Int.]------------- ------------- ------------- ------------- ------------- -------------Intercept -5.5335 1.036 -5.342 0.000 -7.710 -3.357x 2.9369 0.341 8.604 0.000 2.220 3.654==============================================================================Omnibus: 0.100 Durbin-Watson: 2.956Prob(Omnibus): 0.951 Jarque-Bera (JB): 0.322Skew: -0.058 Prob(JB): 0.851Kurtosis: 2.390 Cond. No. 3.03==============================================================================Warnings:[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
术语:
Statsmodels使用统计术语: statsmodel的y变量被称为‘endogenous’而x变量被称为exogenous。更详细的讨论见这里。
为了简化，y (endogenous) 是你要尝试预测的值，而 x (exogenous) 代表用来进行这个预测的特征。
练习
从以上模型中取回估计参数。提示: 使用tab-完成来找到相关的属性。
3.1.3.1.2 类别变量: 比较组或多个类别
让我们回到大脑大小的数据:
In [27]:
data = pandas.read_csv('examples/brain_size.csv', sep=';', na_values=".")
我们可以写一个比较，用线性模型比较男女IQ:
In [28]:
model = ols("VIQ ~ Gender + 1", data).fit()print(model.summary())
OLS Regression Results==============================================================================Dep. Variable: VIQ R-squared: 0.015Model: OLS Adj. R-squared: -0.010Method: Least Squares F-statistic: 0.5969Date: Wed, 18 Nov 2015 Prob (F-statistic): 0.445Time: 17:34:10 Log-Likelihood: -182.42No. Observations: 40 AIC: 368.8Df Residuals: 38 BIC: 372.2Df Model: 1Covariance Type: nonrobust==================================================================================coef std err t P>|t| [95.0% Conf. Int.]------------- ------------- ------------- ------------- ------------- -----------------Intercept 109.4500 5.308 20.619 0.000 98.704 120.196Gender[T.Male] 5.8000 7.507 0.773 0.445 -9.397 20.997==============================================================================Omnibus: 26.188 Durbin-Watson: 1.709Prob(Omnibus): 0.000 Jarque-Bera (JB): 3.703Skew: 0.010 Prob(JB): 0.157Kurtosis: 1.510 Cond. No. 2.62==============================================================================Warnings:[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
特定模型的提示
强制类别: ‘Gender’ 被自动识别为类别变量，因此，它的每一个不同值都被处理为不同的实体。 使用:
In [29]:
model = ols('VIQ ~ C(Gender)', data).fit()
可以将一个整数列强制作为类别处理。
截距: 我们可以在公式中用-1删除截距，或者用+1强制使用截距。
默认，statsmodel将带有K和可能值的类别变量处理为K-1'虚拟变量' (最后一个水平被吸收到截距项中)。在绝大多数情况下，这都是很好的默认选择 - 但是，为类别变量指定不同的编码方式也是可以的 (http://statsmodels.sourceforge.net/devel/contrasts.html)。。)
FSIQ和PIQ差异的t-检验
要比较不同类型的IQ，我们需要创建一个"长形式"的表格，用一个类别变量来标识IQ类型:
In [30]:
data_fisq = pandas.DataFrame({'iq': data['FSIQ'], 'type': 'fsiq'})data_piq = pandas.DataFrame({'iq': data['PIQ'], 'type': 'piq'})data_long = pandas.concat((data_fisq, data_piq))print(data_long)
iq type0 133 fsiq1 140 fsiq2 139 fsiq3 133 fsiq4 137 fsiq5 99 fsiq6 138 fsiq7 92 fsiq8 89 fsiq9 133 fsiq10 132 fsiq11 141 fsiq12 135 fsiq13 140 fsiq14 96 fsiq15 83 fsiq16 132 fsiq17 100 fsiq18 101 fsiq19 80 fsiq20 83 fsiq21 97 fsiq22 135 fsiq23 139 fsiq24 91 fsiq25 141 fsiq26 85 fsiq27 103 fsiq28 77 fsiq29 130 fsiq.. ... ...10 124 piq11 128 piq12 124 piq13 147 piq14 90 piq15 96 piq16 120 piq17 102 piq18 84 piq19 86 piq20 86 piq21 84 piq22 134 piq23 128 piq24 102 piq25 131 piq26 84 piq27 110 piq28 72 piq29 124 piq30 132 piq31 137 piq32 110 piq33 86 piq34 81 piq35 128 piq36 124 piq37 94 piq38 74 piq39 89 piq[80 rows x 2 columns]
In [31]:
model = ols("iq ~ type", data_long).fit()print(model.summary())
OLS Regression Results==============================================================================Dep. Variable: iq R-squared: 0.003Model: OLS Adj. R-squared: -0.010Method: Least Squares F-statistic: 0.2168Date: Wed, 18 Nov 2015 Prob (F-statistic): 0.643Time: 18:16:40 Log-Likelihood: -364.35No. Observations: 80 AIC: 732.7Df Residuals: 78 BIC: 737.5Df Model: 1Covariance Type: nonrobust===============================================================================coef std err t P>|t| [95.0% Conf. Int.]------------- ------------- ------------- ------------- ------------- --------------Intercept 113.4500 3.683 30.807 0.000 106.119 120.781type[T.piq] -2.4250 5.208 -0.466 0.643 -12.793 7.943==============================================================================Omnibus: 164.598 Durbin-Watson: 1.531Prob(Omnibus): 0.000 Jarque-Bera (JB): 8.062Skew: -0.110 Prob(JB): 0.0178Kurtosis: 1.461 Cond. No. 2.62==============================================================================Warnings:[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
我们可以看到我们获得了与前面t-检验相同的值，以及相同的对应iq type的p-值:
In [32]:
stats.ttest_ind(data['FSIQ'], data['PIQ'])
Out[32]:
(0.46563759638096403, 0.64277250094148408)
