比萨铁塔——统计显著性检验

Dataset

比萨斜塔是意大利最大的旅游景点之一。几百年来这座塔慢慢靠向一边,最终达到5.5度的倾斜角度，在顶端水平偏离了近3米。年度数据pisa.csv文件记录了从1975年到1987年测量塔的倾斜，其中lean代表了偏离的角度。在这个任务,我们将尝试使用线性回归来估计倾斜率以及解释其系数和统计数据。

# 读取数据
import pandas
import matplotlib.pyplot as plt
pisa = pandas.DataFrame({"year": range(1975, 1988), 
                         "lean": [2.9642, 2.9644, 2.9656, 2.9667, 2.9673, 2.9688, 2.9696, 
                                  2.9698, 2.9713, 2.9717, 2.9725, 2.9742, 2.9757]})

print(pisa)
‘‘‘
      lean  year
0   2.9642  1975
1   2.9644  1976
2   2.9656  1977
3   2.9667  1978
4   2.9673  1979
5   2.9688  1980
6   2.9696  1981
7   2.9698  1982
8   2.9713  1983
9   2.9717  1984
10  2.9725  1985
11  2.9742  1986
12  2.9757  1987
‘‘‘
plt.scatter(pisa["year"], pisa["lean"])

Fit The Linear Model

从散点图中我们可以看到用曲线可以很好的拟合该数据。在之前我们利用线性回归来分析葡萄酒的质量以及股票市场，但在这个任务中，我们将学习如何理解关键的统计学概念。Statsmodels是Python中进行严格统计分析的一个库，对于线性模型，Statsmodels提供了足够多的统计方法以及适当的评估方法。sm.OLS这个类用于拟合线性模型，采取的优化方法是最小二乘法。OLS()不会自动添加一个截距到模型中，但是可以自己添加一列属性，使其值都是1即可产生截距。

import statsmodels.api as sm

y = pisa.lean # target
X = pisa.year  # features
X = sm.add_constant(X)  # add a column of 1‘s as the constant term

# OLS -- Ordinary Least Squares Fit
linear = sm.OLS(y, X)
# fit model
linearfit = linear.fit()
print(linearfit.summary())
‘‘‘
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                   lean   R-squared:                       0.988
Model:                            OLS   Adj. R-squared:                  0.987
Method:                 Least Squares   F-statistic:                     904.1
Date:                Mon, 25 Apr 2016   Prob (F-statistic):           6.50e-12
Time:                        13:30:20   Log-Likelihood:                 83.777
No. Observations:                  13   AIC:                            -163.6
Df Residuals:                      11   BIC:                            -162.4
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
const          1.1233      0.061     18.297      0.000         0.988     1.258
year           0.0009    3.1e-05     30.069      0.000         0.001     0.001
==============================================================================
Omnibus:                        0.310   Durbin-Watson:                   1.642
Prob(Omnibus):                  0.856   Jarque-Bera (JB):                0.450
Skew:                           0.094   Prob(JB):                        0.799
Kurtosis:                       2.108   Cond. No.                     1.05e+06
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.05e+06. This might indicate that there are
strong multicollinearity or other numerical problems.
‘‘‘

Define A Basic Linear Model

• 打印summary时发现有很多关于模型的信息，为了弄清楚这些统计指标，我们需要仔细研究一下标准的线性回归模型，下面模型中的ei是预测值和真实值的差，意味着我们默认模型的误差是正态分布的，均值为0.：
  

• 计算模型的残差（residuals）：

# Our predicted values of y
yhat = linearfit.predict(X)
print(yhat)
‘‘‘
[ 2.96377802  2.96470989  2.96564176  2.96657363  2.96750549  2.96843736
  2.96936923  2.9703011   2.97123297  2.97216484  2.9730967   2.97402857
  2.97496044]
‘‘‘
residuals = yhat - y
‘‘‘
residuals :Series (<class ‘pandas.core.series.Series‘>)
0    -0.000422
1     0.000310
2     0.000042
3    -0.000126
4     0.000205
5    -0.000363
6    -0.000231
7     0.000501
8    -0.000067
9     0.000465
10    0.000597
11   -0.000171
12   -0.000740
Name: lean, dtype: float64
‘‘‘

Histogram Of Residuals

