Python统计学：批判性实例教程与常见误用解构

Python统计学：批判性实例教程与常见误用解构

作为一名统计学教授，我发现当前许多Python统计学教程存在严重的简化和误导。它们往往忽略了统计学方法的适用条件，导致读者在实际应用中得出错误的结论。本文旨在通过具体的实例，揭示这些常见的误用，并提供更严谨的分析方法。

1. 线性回归的陷阱：数据分布与残差分析

问题背景： 假设我们想研究广告支出与销售额之间的关系。我们收集了一些数据，并使用线性回归模型进行分析。

数据描述： 数据集包含两个变量：广告支出（以千元为单位）和销售额（以万元为单位）。我们使用Python的pandas库加载数据，并使用statsmodels库进行线性回归分析。

import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
import numpy as np
from scipy import stats

# 生成一些示例数据，模拟非线性关系
np.random.seed(0)
x = np.linspace(1, 10, 100)
y = 2 * x + np.sin(x) * 5 + np.random.normal(0, 5, 100)

data = pd.DataFrame({'广告支出': x, '销售额': y})

# 添加截距项
X = sm.add_constant(data['广告支出'])

# 拟合线性回归模型
model = sm.OLS(data['销售额'], X)
results = model.fit()

# 打印回归结果
print(results.summary())

#  残差图
plt.figure(figsize=(10,5))
plt.scatter(results.fittedvalues, results.resid)
plt.xlabel('Fitted values')
plt.ylabel('Residuals')
plt.title('Residual Plot')
plt.axhline(y=0, color='r', linestyle='--')
plt.show()

#  QQ图
fig = sm.qqplot(results.resid, stats.norm, fit=True, line='45')
plt.title('QQ Plot')
plt.show()

#  直方图
plt.hist(results.resid, bins=20)
plt.xlabel('Residuals')
plt.ylabel('Frequency')
plt.title('Histogram of Residuals')
plt.show()

# Shapiro-Wilk检验
shapiro_test = stats.shapiro(results.resid)
print(f'Shapiro-Wilk Test: statistic={shapiro_test[0]:.3f}, p-value={shapiro_test[1]:.3f}')

# Breusch-Pagan检验
bp_test = sm.stats.diagnostic.het_breuschpagan(results.resid, results.model.exog)
print(f'Breusch-Pagan Test: LM statistic={bp_test[0]:.3f}, p-value={bp_test[1]:.3f}')

统计模型选择： 线性回归模型假设因变量和自变量之间存在线性关系。然而，如果数据实际上是非线性的，那么线性回归模型可能会产生误导性的结果。此外，线性回归还假设残差是独立同分布的，如果残差存在异方差性或自相关性，那么模型的推断结果可能不可靠。

结果解释与讨论： 观察回归结果的摘要，重点关注R平方值、系数的显著性水平以及残差的诊断结果。如果残差图显示明显的模式（例如，漏斗形），则表明存在异方差性。如果QQ图显示残差偏离正态分布，则表明数据不满足正态性假设。Shapiro-Wilk和Breusch-Pagan检验可以用来定量地检验正态性和异方差性。如果这些假设不满足，则需要考虑使用其他模型，例如广义线性模型或非线性模型。或者，可以尝试对数据进行转换，使其更符合线性回归的假设。

参考文献：
* Statsmodels
* Pandas

2. A/B测试的偏差：选择偏差与辛普森悖论

问题背景： 一家电商平台希望测试一个新的推荐算法是否能够提高用户的点击率。他们将用户随机分为两组：对照组（使用旧算法）和实验组（使用新算法）。

数据描述： 数据集包含用户的ID、分组（对照组或实验组）和点击率。然而，我们忽略了一个重要的因素：用户类型。实际上，平台的用户可以分为两类：新用户和老用户。新用户对推荐算法的敏感度可能高于老用户。

import pandas as pd
import numpy as np
from scipy import stats

# 模拟数据
np.random.seed(0)
n_samples = 200

