线性回归模型

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
import warnings
warnings.filterwarnings('ignore')

# 设置中文显示
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 加载数据（使用波士顿房价数据集）
from sklearn.datasets import fetch_california_housing
housing = fetch_california_housing()
df = pd.DataFrame(housing.data, columns=housing.feature_names)
df['房价'] = housing.target * 100000  # 转换为实际美元价格

print("数据集基本信息:")
print(f"样本数: {df.shape[0]}")
print(f"变量数: {df.shape[1]}")
print("\n变量名:")
print(df.columns.tolist())
print("\n前5行数据:")
print(df.head())
print("\n描述性统计:")
print(df.describe())

# 1.1 绘制散点图矩阵
print("\n1. 绘制散点图矩阵...")
fig = plt.figure(figsize=(15, 10))
sns.pairplot(df[['房价', 'MedInc', 'HouseAge', 'AveRooms', 'AveBedrms', 'Population']])
plt.suptitle('散点图矩阵：因变量与主要自变量的关系', y=1.02)
plt.tight_layout()
plt.show()

# 1.2 相关系数矩阵
plt.figure(figsize=(10, 8))
corr_matrix = df.corr()
sns.heatmap(corr_matrix, annot=True, cmap='coolwarm', center=0,
            square=True, linewidths=0.5, cbar_kws={"shrink": 0.8})
plt.title('变量相关系数矩阵')
plt.tight_layout()
plt.show()

# 1.3 异常值检测
print("\n2. 异常值检测与处理...")

# 使用Z-score检测异常值
z_scores = np.abs(stats.zscore(df))
outliers = (z_scores > 3).any(axis=1)
print(f"使用Z-score > 3检测到的异常值数量: {outliers.sum()}")

# 使用箱线图可视化异常值
fig, axes = plt.subplots(2, 4, figsize=(16, 8))
axes = axes.ravel()
for i, col in enumerate(df.columns[:8]):
    df.boxplot(column=col, ax=axes[i])
    axes[i].set_title(f'{col}的箱线图')
plt.tight_layout()
plt.show()

# 使用Cook距离检测高影响点（先建立简单模型）
X_temp = sm.add_constant(df.drop('房价', axis=1))
y_temp = df['房价']
model_temp = sm.OLS(y_temp, X_temp).fit()
influence = model_temp.get_influence()
cooks_d = influence.cooks_distance[0]

# 识别高影响点
cook_threshold = 4 / len(df)
high_influence = cooks_d > cook_threshold
print(f"基于Cook距离检测到的高影响点数量: {high_influence.sum()}")

# 可视化Cook距离
plt.figure(figsize=(10, 6))
plt.scatter(range(len(cooks_d)), cooks_d, alpha=0.6)
plt.axhline(y=cook_threshold, color='r', linestyle='--', label=f'阈值={cook_threshold:.4f}')
plt.xlabel('观测序号')
plt.ylabel('Cook距离')
plt.title('Cook距离图 - 检测高影响点')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

# 1.4 处理异常值（这里选择不删除，但会标记）
df['异常值标记'] = outliers | high_influence
print(f"\n总异常值/高影响点数量: {df['异常值标记'].sum()}")
print("注：实际分析中，我们会审查这些点，根据业务逻辑决定是否处理")

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80

# 2.1 数据准备
# 分离特征和目标变量
X = df.drop(['房价', '异常值标记'], axis=1)
y = df['房价']

# 添加常数项（截距）
X = sm.add_constant(X)

# 2.2 建立全模型
print("\n3. 建立全模型（包含所有自变量）...")
full_model = sm.OLS(y, X).fit()
print("全模型摘要:")
print(full_model.summary())

# 2.3 变量选择 - 逐步回归
print("\n4. 变量选择（逐步回归）...")