• 之前我们用过直方图（histograms ）来显示数据的分布，现在我们也可以显示残差的分布，来确认它是否满足正态分布（其实有很多统计测试来检验正态分布）：

plt.hist(residuals, bins=5)

• 由于我们的数据集只有13个样本，因此这样画出来的直方图并没有太大意义，尽管中间最高的有4个样本

Sum Of Squares

许多线性回归模型的统计测量都依赖于三个平方测量值：Error (SSE), Regression Sum of Squares (RSS)以及Total Sum of Squares (TSS).

• Error (SSE)：真实值与预测值的差的平方和

  

• Regression Sum of Squares (RSS) :预测值和真实值的均值的差的平方和，其中的均值是真实值的均值。如果将预测值都设置为观测值的均值，RSS会非常低，但这并没有什么意义。反而是一个大的RSS和一个小的SSE表示一个很好的模型。

  

• Total Sum of Squares (TSS)：观测值与观测值的均值的差的平方和，大概就是训练集的方差。

  

  TSS=RSS+SSE：数据总量的方差 = 模型的方差+残差的方差

import numpy as np

# sum the (predicted - observed) squared
SSE = np.sum((yhat-y.values)**2)
‘‘‘
1.9228571428562889e-06
‘‘‘
# Average y
ybar = np.mean(y.values)

# sum the (mean - predicted) squared
RSS = np.sum((ybar-yhat)**2)
‘‘‘
0.00015804483516480448
‘‘‘

# sum the (mean - observed) squared
TSS = np.sum((ybar-y.values)**2)
‘‘‘
0.00015996769230769499
‘‘‘
print(TSS-RSS-SSE)
‘‘‘
3.42158959043e-17
‘‘‘

R-Squared

• 线性判定（coefficient of determination）也叫R-Squared，是用来测定线性依赖性的。它是一个数字用来告诉我们数据的总方差中模型的方差的占比：

• 前面提到一个低的SSE和一个高的RSS表示一个很好的模型拟合，这个R-Squared就表示了这个意思，介于0到1之间。

R2 = RSS/TSS
print(R2)
‘‘‘
0.987979715684
‘‘‘

T-Distribution

统计测验表明塔的倾斜程度与年份有关系，一个常见的统计显著性测试是student t-test。这个测试的基础是T分布。和正态分布很相似，都是钟型但是峰值较低。T检验是用于小样本（样本容量小于30）的两个平均值差异程度的检验方法。它是用T分布理论来推断差异发生的概率，从而判定两个平均数的差异是否显著。

from scipy.stats import t

# 100 values between -3 and 3
x = np.linspace(-3,3,100)

# Compute the pdf with 3 degrees of freedom
print(t.pdf(x=x, df=3))
‘‘‘
[ 0.02297204  0.02441481  0.02596406  0.02762847  0.0294174   0.031341
  0.03341025  0.03563701  0.03803403  0.04061509  0.04339497  0.04638952
  0.04961567  0.05309149  0.05683617  0.06086996  0.0652142   0.06989116
  0.07492395  0.08033633  0.08615245  0.09239652  0.0990924   0.10626304
  0.11392986  0.12211193  0.13082504  0.14008063  0.14988449  0.16023537
  0.17112343  0.18252859  0.1944188   0.20674834  0.21945618  0.23246464
  0.2456783   0.2589835   0.27224841  0.28532401  0.29804594  0.31023748
  0.32171351  0.33228555  0.34176766  0.34998293  0.35677032  0.36199128
  0.36553585  0.36732769  0.36732769  0.36553585  0.36199128  0.35677032
  0.34998293  0.34176766  0.33228555  0.32171351  0.31023748  0.29804594
  0.28532401  0.27224841  0.2589835   0.2456783   0.23246464  0.21945618
  0.20674834  0.1944188   0.18252859  0.17112343  0.16023537  0.14988449
  0.14008063  0.13082504  0.12211193  0.11392986  0.10626304  0.0990924
  0.09239652  0.08615245  0.08033633  0.07492395  0.06989116  0.0652142
  0.06086996  0.05683617  0.05309149  0.04961567  0.04638952  0.04339497
  0.04061509  0.03803403  0.03563701  0.03341025  0.031341    0.0294174
  0.02762847  0.02596406  0.02441481  0.02297204]
‘‘‘