#  新用户的数据
new_users_control = np.random.normal(0.02, 0.01, n_samples // 2)
new_users_treatment = np.random.normal(0.03, 0.015, n_samples // 2)

#  老用户的数据
old_users_control = np.random.normal(0.05, 0.02, n_samples // 2)
old_users_treatment = np.random.normal(0.04, 0.01, n_samples // 2)

# 创建DataFrame
data = pd.DataFrame({
    '用户类型': ['新用户'] * (n_samples) + ['老用户'] * (n_samples),
    '分组': (['对照组'] * (n_samples // 2) + ['实验组'] * (n_samples // 2)) * 2,
    '点击率': np.concatenate([new_users_control, new_users_treatment, old_users_control, old_users_treatment])
})

# 整体A/B测试
control_overall = data[data['分组'] == '对照组']['点击率']
treatment_overall = data[data['分组'] == '实验组']['点击率']
t_stat_overall, p_val_overall = stats.ttest_ind(control_overall, treatment_overall)
print(f'整体A/B测试 - T统计量: {t_stat_overall:.3f}, P值: {p_val_overall:.3f}')

# 分组A/B测试
for user_type in ['新用户', '老用户']:
    control_group = data[(data['分组'] == '对照组') & (data['用户类型'] == user_type)]['点击率']
    treatment_group = data[(data['分组'] == '实验组') & (data['用户类型'] == user_type)]['点击率']
    t_stat, p_val = stats.ttest_ind(control_group, treatment_group)
    print(f'{user_type} - T统计量: {t_stat:.3f}, P值: {p_val:.3f}')

# 辛普森悖论的可视化
import seaborn as sns

sns.boxplot(x='用户类型', y='点击率', hue='分组', data=data)
plt.title('点击率 by 用户类型 and 分组')
plt.show()

统计模型选择： 通常，我们会使用t检验来比较两组的平均点击率。然而，如果存在选择偏差或辛普森悖论，那么简单的t检验可能会得出错误的结论。选择偏差是指实验组和对照组在某些关键特征上存在差异。辛普森悖论是指在分组数据中观察到的趋势与在总体数据中观察到的趋势相反。

结果解释与讨论： 在这个例子中，我们可能会发现，在总体数据中，实验组的点击率低于对照组。然而，当我们分别分析新用户和老用户时，我们可能会发现，对于新用户，实验组的点击率高于对照组，而对于老用户，实验组的点击率也高于对照组。这就是辛普森悖论。出现这种现象的原因是，实验组中老用户的比例高于对照组，而老用户的点击率普遍较低。因此，我们需要对数据进行分层分析，并考虑用户类型的影响。可以使用更复杂的模型，例如方差分析或回归模型，来控制用户类型的影响。

参考文献：
* Scipy
* 辛普森悖论：统计学中的常见陷阱

3. 时间序列分析的误用：忽略季节性与趋势

问题背景： 假设我们想预测未来一年的月度销售额。我们收集了过去五年的销售数据，并使用时间序列模型进行分析。

数据描述： 数据集包含日期和销售额两个变量。我们使用Python的pandas库加载数据，并使用statsmodels库进行时间序列分析。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.arima.model import ARIMA
from sklearn.metrics import mean_squared_error

# 生成模拟数据
np.random.seed(0)
dates = pd.date_range(start='2017-01-01', end='2021-12-31', freq='MS')
trend = np.linspace(100, 200, len(dates))
seasonal = 50 * np.sin(2 * np.pi * np.arange(len(dates)) / 12)
noise = np.random.normal(0, 10, len(dates))
sales = trend + seasonal + noise

data = pd.DataFrame({'日期': dates, '销售额': sales})
data.set_index('日期', inplace=True)

#  分解时间序列
decomposition = seasonal_decompose(data['销售额'], model='additive')

plt.figure(figsize=(12,8))
plt.subplot(411)
plt.plot(data['销售额'], label='Original')
plt.legend(loc='upper left')
plt.subplot(412)
plt.plot(decomposition.trend, label='Trend')
plt.legend(loc='upper left')
plt.subplot(413)
plt.plot(decomposition.seasonal, label='Seasonal')
plt.legend(loc='upper left')
plt.subplot(414)
plt.plot(decomposition.resid, label='Residuals')
plt.legend(loc='upper left')
plt.tight_layout()
plt.show()

# 划分训练集和测试集
train_data = data[:-12]
test_data = data[-12:]

#  拟合ARIMA模型
model = ARIMA(train_data['销售额'], order=(5, 1, 0)) #  示例参数
model_fit = model.fit()

# 预测
predictions = model_fit.forecast(steps=len(test_data))

# 评估模型
mse = mean_squared_error(test_data['销售额'], predictions)
print(f'Mean Squared Error: {mse:.3f}')

# 预测结果可视化
plt.figure(figsize=(12,6))
plt.plot(test_data['销售额'], label='Actual')
plt.plot(predictions, label='Predicted')
plt.legend(loc='upper left')
plt.title('Sales Forecast')
plt.show()

统计模型选择： 时间序列模型，例如ARIMA模型，可以用来预测未来的销售额。然而，如果时间序列数据存在季节性或趋势，那么简单的ARIMA模型可能无法捕捉到这些模式。在这种情况下，我们需要使用更复杂的模型，例如SARIMA模型，或者先对数据进行季节性调整或趋势分解。

结果解释与讨论： 在这个例子中，我们可以使用seasonal_decompose函数将时间序列分解为趋势、季节性和残差三个部分。通过分析趋势和季节性成分，我们可以更好地理解销售额的变化规律。然后，我们可以使用SARIMA模型来预测未来的销售额，并使用均方误差（MSE）来评估模型的预测精度。此外，我们还需要注意时间序列数据的平稳性。如果数据不平稳，我们需要对其进行差分处理，使其平稳。

参考文献：
* statsmodels时间序列分析
* 时间序列预测：方法与实践

希望这些实例能够帮助读者更严谨地使用Python进行统计分析。记住，统计学不仅仅是运行代码，更重要的是理解统计方法的假设、局限性和潜在偏差。