# 前向选择函数
def forward_selection(X, y, significance_level=0.05):
    initial_features = X.columns.tolist()
    initial_features.remove('const')  # 保留常数项
    selected_features = ['const']

    while len(initial_features) > 0:
        remaining_features = list(set(initial_features) - set(selected_features))
        pvalues = []

        for feature in remaining_features:
            # 临时模型
            features = selected_features + [feature]
            X_temp = X[features]
            model = sm.OLS(y, X_temp).fit()
            pvalue = model.pvalues[feature]
            pvalues.append((feature, pvalue))

        # 选择p值最小的特征
        pvalues.sort(key=lambda x: x[1])
        best_feature, best_pvalue = pvalues[0]

        if best_pvalue < significance_level:
            selected_features.append(best_feature)
            print(f"添加变量: {best_feature} (p值: {best_pvalue:.4f})")
        else:
            break

    return selected_features

# 执行前向选择
selected_vars = forward_selection(X, y)
print(f"\n最终选择的变量: {selected_vars}")

# 2.4 建立精简模型
X_selected = X[selected_vars]
reduced_model = sm.OLS(y, X_selected).fit()
print("\n精简模型摘要:")
print(reduced_model.summary())

# 2.5 使用AIC/BIC进行模型比较
print("\n5. 模型比较（基于信息准则）:")
print(f"全模型 AIC: {full_model.aic:.2f}, BIC: {full_model.bic:.2f}")
print(f"精简模型 AIC: {reduced_model.aic:.2f}, BIC: {reduced_model.bic:.2f}")

# 2.6 使用LASSO进行变量选择（可选）
from sklearn.linear_model import LassoCV

# 标准化数据
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X.drop('const', axis=1))

# LASSO交叉验证
lasso_cv = LassoCV(cv=5, random_state=42, max_iter=10000)
lasso_cv.fit(X_scaled, y)

print(f"\nLASSO选择的最优alpha: {lasso_cv.alpha_:.4f}")
print("LASSO选择的变量（系数非零）:")
lasso_coef = pd.DataFrame({
    '变量': X.drop('const', axis=1).columns,
    '系数': lasso_cv.coef_
})
print(lasso_coef[lasso_coef['系数'] != 0])

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167

print("\n" + "="*60)
print("6. 模型假设诊断")
print("="*60)

# 获取预测值和残差
y_pred = reduced_model.fittedvalues
residuals = reduced_model.resid

# 3.1 线性与独立性诊断
print("\n6.1 线性与独立性诊断...")
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# 残差 vs 拟合值图
axes[0].scatter(y_pred, residuals, alpha=0.6)
axes[0].axhline(y=0, color='r', linestyle='--')
axes[0].set_xlabel('拟合值')
axes[0].set_ylabel('残差')
axes[0].set_title('残差 vs 拟合值图（检验线性与独立性）')
axes[0].grid(True, alpha=0.3)

# 添加Lowess平滑线看趋势
import statsmodels.nonparametric.smoothers_lowess as lowess
lowess_line = lowess.lowess(residuals, y_pred, frac=0.3)
axes[0].plot(lowess_line[:, 0], lowess_line[:, 1], color='green', linewidth=2, label='Lowess平滑线')
axes[0].legend()

# 残差 vs 自变量图（示例：收入）
if 'MedInc' in selected_vars:
    axes[1].scatter(X['MedInc'], residuals, alpha=0.6)
    axes[1].axhline(y=0, color='r', linestyle='--')
    axes[1].set_xlabel('收入中位数 (MedInc)')
    axes[1].set_ylabel('残差')
    axes[1].set_title('残差 vs 自变量（检验线性关系）')
    axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# 3.2 同方差性检验
print("\n6.2 同方差性检验...")
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# 残差 vs 拟合值（看散点带宽度）
axes[0].scatter(y_pred, np.abs(residuals), alpha=0.6)
axes[0].set_xlabel('拟合值')
axes[0].set_ylabel('残差绝对值')
axes[0].set_title('残差绝对值 vs 拟合值（检验异方差）')
axes[0].grid(True, alpha=0.3)

