{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"import statsmodels.api as sm\n",
"from statsmodels.stats.outliers_influence import variance_inflation_factor\n",
"import scipy.stats.stats as scss\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline\n",
"\n",
"# 设置随机数种子，使运行结果可复现\n",
"np.random.seed(2046)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" x1 | \n",
" x2 | \n",
" x3 | \n",
" y | \n",
"
\n",
" \n",
" \n",
" \n",
" | 0 | \n",
" 0 | \n",
" 0 | \n",
" 0.894992 | \n",
" 0.297397 | \n",
"
\n",
" \n",
" | 1 | \n",
" 0 | \n",
" 1 | \n",
" 0.596078 | \n",
" -0.888134 | \n",
"
\n",
" \n",
" | 2 | \n",
" 0 | \n",
" 2 | \n",
" 0.435564 | \n",
" -1.979664 | \n",
"
\n",
" \n",
" | 3 | \n",
" 1 | \n",
" 0 | \n",
" 1.153073 | \n",
" 1.115046 | \n",
"
\n",
" \n",
" | 4 | \n",
" 1 | \n",
" 1 | \n",
" 1.340793 | \n",
" 0.064129 | \n",
"
\n",
" \n",
" | 5 | \n",
" 1 | \n",
" 2 | \n",
" 1.621597 | \n",
" -0.930433 | \n",
"
\n",
" \n",
" | 6 | \n",
" 2 | \n",
" 0 | \n",
" 2.139849 | \n",
" 2.060672 | \n",
"
\n",
" \n",
" | 7 | \n",
" 2 | \n",
" 1 | \n",
" 2.488654 | \n",
" 1.067725 | \n",
"
\n",
" \n",
" | 8 | \n",
" 2 | \n",
" 2 | \n",
" 2.049478 | \n",
" -0.115373 | \n",
"
\n",
" \n",
" | 9 | \n",
" 0 | \n",
" 0 | \n",
" 0.564523 | \n",
" 0.206383 | \n",
"
\n",
" \n",
" | 10 | \n",
" 0 | \n",
" 1 | \n",
" 0.270165 | \n",
" -1.018755 | \n",
"
\n",
" \n",
" | 11 | \n",
" 0 | \n",
" 2 | \n",
" 0.596075 | \n",
" -2.018875 | \n",
"
\n",
" \n",
" | 12 | \n",
" 1 | \n",
" 0 | \n",
" 1.291240 | \n",
" 1.134136 | \n",
"
\n",
" \n",
" | 13 | \n",
" 1 | \n",
" 1 | \n",
" 1.251871 | \n",
" 0.026758 | \n",
"
\n",
" \n",
" | 14 | \n",
" 1 | \n",
" 2 | \n",
" 1.110560 | \n",
" -1.144212 | \n",
"
\n",
" \n",
" | 15 | \n",
" 2 | \n",
" 0 | \n",
" 2.805694 | \n",
" 2.271965 | \n",
"
\n",
" \n",
" | 16 | \n",
" 2 | \n",
" 1 | \n",
" 2.465100 | \n",
" 1.058364 | \n",
"
\n",
" \n",
" | 17 | \n",
" 2 | \n",
" 2 | \n",
" 2.258865 | \n",
" -0.081940 | \n",
"
\n",
" \n",
"
\n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# 为在Matplotlib中显示中文，设置特殊字体\n",
"plt.rcParams[\"font.sans-serif\"] = [\"SimHei\"]\n",
"# 正确显示负号\n",
"plt.rcParams['axes.unicode_minus'] = False\n",
"plt.rcParams.update({'font.size': 20})\n",
"# 创建一个图形框\n",
"fig = plt.figure(figsize=(18, 6), dpi=100)\n",
"# 计算相关系数\n",
"corr = data.corr(method='pearson')\n",
"var = [('x1', 'x2'), ('x1', 'x3'), ('x2', 'x3')]\n",
"for i, (v1, v2) in enumerate(var):\n",
" ax = fig.add_subplot(1, 3, i + 1)\n",
" ax.set_xlabel(f'${v1}$')\n",
" ax.set_ylabel(f'${v2}$')\n",
" ax.scatter(data[v1], data[v2])\n",
" ax.annotate(f'相关系数：{corr[v1][v2]:.2f}', (0.5, 0.9),\n",
" xycoords='axes fraction', ha='center', va='center')\n",
" ax.set_xticks([])\n",
" ax.set_yticks([])\n",
"plt.savefig('multicollinearity_data.png', dpi=200)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def train_model(X, Y):\n",
" '''\n",
" 训练线性回归模型\n",
" '''\n",
" model = sm.OLS(Y, X)\n",
" re = model.fit()\n",
" return re"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: y R-squared: 0.431\n",
"Model: OLS Adj. R-squared: 0.396\n",
"Method: Least Squares F-statistic: 12.14\n",
"Date: Tue, 07 Nov 2023 Prob (F-statistic): 0.00307\n",
"Time: 17:46:01 Log-Likelihood: -23.867\n",
"No. Observations: 18 AIC: 51.73\n",
"Df Residuals: 16 BIC: 53.51\n",
"Df Model: 1 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"const -0.9094 0.360 -2.525 0.023 -1.673 -0.146\n",
"x1 0.9719 0.279 3.484 0.003 0.380 1.563\n",
"==============================================================================\n",
"Omnibus: 6.365 Durbin-Watson: 2.636\n",
"Prob(Omnibus): 0.041 Jarque-Bera (JB): 1.638\n",
"Skew: 0.020 Prob(JB): 0.441\n",
"Kurtosis: 1.523 Cond. No. 2.92\n",
"==============================================================================\n",
"\n",
"Notes:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/tgbaggio/opt/anaconda3/lib/python3.8/site-packages/scipy/stats/stats.py:1603: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=18\n",
" warnings.warn(\"kurtosistest only valid for n>=20 ... continuing \"\n"
]
}
],
"source": [
"# 用x1搭建回归模型\n",
"Y = data['y']\n",
"X = sm.add_constant(data['x1'])\n",
"re = train_model(X, Y)\n",
"print(re.summary())"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: y R-squared: 0.566\n",
"Model: OLS Adj. R-squared: 0.539\n",
"Method: Least Squares F-statistic: 20.84\n",
"Date: Tue, 07 Nov 2023 Prob (F-statistic): 0.000318\n",
"Time: 17:46:01 Log-Likelihood: -21.442\n",
"No. Observations: 18 AIC: 46.88\n",
"Df Residuals: 16 BIC: 48.66\n",
"Df Model: 1 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"const 1.1755 0.315 3.734 0.002 0.508 1.843\n",
"x2 -1.1130 0.244 -4.565 0.000 -1.630 -0.596\n",
"==============================================================================\n",
"Omnibus: 6.451 Durbin-Watson: 0.642\n",
"Prob(Omnibus): 0.040 Jarque-Bera (JB): 1.649\n",
"Skew: 0.035 Prob(JB): 0.438\n",
"Kurtosis: 1.519 Cond. No. 2.92\n",
"==============================================================================\n",
"\n",
"Notes:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/tgbaggio/opt/anaconda3/lib/python3.8/site-packages/scipy/stats/stats.py:1603: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=18\n",
" warnings.warn(\"kurtosistest only valid for n>=20 ... continuing \"\n"
]
}
],
"source": [
"# 用x2搭建回归模型\n",
"X1 = sm.add_constant(data['x2'])\n",
"re1 = train_model(X1, Y)\n",
"print(re1.summary())"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: y R-squared: 0.997\n",
"Model: OLS Adj. R-squared: 0.997\n",
"Method: Least Squares F-statistic: 2488.\n",
"Date: Tue, 07 Nov 2023 Prob (F-statistic): 1.21e-19\n",
"Time: 17:46:01 Log-Likelihood: 23.319\n",
"No. Observations: 18 AIC: -40.64\n",
"Df Residuals: 15 BIC: -37.97\n",
"Df Model: 2 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"const 0.2036 0.034 5.952 0.000 0.131 0.277\n",
"x1 0.9719 0.021 46.398 0.000 0.927 1.017\n",
"x2 -1.1130 0.021 -53.133 0.000 -1.158 -1.068\n",
"==============================================================================\n",
"Omnibus: 0.411 Durbin-Watson: 2.369\n",
"Prob(Omnibus): 0.814 Jarque-Bera (JB): 0.514\n",
"Skew: 0.278 Prob(JB): 0.774\n",
"Kurtosis: 2.387 Cond. No. 4.26\n",
"==============================================================================\n",
"\n",
"Notes:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/tgbaggio/opt/anaconda3/lib/python3.8/site-packages/scipy/stats/stats.py:1603: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=18\n",
" warnings.warn(\"kurtosistest only valid for n>=20 ... continuing \"\n"
]
}
],
"source": [
"# 用不相关的x1，x2搭建回归模型\n",
"X2 = sm.add_constant(data[['x1', 'x2']])\n",
"re2 = train_model(X2, Y)\n",
"print(re2.summary())"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/tgbaggio/opt/anaconda3/lib/python3.8/site-packages/scipy/stats/stats.py:1603: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=18\n",
" warnings.warn(\"kurtosistest only valid for n>=20 ... continuing \"\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: y R-squared: 0.431\n",
"Model: OLS Adj. R-squared: 0.396\n",
"Method: Least Squares F-statistic: 12.14\n",
"Date: Tue, 07 Nov 2023 Prob (F-statistic): 0.00307\n",
"Time: 17:46:01 Log-Likelihood: -23.867\n",
"No. Observations: 18 AIC: 51.73\n",
"Df Residuals: 16 BIC: 53.51\n",
"Df Model: 1 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"const -0.9094 0.360 -2.525 0.023 -1.673 -0.146\n",
"x1 0.9719 0.279 3.484 0.003 0.380 1.563\n",
"==============================================================================\n",
"Omnibus: 6.365 Durbin-Watson: 2.636\n",
"Prob(Omnibus): 0.041 Jarque-Bera (JB): 1.638\n",
"Skew: 0.020 Prob(JB): 0.441\n",
"Kurtosis: 1.523 Cond. No. 2.92\n",
"==============================================================================\n",
"\n",
"Notes:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
]
}
],
"source": [
"# 用强相关的x1，x3搭建模型\n",
"Y = data['y']\n",
"X = sm.add_constant(data['x1'])\n",
"re = train_model(X, Y)\n",
"print(re.summary())"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: y R-squared: 0.483\n",
"Model: OLS Adj. R-squared: 0.451\n",
"Method: Least Squares F-statistic: 14.96\n",
"Date: Tue, 07 Nov 2023 Prob (F-statistic): 0.00136\n",
"Time: 17:46:01 Log-Likelihood: -23.007\n",
"No. Observations: 18 AIC: 50.01\n",
"Df Residuals: 16 BIC: 51.79\n",
"Df Model: 1 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"const -1.4717 0.452 -3.254 0.005 -2.431 -0.513\n",
"x3 1.0901 0.282 3.867 0.001 0.493 1.688\n",
"==============================================================================\n",
"Omnibus: 4.131 Durbin-Watson: 2.709\n",
"Prob(Omnibus): 0.127 Jarque-Bera (JB): 1.375\n",
"Skew: 0.060 Prob(JB): 0.503\n",
"Kurtosis: 1.652 Cond. No. 4.41\n",
"==============================================================================\n",
"\n",
"Notes:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/tgbaggio/opt/anaconda3/lib/python3.8/site-packages/scipy/stats/stats.py:1603: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=18\n",
" warnings.warn(\"kurtosistest only valid for n>=20 ... continuing \"\n"
]
}
],
"source": [
"X1 = sm.add_constant(data['x3'])\n",
"re1 = train_model(X1, Y)\n",
"print(re1.summary())"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/tgbaggio/opt/anaconda3/lib/python3.8/site-packages/scipy/stats/stats.py:1603: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=18\n",
" warnings.warn(\"kurtosistest only valid for n>=20 ... continuing \"\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: y R-squared: 0.484\n",
"Model: OLS Adj. R-squared: 0.415\n",
"Method: Least Squares F-statistic: 7.040\n",
"Date: Tue, 07 Nov 2023 Prob (F-statistic): 0.00698\n",
"Time: 17:46:01 Log-Likelihood: -22.989\n",
"No. Observations: 18 AIC: 51.98\n",
"Df Residuals: 15 BIC: 54.65\n",
"Df Model: 2 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"const -1.5437 0.622 -2.480 0.025 -2.870 -0.217\n",
"x1 -0.1677 0.959 -0.175 0.864 -2.213 1.877\n",
"x3 1.2604 1.017 1.240 0.234 -0.907 3.428\n",
"==============================================================================\n",
"Omnibus: 3.209 Durbin-Watson: 2.726\n",
"Prob(Omnibus): 0.201 Jarque-Bera (JB): 1.239\n",
"Skew: 0.079 Prob(JB): 0.538\n",
"Kurtosis: 1.724 Cond. No. 14.5\n",
"==============================================================================\n",
"\n",
"Notes:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
]
}
],
"source": [
"X2 = sm.add_constant(data[['x1', 'x3']])\n",
"re2 = train_model(X2, Y)\n",
"print(re2.summary())"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"检测假设x1和x3同时不显著：\n",
"\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" VIF Factor | \n",
" features | \n",
"
\n",
" \n",
" \n",
" \n",
" | 0 | \n",
" 11.049590 | \n",
" const | \n",
"
\n",
" \n",
" | 1 | \n",
" 12.907932 | \n",
" x1 | \n",
"
\n",
" \n",
" | 2 | \n",
" 1.061102 | \n",
" x2 | \n",
"
\n",
" \n",
" | 3 | \n",
" 12.969034 | \n",
" x3 | \n",
"
\n",
" \n",
"
\n",
"