import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import scipy.stats as stats # 修改为完整导入，避免命名冲突
import seaborn as sns
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
import statsmodels.api as sm
import os
import warnings
def analyze_enhanced_data(file_path):
"""
分析增强相关性后的数字乡村与农民收入数据
参数:
file_path: 增强后的Excel文件路径
"""
# 忽略警告
warnings.filterwarnings('ignore')
# 设置中文字体
plt.rcParams['font.sans-serif'] = ['Arial Unicode MS', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
# 读取Excel文件
print(f"正在读取数据: {file_path}")
df = pd.read_excel(file_path)
# 显示数据基本信息
print(f"数据维度: {df.shape}")
print("数据前5行:")
print(df.head())
# 农民收入列名
income_col = "农村居民人均可支配收入(元)"
# 检查列是否存在
if income_col not in df.columns:
print(f"警告: '{income_col}'列不存在，请确认正确的列名")
return None
# 创建结果目录
results_dir = "数字乡村分析结果_增强版"
os.makedirs(results_dir, exist_ok=True)
# 整体相关性分析
print("\n整体数据分析:")
overall_corr = df['数字乡村指数'].corr(df[income_col])
print(f"总体相关系数: {overall_corr:.4f}")
# 进行显著性检验
r, p_value = stats.pearsonr(df['数字乡村指数'], df[income_col])
print(f"皮尔逊相关系数: {r:.4f}, p值: {p_value:.4f}")
if p_value < 0.05:
print("结论: 数字乡村指数与农民收入存在显著相关关系")
else:
print("结论: 数字乡村指数与农民收入无显著相关关系")
# 绘制整体散点图和回归线
plt.figure(figsize=(10, 6))
sns.regplot(x='数字乡村指数', y=income_col, data=df, scatter_kws={'alpha': 0.5}, line_kws={'color': 'red'})
plt.title('Digital Rural Index and Rural Income Relationship')
plt.xlabel('Digital Rural Index')
plt.ylabel('Rural Income per Capita (Yuan)')
plt.grid(True, linestyle='--', alpha=0.7)
plt.savefig(f'{results_dir}/overall_relationship.png', dpi=300, bbox_inches='tight')
plt.close()
# 每年的分析
years = sorted(df['年份'].unique())
results_by_year = []
print("\n按年份分析数字乡村指数与农民收入的关系:")
for year in years:
year_data = df[df['年份'] == year]
if len(year_data) < 5: # 样本太小，跳过
continue
corr = year_data['数字乡村指数'].corr(year_data[income_col])
results_by_year.append({'年份': year, '相关系数': corr})
print(f"{year}年: 相关系数 = {corr:.4f}, 样本数 = {len(year_data)}")
# 为每年生成散点图
plt.figure(figsize=(8, 6))
sns.regplot(x='数字乡村指数', y=income_col, data=year_data,
scatter_kws={'alpha': 0.6}, line_kws={'color': 'red'})
plt.title(f'Relationship between Digital Rural Index and Income ({year})')
plt.xlabel('Digital Rural Index')
plt.ylabel('Rural Income per Capita (Yuan)')
plt.grid(True, linestyle='--', alpha=0.7)
plt.savefig(f'{results_dir}/annual_relationship_{year}.png', dpi=300, bbox_inches='tight')
plt.close()
# 生成相关系数热力图
corr_by_year_df = pd.DataFrame(results_by_year).set_index('年份')
plt.figure(figsize=(12, 6))
sns.heatmap(pd.DataFrame(corr_by_year_df['相关系数']).T,
cmap='RdYlGn', annot=True, fmt='.3f', vmin=0, vmax=1)
plt.title('Correlation Coefficient by Year')
plt.savefig(f'{results_dir}/correlation_heatmap.png', dpi=300, bbox_inches='tight')
plt.close()
# 汇总每年的相关系数变化
results_df = pd.DataFrame(results_by_year)
if not results_df.empty:
plt.figure(figsize=(10, 6))
plt.plot(results_df['年份'], results_df['相关系数'], marker='o', linestyle='-', color='blue')
plt.axhline(y=0, color='gray', linestyle='--', alpha=0.7)
plt.title('Annual Changes in Correlation Coefficient')
plt.xlabel('Year')
plt.ylabel('Correlation Coefficient')
plt.grid(True, alpha=0.7)
# 添加趋势线
z = np.polyfit(results_df['年份'], results_df['相关系数'], 1)
p = np.poly1d(z)
plt.plot(results_df['年份'], p(results_df['年份']), "r--", alpha=0.7,
label=f"Trend line (slope: {z[0]:.4f})")
plt.legend()
plt.savefig(f'{results_dir}/correlation_coefficient_changes.png', dpi=300, bbox_inches='tight')
plt.close()
# 构建回归模型
print("\n线性回归分析:")
X = df[['数字乡村指数']]
y = df[income_col]
# 添加常数项
X_with_const = sm.add_constant(X)
model = sm.OLS(y, X_with_const).fit()
print(model.summary())
# 提取关键结果
r_squared = model.rsquared
adj_r_squared = model.rsquared_adj
f_value = model.fvalue
f_pvalue = model.f_pvalue
coefficients = model.params
print(f"\n决定系数(R²): {r_squared:.4f}")
print(f"调整后的决定系数: {adj_r_squared:.4f}")
print(f"F统计量: {f_value:.4f}, p值: {f_pvalue:.4f}")
print(f"回归方程: {income_col} = {coefficients['const']:.4f} + {coefficients['数字乡村指数']:.4f} × 数字乡村指数")
# 计算弹性系数(半弹性)
mean_digital = df['数字乡村指数'].mean()
mean_income = df[income_col].mean()
elasticity = coefficients['数字乡村指数'] * (mean_digital / mean_income)
print(f"\n弹性系数: {elasticity:.4f}")
print(f"解释: 数字乡村指数增加1%时，农民人均收入平均增加约{elasticity:.4f}%")
# 分年度回归分析
print("\n分年度回归分析(每5年一组):")
year_groups = []
for i in range(0, len(years), 5):
year_groups.append(years[i:i + 5])
group_results = []
for group in year_groups:
if len(group) > 0:
group_data = df[df['年份'].isin(group)]
if len(group_data) > 30: # 确保样本充足
X_group = group_data[['数字乡村指数']]
X_group = sm.add_constant(X_group)
y_group = group_data[income_col]
try:
group_model = sm.OLS(y_group, X_group).fit()
years_str = f"{min(group)}-{max(group)}"
print(f"\n{years_str}年度组回归结果:")
print(f"R²: {group_model.rsquared:.4f}, p值: {group_model.f_pvalue:.4f}")
print(
f"系数: {group_model.params['数字乡村指数']:.4f}, p值: {group_model.pvalues['数字乡村指数']:.4f}")
group_results.append({
'年份组': years_str,
'R²': group_model.rsquared,
'系数': group_model.params['数字乡村指数'],
'p值': group_model.pvalues['数字乡村指数']
})
except Exception as e:
print(f"年度组{group}回归失败: {e}")
# 绘制分组回归结果比较图
if group_results:
group_df = pd.DataFrame(group_results)
fig, ax1 = plt.figure(figsize=(10, 6)), plt.gca()
# 绘制R²变化
ax1.plot(group_df['年份组'], group_df['R²'], 'o-', color='blue', label='R²')
ax1.set_xlabel('Year Groups')
ax1.set_ylabel('R²', color='blue')
ax1.tick_params(axis='y', labelcolor='blue')
# 创建第二个y轴
ax2 = ax1.twinx()
ax2.plot(group_df['年份组'], group_df['系数'], 'o-', color='red', label='Coefficient')
ax2.set_ylabel('Coefficient', color='red')
ax2.tick_params(axis='y', labelcolor='red')
plt.title('Changes in R² and Coefficient over Time')
fig.tight_layout()
plt.savefig(f'{results_dir}/group_regression_results.png', dpi=300, bbox_inches='tight')
plt.close()
# 绘制预测图
plt.figure(figsize=(10, 6))
# 散点图
plt.scatter(df['数字乡村指数'], df[income_col], alpha=0.5, color='blue', label='Observed Data')
# 排序预测值以便绘制平滑曲线
sorted_idx = np.argsort(df['数字乡村指数'].values)
sorted_x = df['数字乡村指数'].values[sorted_idx]
predicted_y = model.predict(X_with_const).values[sorted_idx]
plt.plot(sorted_x, predicted_y, color='red', linewidth=2, label='Regression Line')
# 添加置信区间
try:
# 使用预测区间代替置信区间
pred_ols = model.get_prediction(X_with_const.iloc[sorted_idx])
iv_l = pred_ols.summary_frame()['obs_ci_lower']
iv_u = pred_ols.summary_frame()['obs_ci_upper']
plt.fill_between(sorted_x, iv_l, iv_u, color='red', alpha=0.1, label='95% Prediction Interval')
except Exception as e:
print(f"无法计算置信区间: {e}")
plt.title('Digital Rural Index and Rural Income with Prediction')
plt.xlabel('Digital Rural Index')
plt.ylabel('Rural Income per Capita (Yuan)')
plt.grid(True, linestyle='--', alpha=0.7)
plt.legend()
plt.savefig(f'{results_dir}/regression_prediction.png', dpi=300, bbox_inches='tight')
plt.close()
# 绘制残差分析图
plt.figure(figsize=(12, 10))
# 残差与拟合值
plt.subplot(2, 2, 1)
plt.scatter(model.fittedvalues, model.resid, alpha=0.5)
plt.axhline(y=0, color='r', linestyle='--')
plt.xlabel('Fitted Values')
plt.ylabel('Residuals')
plt.title('Residuals vs Fitted')
plt.grid(True, linestyle='--', alpha=0.7)
# 残差QQ图
plt.subplot(2, 2, 2)
# 这里使用已导入的stats模块，避免重复导入
stats.probplot(model.resid, plot=plt)
plt.title('Q-Q Plot of Residuals')
# 残差直方图
plt.subplot(2, 2, 3)
plt.hist(model.resid, bins=30, alpha=0.7, color='blue', edgecolor='black')
plt.xlabel('Residuals')
plt.ylabel('Frequency')
plt.title('Histogram of Residuals')
plt.grid(True, linestyle='--', alpha=0.7)
# 残差与自变量
plt.subplot(2, 2, 4)
plt.scatter(df['数字乡村指数'], model.resid, alpha=0.5)
plt.axhline(y=0, color='r', linestyle='--')
plt.xlabel('Digital Rural Index')
plt.ylabel('Residuals')
plt.title('Residuals vs Digital Rural Index')
plt.grid(True, linestyle='--', alpha=0.7)
plt.tight_layout()
plt.savefig(f'{results_dir}/residual_analysis.png', dpi=300, bbox_inches='tight')
plt.close()
# 总结结论
print("\n研究结论:")
if r_squared > 0.3:
print("1. 数据分析表明数字乡村建设对农民收入有显著的强相关影响。")
print(f"2. 数字乡村指数每提高1个单位，农民人均收入平均增加{coefficients['数字乡村指数']:.2f}元。")
print(f"3. 数字乡村发展水平能解释农民收入变化的{r_squared * 100:.2f}%，表明数字乡村是影响农民收入的重要因素。")
print("4. 建议：大力推进数字乡村建设，加大信息基础设施投入，促进农村数字化转型，有助于提高农民收入水平。")
elif r_squared > 0.1:
print("1. 数据分析表明数字乡村建设对农民收入有明显的积极影响。")
print(f"2. 数字乡村指数每提高1个单位，农民人均收入平均增加{coefficients['数字乡村指数']:.2f}元。")
print(f"3. 数字乡村发展水平能解释农民收入变化的{r_squared * 100:.2f}%。")
print("4. 建议：继续推进数字乡村建设，结合其他农村发展政策，综合促进农民收入增长。")
else:
print("1. 数据分析表明数字乡村建设与农民收入存在一定相关关系。")
print(f"2. 数字乡村指数每提高1个单位，农民人均收入平均增加{coefficients['数字乡村指数']:.2f}元。")
print(f"3. 虽然解释力有限，但存在统计显著性，说明数字乡村对农民收入有一定积极影响。")
print("4. 建议：将数字乡村建设作为提升农民收入的辅助手段，需配合其他更直接的增收措施。")
# 保存分析结果
summary_file = f"{results_dir}/分析报告.txt"
with open(summary_file, 'w', encoding='utf-8') as f:
f.write("数字乡村与农民收入关系分析报告\n")
f.write("=" * 50 + "\n\n")
f.write(f"1. 总体相关系数: {overall_corr:.4f}\n")
f.write(f"2. 显著性检验: p值 = {p_value:.4f}\n")
f.write(f"3. 决定系数(R²): {r_squared:.4f}\n")
f.write(
f"4. 回归方程: {income_col} = {coefficients['const']:.4f} + {coefficients['数字乡村指数']:.4f} × 数字乡村指数\n")
f.write(f"5. 弹性系数: {elasticity:.4f}\n\n")
if p_value < 0.05:
f.write("结论: 数字乡村与农民收入之间存在显著相关关系\n")
f.write(f"数字乡村指数每提高1个单位，农民人均收入平均增加{coefficients['数字乡村指数']:.2f}元\n")
else:
f.write("结论: 未能证明数字乡村与农民收入之间存在显著关系\n")
# 添加年度分析摘要
f.write("\n年度相关系数分析:\n")
for result in results_by_year:
f.write(f"{result['年份']}年: {result['相关系数']:.4f}\n")
# 添加年度组回归结果
if group_results:
f.write("\n年度组回归分析:\n")
f.write("年份组\tR²\t系数\tp值\n")
for result in group_results:
f.write(f"{result['年份组']}\t{result['R²']:.4f}\t{result['系数']:.4f}\t{result['p值']:.4f}\n")
f.write("\n政策建议:\n")
if r_squared > 0.3:
f.write("1. 大力推进数字乡村建设，加大农村信息基础设施投入\n")
f.write("2. 完善农村互联网应用体系，提升农民数字技能\n")
f.write("3. 发展农村电子商务，拓宽农产品销售渠道\n")
f.write("4. 建设智慧农业系统，提高农业生产效率\n")
f.write("5. 数字金融赋能，解决农村融资难问题\n")
elif r_squared > 0.1:
f.write("1. 有序推进数字乡村建设，提高农村数字化水平\n")
f.write("2. 加强农民数字技能培训，缩小城乡数字鸿沟\n")
f.write("3. 结合产业政策，利用数字技术助力农业升级\n")
f.write("4. 发展适合本地特色的数字应用，提高实用性\n")
else:
f.write("1. 将数字乡村建设作为农村发展的辅助手段\n")
f.write("2. 优先发展更直接影响农民收入的产业与政策\n")
f.write("3. 探索适合不同地区的数字乡村发展模式\n")
f.write("4. 研究数字乡村对农民收入的长期影响机制\n")
print(f"\n分析完成，结果已保存到 {results_dir} 目录")
# 返回分析结果
results = {
'df': df,
'correlation': r,
'r_squared': r_squared,
'model': model,
'year_correlations': results_by_year,
'group_results': group_results if 'group_results' in locals() else None
}
return results
if __name__ == "__main__":
# 分析增强后的数据
file_path = "统计年鉴_增强相关性_方法2.xlsx" # 替换为实际文件路径
results = analyze_enhanced_data(file_path)
# 也可以分析方法2增强的数据
# file_path2 = "统计年鉴_增强相关性_方法2.xlsx"
# results2 = analyze_enhanced_data(file_path2)这个代码要加入变量分析应该加入到哪里
最新发布
本文介绍了Python中head()函数的使用，包括其语法、工作原理和应用示例。head()函数通过生成器表达式快速读取文件的前n行，避免遍历整个文件，提高效率。文中提供了一个示例，展示如何读取文件的前三行。
扫一扫
516
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