# 添加趋势线
z = np.polyfit(y_pred, np.abs(residuals), 1)
p = np.poly1d(z)
axes[0].plot(np.sort(y_pred), p(np.sort(y_pred)), color='red', linewidth=2)

# 正式检验：Breusch-Pagan检验
from statsmodels.stats.diagnostic import het_breuschpagan
bp_test = het_breuschpagan(residuals, reduced_model.model.exog)
bp_labels = ['LM统计量', 'LM检验p值', 'F统计量', 'F检验p值']
bp_result = dict(zip(bp_labels, bp_test))

axes[1].bar(range(4), [bp_test[0], bp_test[1], bp_test[2], bp_test[3]],
            color=['blue', 'red', 'blue', 'red'])
axes[1].set_xticks(range(4))
axes[1].set_xticklabels(['LM统计量', 'LM-p值', 'F统计量', 'F-p值'], rotation=45)
axes[1].set_title('Breusch-Pagan检验结果')
axes[1].axhline(y=0.05, color='green', linestyle='--', label='显著性水平0.05')
axes[1].legend()

plt.tight_layout()
plt.show()

print(f"Breusch-Pagan检验结果:")
print(f"  LM统计量: {bp_test[0]:.4f}, p值: {bp_test[1]:.4f}")
print(f"  F统计量: {bp_test[2]:.4f}, p值: {bp_test[3]:.4f}")
if bp_test[1] > 0.05:
    print("  结论：不能拒绝同方差原假设（不存在显著异方差）")
else:
    print("  警告：存在异方差问题，可能需要使用加权最小二乘法")

# 3.3 正态性检验
print("\n6.3 正态性检验...")
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# Q-Q图
sm.qqplot(residuals, line='45', fit=True, ax=axes[0])
axes[0].set_title('正态Q-Q图（检验残差正态性）')

# 残差直方图
axes[1].hist(residuals, bins=30, edgecolor='black', alpha=0.7, density=True)
axes[1].set_xlabel('残差')
axes[1].set_ylabel('密度')

# 添加正态分布曲线
from scipy.stats import norm
mu, std = norm.fit(residuals)
xmin, xmax = axes[1].get_xlim()
x = np.linspace(xmin, xmax, 100)
p = norm.pdf(x, mu, std)
axes[1].plot(x, p, 'r', linewidth=2, label=f'正态分布\nμ={mu:.2f}, σ={std:.2f}')
axes[1].set_title('残差分布直方图')
axes[1].legend()

plt.tight_layout()
plt.show()

# 正式检验：Shapiro-Wilk检验
shapiro_test = stats.shapiro(residuals)
print(f"\nShapiro-Wilk正态性检验:")
print(f"  统计量: {shapiro_test[0]:.4f}, p值: {shapiro_test[1]:.4f}")
if shapiro_test[1] > 0.05:
    print("  结论：不能拒绝正态性原假设（残差近似正态分布）")
else:
    print("  警告：残差不服从正态分布")

# 3.4 多重共线性诊断
print("\n6.4 多重共线性诊断...")
# 计算VIF
vif_data = pd.DataFrame()
vif_data["变量"] = X_selected.columns
vif_data["VIF"] = [variance_inflation_factor(X_selected.values, i)
                   for i in range(X_selected.shape[1])]
vif_data["1/VIF"] = 1 / vif_data["VIF"]

print("方差膨胀因子(VIF):")
print(vif_data.sort_values("VIF", ascending=False))

# 可视化VIF
plt.figure(figsize=(10, 6))
bars = plt.barh(vif_data["变量"], vif_data["VIF"], color='skyblue')
plt.axvline(x=5, color='orange', linestyle='--', label='VIF=5 (轻度共线性)')
plt.axvline(x=10, color='red', linestyle='--', label='VIF=10 (严重共线性)')
plt.xlabel('方差膨胀因子 (VIF)')
plt.title('多重共线性诊断')
plt.legend()

