在数据分析的旅程中，我们曾一同探索Statsmodels库的基本用法，如何进行简单线性回归分析。然而，现在是时候深入了解更为复杂、多变的数据关系了。本文将引领你进入Statsmodels库的进阶用法，涵盖多元线性回归、Logistic回归以及时间序列分析，让我们更深刻地认识数据背后的关系。
多元线性回归是一种用于建模和分析多个自变量与一个因变量之间关系的统计技术。与简单线性回归不同，多元线性回归考虑了多个自变量对因变量的影响。模型的基本形式为：
Y = β 0 + β 1 X 1 + β 2 X 2 + … + β n X n + ϵ Y = \beta_0 + \beta_1X_1 + \beta_2X_2 + \ldots + \beta_nX_n + \epsilon Y = β 0 + β 1 X 1 + β 2 X 2 + … + β n X n + ϵ
其中， Y Y Y 是因变量， X 1 , X 2 , … , X n X_1, X_2, \ldots, X_n X 1 , X 2 , … , X n 是自变量， β 0 , β 1 , … , β n \beta_0, \beta_1, \ldots, \beta_n β 0 , β 1 , … , β n 是回归系数， ϵ \epsilon ϵ 是误差。
假设我们有一份数据集，其中包含房屋的面积（SquareFeet）、卧室的数量（Bedrooms）以及距离市中心的距离（Distance）等多个特征。我们想要建立一个模型，预测房屋的价格（Price）。这是一个典型的多元线性回归问题。
我们利用stasmodels进行建模分析：
import numpy as np
import pandas as pd
import statsmodels.api as sm
# 生成模拟数据
np.random.seed(12)
data = {
'SquareFeet': np.random.randint(800, 2500, 50),
'Bedrooms': np.random.randint(1, 5, 50),
'Distance': np.random.uniform(1, 20, 50)
}
df = pd.DataFrame(data)
df['Price'] = df['SquareFeet'] * 50 + df['Bedrooms'] * 10000 - df['Distance'] * 3000 + (np.random.rand(50,1) * 100).reshape(50) * 200
# 添加常数列
X = sm.add_constant(df[['SquareFeet', 'Bedrooms', 'Distance']])
y = df['Price']
# 拟合多元线性回归模型
model = sm.OLS(y, X).fit()
# 打印模型摘要
print(model.summary())
运行上述代码后，我们可以得到下面的结果：
根据输出信息，逐项来评估模型的拟合效果：
R-squared (R²) :
Adjusted R-squared (调整后R²) :
F-statistic (F统计量) :
各个自变量的系数（coef） :
常数项（const） :
诊断信息 :
综合来看，这个模型在统计上表现出色，大部分自变量对因变量有显著影响。但需要进一步检验误差的正态性、共线性等问题，并且注意到常数项可能不是显著的。如果误差项不符合模型假设，或者存在共线性，可能需要对模型进行修正或改进。
预测房价
假设现在我们利用上述的模型，进行房价预测，那么我么可以利用下面的代码进行实现：
import numpy as np
import pandas as pd
import statsmodels.api as sm
# 要预测的数据
new_data = {
'SquareFeet': [2000,1000], # 新的面积
'Bedrooms': [3,2], # 卧室数量
'Distance': [15,5] # 距离
}
# 转换为DataFrame
new_df = pd.DataFrame(new_data)
# 添加常数列
new_df = sm.add_constant(new_df)
# 确保新数据列的顺序与模型训练时的列顺序相同
new_df = new_df[['const', 'SquareFeet', 'Bedrooms', 'Distance']]
# 使用模型进行预测
predicted_price = model.predict(new_df)
print("预测价格为:", predicted_price)
利用上述代码后，我们可以预测未来的房价情况。
Logistic回归是一种用于解决分类问题的机器学习算法。尽管名字中带有"回归"，但实际上它是一种分类算法，用于预测离散型的输出，通常是二元分类（有两个类别）问题。
这个算法的基本思想是利用一个逻辑函数（Logistic函数）来建立一个模型，该模型能够估计特征与一个事件发生之间的概率关系。Logistic函数能够把任意实数值映射到0和1之间，其形式为 σ ( z ) = 1 1 + e − z \sigma(z) = \frac{1}{1 + e^{-z}} σ ( z ) = 1 + e − z 1 其中，z是输入特征的线性组合，通过特征的加权求和得到。
Logistic回归的训练过程就是通过拟合这个逻辑函数来找到最佳的参数（权重），使得模型能够最好地描述已知数据的特征与类别之间的关系。通常使用梯度下降等优化算法来最小化损失函数，以找到最优的参数值。
这种算法适用于二分类问题，比如预测学生是否通过考试、邮件是否为垃圾邮件等。虽然最初用于处理线性可分问题，但它也可以通过特征工程或引入多项式特征来解决非线性分类问题。
当建立逻辑回归模型时，首先需要将数据整理成适合建模的格式，并对分类变量进行编码（比如独热编码）。
独热编码（One-Hot Encoding)
是一种用来表示分类变量的技术。在机器学习和数据处理中，经常会遇到分类数据，例如颜色、类别、标签等。独热编码将这些分类数据转换为向量形式，其中每个分类都表示为一个向量，该向量的长度等于分类的总数，其中只有一个元素为1，其他元素为0。这个1所在的位置表示原始数据中的分类。
举个例子，假设有三种水果：苹果、橙子和香蕉。使用独热编码，苹果可以表示为[1, 0, 0]，橙子可以表示为[0, 1, 0]，香蕉可以表示为[0, 0, 1]。这样表示的好处是，每个分类之间彼此独立，不会产生数值大小上的误解或关联。
下面是一个示例，假设你已经有了包含广告 点击信息(Clicked) 、 广告质量(AdQuality) 、 广告位置(AdPlacement) 、 受众年龄(AudienceAge) 、 用户设备(UserDevice) 、 用户平台(UserPlatform) 这些信息，建立一个模型用来预测用户对广告的点击情况。
import pandas as pd
import numpy as np
import statsmodels.api as sm
# 生成示例数据
np.random.seed(15)
# 创建一个包含广告点击信息和其他特征的示例数据集
num_samples = 1000
data = {
'AdQuality': np.random.randint(1, 11, num_samples),
'AdPlacement': np.random.choice(['Header', 'Sidebar', 'Footer'], num_samples),
'AudienceAge': np.random.choice(['18-25', '25-35', '35-45'], num_samples),
'UserDevice': np.random.choice(['Mobile', 'Desktop'], num_samples),
'UserPlatform': np.random.choice(['iOS', 'Android', 'Windows'], num_samples),
'Clicked': np.random.randint(0, 2, num_samples) # 二元点击情况，0表示未点击，1表示点击
}
df = pd.DataFrame(data)
# 对分类变量进行独热编码前的列名
original_columns = df.columns.tolist()
# 对分类变量进行独热编码，并将数据类型转换为整数型
df_encoded = pd.get_dummies(df, columns=['AdPlacement', 'AudienceAge', 'UserDevice', 'UserPlatform'], drop_first=True)
df_encoded = df_encoded.astype(int) # 转换数据类型为整数型
# 对分类变量进行独热编码后的列名
encoded_columns = df_encoded.columns.tolist()
# 准备数据
X = df_encoded.drop('Clicked', axis=1) # 自变量
y = df_encoded['Clicked'] # 因变量
# 添加常数列
X = sm.add_constant(X)
# 构建逻辑回归模型
logit_model = sm.Logit(y, X)
# 拟合模型
result = logit_model.fit()
# 打印模型摘要
print(result.summary())
# 生成预测数据
new_data = {
'AdQuality': [8, 5, 9, 6, 7],
'AdPlacement': ['Header', 'Sidebar', 'Footer', 'Header', 'Footer'],
'AudienceAge': ['25-35', '18-25', '35-45', '18-25', '25-35'],
'UserDevice': ['Mobile', 'Desktop', 'Mobile', 'Desktop', 'Desktop'],
'UserPlatform': ['Android', 'iOS', 'Windows', 'Android', 'Windows']
}
# 将测试数据转换为DataFrame并进行独热编码
new_df = pd.DataFrame(new_data)
new_df_encoded = pd.get_dummies(new_df, columns=original_columns[1:-1], drop_first=True)
new_df_encoded = new_df_encoded.astype(int) # 转换数据类型为整数型
# 输出独热编码前后的列名对照关系
# print("独热编码前的列名：", original_columns)
# print("独热编码后的列名：", new_df_encoded.columns.tolist())
# 添加常数列
new_df_encoded = sm.add_constant(new_df_encoded)
# 使用模型进行预测
predicted_click = result.predict(new_df_encoded)
print("预测的点击概率为:", predicted_click)
运行上述结果后，我们可以得到下面的结果：
根据上方输出的逻辑回归结果，我们可以做出以下结论：
模型拟合评估 ：
各变量对点击率的影响 ：
预测结果 ：
结论
：
综上所述，根据模型拟合结果和预测效果来看，该模型并未能很好地捕捉到广告点击与其他特征之间的关系。模型的预测能力较差，无法有效区分点击和未点击的情况。可能需要更多的数据或者重新选择和提取特征来改善模型的预测效果。
由于数据是随机生成的，所以建模效果不理想也在理想之中。如果我们在实际场景中，如果出现模型效果不理想的情况，可以考虑下面的方案：
特征工程 ：可能需要重新选择、组合或提取特征，以更好地捕获广告点击的规律。这可能包括：
模型调参 ：可能当前的模型参数配置并不是最优的，可以尝试调整模型的超参数，或者尝试不同的模型类型，比如随机森林、梯度提升树等，以及集成学习模型，如XGBoost、LightGBM等。
解决样本不平衡问题 ：如果点击和未点击的样本不平衡，可以考虑使用过采样或欠采样等技术来平衡数据集，避免模型对多数类别过度拟合。
交叉验证与验证集 ：使用交叉验证来评估模型的稳定性和泛化能力，并使用独立的验证集来验证模型在未见过的数据上的表现。
特征交互与多项式特征 ：尝试将特征进行交互或者多项式展开，以捕捉特征之间的复杂关系。
集成学习方法 ：尝试使用集成学习方法（如随机森林、Adaboost、Gradient Boosting等）来组合多个基本模型，以提高整体的预测能力。
数据增强 ：如果可能，尝试收集更多的数据，特别是具有代表性的样本，以提高模型的泛化能力。
特征缩放和正则化 ：某些模型可能对特征的尺度敏感，尝试对特征进行缩放，同时对于线性模型，考虑加入正则化项来防止过拟合。
Statsmodels 在时间序列分析中提供了丰富的功能和模型，主要用于建模、分析和预测时间序列数据。
ARIMA 模型 ：
单位根检验（ADF 检验） ：
自相关函数和偏自相关函数 ：
时间序列拟合和预测 ：
季节性调整 ：
时间序列分解 ：
平稳性和季节性检验 ：
时间序列数据可视化 ：
这些功能和工具可以帮助分析者对时间序列数据进行探索性分析、建模和预测，更好地理解数据的结构和特征，并做出合理的预测和决策。
时序数据包含着时间的信息，通过Statsmodels库，我们可以进行更为细致的时间序列分析。本节将以实例为例，展示如何使用Statsmodels库绘制时间序列图、识别趋势，并解读时序数据中的关键特征。
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
# 生成销售数据
np.random.seed(12)
dates = pd.date_range('2018-01-01', periods=72, freq='M')
sales = 1000 + np.random.normal(0, 200, 72) + 100 * np.sin(np.arange(72) * 2 * np.pi / 12)
sales_data = pd.DataFrame({'Date': dates, 'Sales': sales})
sales_data = sales_data.set_index('Date')
# 绘制销售数据的时间序列图
plt.figure(figsize=(10, 6))
plt.plot(sales_data, label='Sales Data')
plt.xlabel('Date')
plt.ylabel('Sales')
plt.title('Sales Time Series')
plt.legend()
plt.show()
# SARIMA 模型拟合
model = sm.tsa.statespace.SARIMAX(sales_data, order=(1, 1, 1), seasonal_order=(1, 1, 1, 12))
results = model.fit()
# 打印模型摘要
print(results.summary())
# 绘制拟合值和原始数据的对比图
plt.figure(figsize=(10, 6))
plt.plot(sales_data, label='Sales Data')
plt.plot(results.fittedvalues, color='red', label='Fitted Values')
plt.xlabel('Date')
plt.ylabel('Sales')
plt.title('Sales Data vs Fitted SARIMA Model')
plt.legend()
plt.show()
# 进行未来销售值的预测
forecast_steps = 12 # 预测未来 12 个月的销售值
forecast = results.get_forecast(steps=forecast_steps)
# 获取预测结果
forecast_values = forecast.predicted_mean
forecast_ci = forecast.conf_int()
# 打印预测值和置信区间
print("Forecasted Sales for the next 12 months:")
print(forecast_values)
print("\nForecast Confidence Intervals:")
print(forecast_ci)
运行上述代码后，结果如下：
从上述模型的结果来看：
系数显著性 ：
模型拟合效果 ：
残差诊断 ：
异方差性 ：
综合来看，这个 SARIMA 模型在拟合季节性时间序列数据方面表现得相对良好。虽然其中一些系数并不显著，但模型整体的拟合效果仍然较为可靠，能够较好地捕捉数据的季节性特征和趋势。
在本文中，我们通过Statsmodels库深度挖掘了数据之间的关系。多元线性回归让我们能够同时考虑多个自变量，更全面地理解因变量的变化。Logistic回归为我们打开了处理二分类问题的大门，预测概率在实际场景中具有广泛的应用。最后，通过时间序列分析，我们更深刻地认识了时序数据的动态特性。
Statsmodels库在不同数据关系的探索中表现出色，无论是线性关系、概率估计还是时间序列分析，都能得心应手。在你的数据旅程中，不妨深入研究Statsmodels库的更多用法，挖掘数据更丰富的内涵。希望本文能为你在数据分析的路上指明方向，让你更加游刃有余地驾驭复杂的数据关系。