Python 的 matplotlib + seaborn 实现
举例：
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import numpy as np

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei']
plt.rcParams['axes.unicode_minus'] = False

# 假设你的数据格式如下
# df 包含: x变量, y变量(城镇化率), treat(0或1)

# 示例数据（替换成你的真实数据）
np.random.seed(42)
df = pd.DataFrame({
    '幼年占比': np.concatenate([np.random.uniform(15, 45, 30), np.random.uniform(10, 35, 30)]),
    '壮年占比': np.concatenate([np.random.uniform(50, 75, 30), np.random.uniform(55, 80, 30)]),
    '老年占比': np.concatenate([np.random.uniform(2, 15, 30), np.random.uniform(5, 20, 30)]),
    '城镇化率': np.concatenate([np.random.uniform(35, 75, 30), np.random.uniform(50, 90, 30)]),
    'treat': [0]*30 + [1]*30
})

# 创建三个子图
fig, axes = plt.subplots(1, 3, figsize=(14, 4))

x_vars = ['幼年占比', '壮年占比', '老年占比']
colors = {0: 'steelblue', 1: 'indianred'}
labels = {0: 'treat=0，亚洲国家', 1: 'treat=1，非亚洲国家'}

for ax, x_var in zip(axes, x_vars):
    for treat_val in [0, 1]:
        subset = df[df['treat'] == treat_val]
        # 散点
        ax.scatter(subset[x_var], subset['城镇化率'],
                   c=colors[treat_val], label=labels[treat_val], alpha=0.7, s=30)
        # 拟合线
        z = np.polyfit(subset[x_var], subset['城镇化率'], 1)
        p = np.poly1d(z)
        x_line = np.linspace(subset[x_var].min(), subset[x_var].max(), 100)
        ax.plot(x_line, p(x_line), c=colors[treat_val], linewidth=1.5)
   
    ax.set_xlabel(f'{x_var}(%)')
    ax.set_ylabel('城镇化率(%)')
    ax.legend(fontsize=8)
    ax.axhline(y=60, color='gray', linestyle='--', linewidth=0.8)  # 参考线

plt.suptitle('"一带一路"沿线国家人口年龄结构对城镇化影响的亚洲国家和非亚洲国家处理效应差异', fontsize=12)
plt.tight_layout()
plt.savefig('treatment_effect_plot.png', dpi=300, bbox_inches='tight')
plt.show()