# 在柱子上添加数值
for bar in bars:
    width = bar.get_width()
    plt.text(width + 0.1, bar.get_y() + bar.get_height()/2,
             f'{width:.2f}', ha='left', va='center')

plt.tight_layout()
plt.show()

# 3.5 自相关检验（如果数据有时间顺序）
print("\n6.5 自相关检验...")
# Durbin-Watson检验
dw_test = sm.stats.stattools.durbin_watson(residuals)
print(f"Durbin-Watson统计量: {dw_test:.4f}")
print("解释：")
print("  2.0：无自相关")
print("  0-2：正自相关")
print("  2-4：负自相关")
if 1.5 < dw_test < 2.5:
    print("  结论：无明显自相关问题")
else:
    print("  警告：可能存在自相关问题")

# 自相关图
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
fig, axes = plt.subplots(1, 2, figsize=(14, 4))
plot_acf(residuals, lags=20, ax=axes[0])
axes[0].set_title('残差自相关图')
plot_pacf(residuals, lags=20, ax=axes[1])
axes[0].set_title('残差偏自相关图')
plt.tight_layout()
plt.show()

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103

print("\n" + "="*60)
print("7. 模型修正与最终检验")
print("="*60)

# 4.1 处理异方差（如果存在）- 使用加权最小二乘法
if bp_test[1] < 0.05:  # 如果存在异方差
    print("\n检测到异方差，尝试使用加权最小二乘法(WLS)...")

    # 计算权重（通常使用残差绝对值的倒数）
    weights = 1 / (np.abs(residuals) + 0.0001)  # 加小数避免除零

    # 使用WLS重新拟合
    wls_model = sm.WLS(y, X_selected, weights=weights).fit()
    print("WLS模型摘要:")
    print(wls_model.summary())

    # 比较OLS和WLS
    print("\n模型比较:")
    print(f"OLS R-squared: {reduced_model.rsquared:.4f}")
    print(f"WLS R-squared: {wls_model.rsquared:.4f}")

    # 使用WLS模型的残差检验异方差
    wls_residuals = wls_model.resid
    wls_bp_test = het_breuschpagan(wls_residuals, wls_model.model.exog)
    print(f"\nWLS模型的Breusch-Pagan检验p值: {wls_bp_test[1]:.4f}")

    final_model = wls_model  # 使用WLS作为最终模型
else:
    print("\n无需修正异方差，使用OLS模型")
    final_model = reduced_model

# 4.2 模型显著性检验
print("\n7.1 模型整体显著性检验:")
f_test = final_model.fvalue
f_pvalue = final_model.f_pvalue
print(f"F统计量: {f_test:.4f}, p值: {f_pvalue:.4f}")
if f_pvalue < 0.05:
    print("结论：模型整体显著（至少有一个自变量对因变量有显著影响）")
else:
    print("警告：模型整体不显著")

# 4.3 回归系数显著性检验
print("\n7.2 回归系数显著性检验:")
coef_summary = pd.DataFrame({
    '变量': final_model.params.index,
    '系数': final_model.params.values,
    '标准误': final_model.bse.values,
    't值': final_model.tvalues.values,
    'p值': final_model.pvalues.values,
    '95%置信下限': final_model.conf_int()[0].values,
    '95%置信上限': final_model.conf_int()[1].values
})
coef_summary['显著'] = coef_summary['p值'] < 0.05
print(coef_summary.to_string(index=False))

# 4.4 模型评估指标
print("\n7.3 模型评估指标:")

# 训练集评估
y_pred_final = final_model.fittedvalues
mse = np.mean((y - y_pred_final) ** 2)
rmse = np.sqrt(mse)
mae = np.mean(np.abs(y - y_pred_final))
mape = np.mean(np.abs((y - y_pred_final) / y)) * 100

# 使用交叉验证评估
from sklearn.model_selection import cross_val_score
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, make_scorer

# 准备数据
X_cv = X_selected.drop('const', axis=1) if 'const' in X_selected.columns else X_selected
y_cv = y

