OLS回归-python实现

背景说明：帕累托系数是用来衡量城市规模位序分布的统计量，回归的目的是计算全国31个省份每个省份每一年的帕累托系数，每一个帕累托系数由某个省份某年中所有城市的观测值回归得到

1.回归方程

ln(rank−1/2)=A−ξ×lnsize+μ

变量说明：size是指城市的人口数据，rank是值城市的位序（desc），A是截距，μ是误差

2.数据说明

输入数据：省份-年份-城市-人口                   面板数据（excel文件形式）

期望输出数据：省份-年份-帕累托系数          面板数据（excel文件形式）

3.主要使用的第三方库

statsmodels 是一个 Python 库，用于拟合统计模型、进行统计测试以及统计数据探索和可视化。它是一个非常强大的工具，可以用于众多的统计数据分析任务，从简单的线性回归到复杂的时间序列分析都能胜任

ols回归

import statsmodels.api as sm    
# 进行OLS回归
X = sm.add_constant(x)
model = sm.OLS(Y, X)    # 生成模型
results = model.fit()    # 进行回归并得到回归结果

其中x是自变量，add_constant是添加截距的，这是因为常用的统计库都默认没有截距项（只接受x，y两个参数），所以添加截距项的时候是要给自变量加的，然后生成一个新的自变量，本质是给自变量的旁边再加一列截距列

回归结果的获取
 

results = model.fit()

# 输出模型的系数result.params有结果的各种系数（包括截距和ln_size的系数），
print("截距:", results.params['const'])  # 模型的截距
print("ln_size的系数:", results.params['ln_size'])  # ln_size的系数

# 输出R方值（模型拟合度的衡量指标）
print("R方值:", results.rsquared)

# 输出调整后的R方值（考虑了自变量数量的拟合度衡量）
print("调整后的R方值:", results.rsquared_adj)

# 输出系数的p值（用于检验系数的统计显著性）
print("截距的p值:", results.pvalues['const'])
print("ln_size的p值:", results.pvalues['ln_size'])

# 输出系数的置信区间
print("ln_size的置信区间:", results.conf_int().loc['ln_size'])

# 输出模型的整体统计摘要
print(results.summary())

# 输出模型残差
print(results.resid)

# 使用模型进行预测（假设有新的数据 new_x）
new_x = [/* 新数据 */]
new_data = sm.add_constant(new_x)
predictions = results.predict(new_data)
print("预测值:", predictions)

4.实现代码

import pandas as pd
import numpy as np
import statsmodels.api as sm

def calculate_pareto_coefficient(df):
    """
    计算帕累托系数
    :param df: DataFrame，包含某个省份某年的所有城市及其人口数据
    :return: 帕累托系数
    """
    # 按人口数量降序排序并计算排名
    df = df.sort_values(by='城市区域人口数量', ascending=False)
    df['rank'] = range(1, len(df) + 1)

    # 应用帕累托回归公式
    df['ln_rank'] = np.log(df['rank'] - 0.5)
    df['ln_size'] = np.log(df['城市区域人口数量'])

    # 进行OLS回归
    X = sm.add_constant(df['ln_size'])
    model = sm.OLS(df['ln_rank'], X)
    results = model.fit()

    # 提取帕累托系数
    return -results.params['ln_size']

if __name__ == '__main__':
    # 加载数据
    data = pd.read_excel('省份_年份_城市_人口_panel_data.xlsx')

    # 计算每个省份每年的帕累托系数
    pareto_coefficients = data.groupby(['province', '年份']).apply(calculate_pareto_coefficient)

    # 将Series转换为DataFrame
    pareto_coefficients_df = pareto_coefficients.reset_index()
    pareto_coefficients_df.columns = ['Province', 'Year', 'Pareto_Coefficient']

    # 将结果保存到新的Excel文件
    pareto_coefficients_df.to_excel('pareto_coefficients.xlsx', index=False)

5.代码重点解析

对于代码中的calculate_pareto_coefficient中的回归部分已经解释的比较清楚了，调用函数谁不会啊，最麻烦的还是回归前的数据整理部分。下面详细解释

由于数据是面板数据，所以没有办法直接用dataframe的index和column来直接索引出想要的数据，每一次的回归需要的都是某一省份某一年的所有城市数据，那么我们每一次回归都要筛选出这种维度的数据来，所以要对原始的df进行两次分组，先按省份再按年份，之后每一组都是某个省份某个年份的数据，且分组后每组中包含了匹配的所有数据，形成一个新的df，并不会改变原数据的结构，每组数据都会作为实参传入calculate_pareto_coefficient函数中

# 计算每个省份每年的帕累托系数
pareto_coefficients = data.groupby(['province', '年份']).apply(calculate_pareto_coefficient)