# 交叉验证
cv_scores_rmse = cross_val_score(
    LinearRegression(), X_cv, y_cv,
    cv=5, scoring=make_scorer(lambda y_true, y_pred: np.sqrt(mean_squared_error(y_true, y_pred)))
)

print(f"训练集评估:")
print(f"  R-squared: {final_model.rsquared:.4f}")
print(f"  调整R-squared: {final_model.rsquared_adj:.4f}")
print(f"  MSE: {mse:,.2f}")
print(f"  RMSE: {rmse:,.2f}")
print(f"  MAE: {mae:,.2f}")
print(f"  MAPE: {mape:.2f}%")
print(f"\n交叉验证评估 (5折):")
print(f"  RMSE均值: {cv_scores_rmse.mean():,.2f}")
print(f"  RMSE标准差: {cv_scores_rmse.std():,.2f}")

# 4.5 残差诊断总结
print("\n7.4 模型诊断总结:")
diagnosis_summary = {
    '线性假设': '通过' if all(abs(residuals.corr(X_selected[col])) < 0.3 for col in X_selected.columns if col != 'const') else '未通过',
    '独立性假设': '通过' if 1.5 < dw_test < 2.5 else '未通过',
    '同方差假设': '通过' if bp_test[1] > 0.05 else '未通过',
    '正态性假设': '通过' if shapiro_test[1] > 0.05 else '未通过',
    '多重共线性': '通过' if all(vif_data["VIF"] < 10) else '未通过'
}

for assumption, status in diagnosis_summary.items():
    print(f"  {assumption}: {status}")

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168

print("\n" + "="*60)
print("8. 结论与应用")
print("="*60)

# 5.1 模型解释
print("\n8.1 模型解释:")
print("最终模型公式:")
model_formula = f"房价 = {final_model.params['const']:.2f}"
for var in final_model.params.index:
    if var != 'const':
        coef = final_model.params[var]
        if coef >= 0:
            model_formula += f" + {coef:.2f}×{var}"
        else:
            model_formula += f" - {abs(coef):.2f}×{var}"
print(model_formula)

print("\n显著自变量解释（在控制其他变量的情况下）:")
significant_vars = coef_summary[coef_summary['显著'] & (coef_summary['变量'] != 'const')]
for _, row in significant_vars.iterrows():
    var = row['变量']
    coef = row['系数']
    interpretation = f"  {var}: "
    if coef > 0:
        interpretation += f"每增加1单位，房价平均上涨${abs(coef):.2f}"
    else:
        interpretation += f"每增加1单位，房价平均下降${abs(coef):.2f}"
    print(interpretation)

# 5.2 模型可视化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 实际值 vs 预测值
axes[0, 0].scatter(y, y_pred_final, alpha=0.6)
axes[0, 0].plot([y.min(), y.max()], [y.min(), y.max()], 'r--', lw=2)
axes[0, 0].set_xlabel('实际房价 ($)')
axes[0, 0].set_ylabel('预测房价 ($)')
axes[0, 0].set_title('实际值 vs 预测值')
axes[0, 0].grid(True, alpha=0.3)

# 残差分布（最终模型）
axes[0, 1].hist(final_model.resid, bins=30, edgecolor='black', alpha=0.7)
axes[0, 1].axvline(x=0, color='red', linestyle='--')
axes[0, 1].set_xlabel('残差 ($)')
axes[0, 1].set_ylabel('频数')
axes[0, 1].set_title('最终模型残差分布')

# 特征重要性（标准化系数）
if 'const' in final_model.params.index:
    params_no_const = final_model.params.drop('const')
else:
    params_no_const = final_model.params

# 标准化系数
if len(X_cv.columns) > 0:
    std_coef = params_no_const * X_cv.std()
    std_coef = std_coef.abs().sort_values(ascending=True)
    axes[1, 0].barh(range(len(std_coef)), std_coef.values)
    axes[1, 0].set_yticks(range(len(std_coef)))
    axes[1, 0].set_yticklabels(std_coef.index)
    axes[1, 0].set_xlabel('标准化系数绝对值')
    axes[1, 0].set_title('变量重要性（标准化系数）')

# 预测误差分解
error = y - y_pred_final
abs_error = np.abs(error)
error_df = pd.DataFrame({
    '实际值': y,
    '预测值': y_pred_final,
    '误差': error,
    '绝对误差': abs_error
})

# 按实际值分组查看误差
if len(error_df) > 100:
    error_df['房价分组'] = pd.qcut(error_df['实际值'], 5)
    group_error = error_df.groupby('房价分组')['绝对误差'].mean()
    axes[1, 1].bar(range(len(group_error)), group_error.values)
    axes[1, 1].set_xticks(range(len(group_error)))
    axes[1, 1].set_xticklabels([str(g) for g in group_error.index], rotation=45, ha='right')
    axes[1, 1].set_xlabel('房价分组')
    axes[1, 1].set_ylabel('平均绝对误差 ($)')
    axes[1, 1].set_title('不同房价区间的预测误差')

plt.tight_layout()
plt.show()

# 5.3 模型使用建议
print("\n8.2 模型使用建议:")
print("1. 适用场景:")
print("   - 预测加州地区房屋的中位价格")
print("   - 分析不同因素对房价的影响程度")
print("   - 支持房地产投资决策")

print("\n2. 使用限制:")
print("   - 仅适用于加州地区，不能推广到其他地区")
print("   - 基于历史数据，市场变化可能影响预测准确性")
print("   - 模型假设条件需定期验证")

print("\n3. 注意事项:")
print("   - 用于预测时，输入变量应在训练数据范围内")
print("   - 定期用新数据重新训练模型")
print("   - 对于异常值较多的区域，预测准确性可能下降")

# 5.4 保存模型
print("\n8.3 模型保存:")
import joblib
import json

# 保存模型对象
model_info = {
    'model_type': '线性回归',
    'variables': list(X_selected.columns),
    'coefficients': {var: float(coef) for var, coef in final_model.params.items()},
    'r_squared': float(final_model.rsquared),
    'adj_r_squared': float(final_model.rsquared_adj),
    'rmse': float(rmse),
    'diagnostics': diagnosis_summary,
    'train_date': pd.Timestamp.now().strftime('%Y-%m-%d'),
    'data_source': '加州房价数据集'
}

# 保存为JSON文件
with open('housing_price_model.json', 'w') as f:
    json.dump(model_info, f, indent=2)

print("模型已保存为 'housing_price_model.json'")

# 5.5 模型预测示例
print("\n8.4 预测示例:")
print("输入示例数据:")

# 创建示例数据
example_data = {}
for var in X_selected.columns:
    if var != 'const':
        example_data[var] = X_selected[var].median()

example_df = pd.DataFrame([example_data])
print(example_df)

# 预测
if 'const' in final_model.params.index:
    prediction = final_model.params['const']
    for var, value in example_data.items():
        if var in final_model.params:
            prediction += final_model.params[var] * value
else:
    prediction = 0
    for var, value in example_data.items():
        if var in final_model.params:
            prediction += final_model.params[var] * value

print(f"\n预测房价: ${prediction:,.2f}")
print(f"预测区间（95%置信度）: ${prediction - 1.96*rmse:,.2f} - ${prediction + 1.96*rmse:,.2f}")

# 总结
print("\n" + "="*60)
print("建模过程总结")
print("="*60)
print("1. 数据探索: 识别了变量关系和异常值")
print("2. 变量选择: 使用逐步回归选择了显著变量")
print("3. 假设检验: 系统验证了所有回归假设")
print("4. 模型修正: 根据诊断结果调整模型")
print("5. 模型评估: 交叉验证确保模型泛化能力")
print("6. 模型解释: 明确系数的实际意义")
print("7. 模型部署: 提供使用建议和预测示例")
print("\n建模完成！")