# 导入必要的库
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import statsmodels.api as sm
from statsmodels.formula.api import ols
from statsmodels.stats.anova import anova_lm
import scipy.stats as stats
from sklearn.preprocessing import PolynomialFeatures
import warnings
warnings.filterwarnings('ignore')
# 设置中文字体 - 确保正确显示标题
sns.set(style="whitegrid", palette="muted")
plt.rcParams['font.sans-serif'] = ['SimHei'] # 解决中文显示问题
plt.rcParams['axes.unicode_minus'] = False # 解决负号显示问题
plt.rcParams['figure.dpi'] = 100 # 提高分辨率
# 创建数据框
data = {
'温度': [30, 30, 30, 30, 30, 30, 30, 30, 30, 40, 40, 40, 40, 40, 40, 40, 40, 40, 50, 50, 50, 50, 50, 50, 50, 50, 50],
'湿度': [50, 50, 50, 70, 70, 70, 90, 90, 90, 50, 50, 50, 70, 70, 70, 90, 90, 90, 50, 50, 50, 70, 70, 70, 90, 90, 90],
'固含量': [6, 6, 6, 8, 8, 8, 10, 10, 10, 8, 8, 8, 10, 10, 10, 6, 6, 6, 10, 10, 10, 6, 6, 6, 8, 8, 8],
'孔面积占比': [18.09, 17.04, 19.45, 9.11, 9.4, 8.85, 9.35, 9.4, 9.08,
30.33, 31.2, 29.2, 24.12, 24.34, 23.11, 36.01, 35.35, 35.66,
6.41, 6.62, 6.73, 8.38, 8.3, 8.29, 8.38, 8.3, 8.29]
}
df = pd.DataFrame(data)
# 添加随机误差项 - 解释同一条件下的变异
np.random.seed(42)
df['随机误差'] = np.random.normal(0, 0.5, len(df))
df['孔面积占比'] = df['孔面积占比'] + df['随机误差']
# 将分类变量转换为类别类型
df['温度'] = df['温度'].astype('category')
df['湿度'] = df['湿度'].astype('category')
df['固含量'] = df['固含量'].astype('category')
# 进行三因素ANOVA分析
model = ols('孔面积占比 ~ C(温度) + C(湿度) + C(固含量) + C(温度):C(湿度) + C(温度):C(固含量) + C(湿度):C(固含量) + C(温度):C(湿度):C(固含量)', data=df).fit()
anova_results = anova_lm(model, typ=2)
# 修复NaN值问题:正确计算效应量η²和偏η²
total_ss = anova_results['sum_sq'].sum()
residual_ss = anova_results.loc['Residual', 'sum_sq']
anova_results['η²'] = anova_results['sum_sq'] / total_ss
anova_results['偏η²'] = anova_results['sum_sq'] / (anova_results['sum_sq'] + residual_ss)
# 处理p=0的情况:科学计数法显示极小值
anova_results['PR(>F)'] = anova_results['PR(>F)'].apply(
lambda x: f"{x:.2e}" if x < 0.001 else f"{x:.4f}"
)
# 对于残差行,偏η²设为NaN
anova_results.at['Residual', '偏η²'] = np.nan
print("\nANOVA分析结果:")
print(anova_results)
# 设置图形风格 - 确保标题显示
sns.set(style="whitegrid", font_scale=1.2)
plt.figure(figsize=(15, 12))
# 1. 主效应图
plt.subplot(2, 2, 1)
sns.pointplot(x='温度', y='孔面积占比', data=df, ci=95)
plt.title('温度对孔面积占比的主效应', fontsize=14, pad=20) # 增加标题与图的间距
plt.subplot(2, 2, 2)
sns.pointplot(x='湿度', y='孔面积占比', data=df, ci=95)
plt.title('湿度对孔面积占比的主效应', fontsize=14, pad=20)
plt.subplot(2, 2, 3)
sns.pointplot(x='固含量', y='孔面积占比', data=df, ci=95)
plt.title('固含量对孔面积占比的主效应', fontsize=14, pad=20)
# 2. 交互作用图 - 温度与湿度
plt.subplot(2, 2, 4)
sns.lineplot(x='温度', y='孔面积占比', hue='湿度', data=df, ci=95, marker='o')
plt.title('温度与湿度的交互作用', fontsize=14, pad=20)
plt.tight_layout(pad=3.0) # 增加子图间距
plt.savefig('主效应与交互作用.png', bbox_inches='tight') # 保存图像确保显示
plt.show()
# 更多交互作用图
plt.figure(figsize=(15, 5))
plt.subplot(1, 3, 1)
sns.lineplot(x='温度', y='孔面积占比', hue='固含量', data=df, ci=95, marker='o')
plt.title('温度与固含量的交互作用', fontsize=14, pad=15)
plt.subplot(1, 3, 2)
sns.lineplot(x='湿度', y='孔面积占比', hue='固含量', data=df, ci=95, marker='o')
plt.title('湿度与固含量的交互作用', fontsize=14, pad=15)
plt.subplot(1, 3, 3)
heatmap_data = df.groupby(['温度', '湿度'])['孔面积占比'].mean().unstack()
sns.heatmap(heatmap_data, annot=True, fmt=".1f", cmap="YlOrRd")
plt.title('温度与湿度组合下的平均孔面积占比', fontsize=14, pad=15)
plt.tight_layout(pad=2.0)
plt.savefig('交互作用热图.png', bbox_inches='tight')
plt.show()
# 将变量转换为数值型以进行回归分析
df_numeric = df.copy()
df_numeric['温度'] = df_numeric['温度'].astype(int)
df_numeric['湿度'] = df_numeric['湿度'].astype(int)
df_numeric['固含量'] = df_numeric['固含量'].astype(int)
# 线性模型
linear_model = ols('孔面积占比 ~ 温度 + 湿度 + 固含量', data=df_numeric).fit()
linear_r2 = linear_model.rsquared
# 多项式模型(包含二次项) - 使用矩阵方法避免公式解析问题
poly = PolynomialFeatures(degree=2, include_bias=False)
X_poly = poly.fit_transform(df_numeric[['温度', '湿度', '固含量']])
poly_feature_names = poly.get_feature_names_out(['温度', '湿度', '固含量'])
# 修复特征名称:替换空格和特殊字符
poly_feature_names = [name.replace(' ', '_').replace('^', '_squared') for name in poly_feature_names]
# 创建多项式特征DataFrame
df_poly = pd.DataFrame(X_poly, columns=poly_feature_names)
y = df_numeric['孔面积占比']
# 添加截距项
X_poly_with_const = sm.add_constant(df_poly)
# 构建多项式模型
poly_model = sm.OLS(y, X_poly_with_const).fit()
poly_r2 = poly_model.rsquared
# 比较模型拟合
print(f"\n线性模型R²: {linear_r2:.4f}")
print(f"多项式模型R²: {poly_r2:.4f}")
# 残差分析 - 确保图表标题显示
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.scatter(linear_model.predict(), linear_model.resid)
plt.axhline(y=0, color='r', linestyle='--')
plt.xlabel('预测值', fontsize=12)
plt.ylabel('残差', fontsize=12)
plt.title('线性模型残差图', fontsize=14, pad=15)
plt.subplot(1, 2, 2)
plt.scatter(poly_model.predict(), poly_model.resid)
plt.axhline(y=0, color='r', linestyle='--')
plt.xlabel('预测值', fontsize=12)
plt.ylabel('残差', fontsize=12)
plt.title('多项式模型残差图', fontsize=14, pad=15)
plt.tight_layout(pad=2.0)
plt.savefig('残差分析.png', bbox_inches='tight')
plt.show()
# 最佳参数组合分析
grouped = df.groupby(['温度', '湿度', '固含量'])['孔面积占比'].mean().reset_index()
max_porosity = grouped.loc[grouped['孔面积占比'].idxmax()]
print("\n最佳参数组合:")
print(f"温度: {max_porosity['温度']}°C")
print(f"湿度: {max_porosity['湿度']}%")
print(f"固含量: {max_porosity['固含量']}%")
print(f"孔面积占比: {max_porosity['孔面积占比']:.2f}%")
# 绘制三维散点图 - 确保标题显示
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
scatter = ax.scatter(
df_numeric['温度'],
df_numeric['湿度'],
df_numeric['固含量'],
c=df_numeric['孔面积占比'],
cmap='viridis',
s=100
)
ax.set_xlabel('温度 (°C)', fontsize=12, labelpad=10)
ax.set_ylabel('湿度 (%)', fontsize=12, labelpad=10)
ax.set_zlabel('固含量 (%)', fontsize=12, labelpad=10)
ax.set_title('参数空间中的孔面积占比分布', fontsize=14, pad=20)
cbar = plt.colorbar(scatter, pad=0.15)
cbar.set_label('孔面积占比 (%)', fontsize=12)
plt.savefig('三维散点图.png', bbox_inches='tight')
plt.show()
# 输出ANOVA结果的详细解释
print("\nANOVA结果详细解释:")
print("="*50)
for factor in anova_results.index:
if factor != 'Residual':
p_value = anova_results.loc[factor, 'PR(>F)']
eta_sq = anova_results.loc[factor, 'η²']
significance = "高度显著" if isinstance(p_value, str) or p_value < 0.001 else "显著" if p_value < 0.05 else "不显著"
print(f"{factor}: p值 = {p_value}, η² = {eta_sq:.4f} ({significance})")
# 输出模型比较结果
print("\n模型比较:")
print("="*50)
print(f"线性模型R²: {linear_r2:.4f}")
print(f"多项式模型R²: {poly_r2:.4f}")
if poly_r2 > linear_r2:
improvement = ((poly_r2 - linear_r2) / linear_r2) * 100
print(f"多项式模型比线性模型拟合更好,改进幅度: {improvement:.2f}%")
else:
print("线性模型比多项式模型拟合更好")
# 检查残差的正态性
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
stats.probplot(linear_model.resid, plot=plt)
plt.title('线性模型残差Q-Q图', fontsize=14, pad=15)
plt.subplot(1, 2, 2)
stats.probplot(poly_model.resid, plot=plt)
plt.title('多项式模型残差Q-Q图', fontsize=14, pad=15)
plt.tight_layout(pad=2.0)
plt.savefig('残差QQ图.png', bbox_inches='tight')
plt.show()
# Shapiro-Wilk正态性检验
_, p_linear = stats.shapiro(linear_model.resid)
_, p_poly = stats.shapiro(poly_model.resid)
print(f"\n残差正态性检验 (Shapiro-Wilk):")
print(f"线性模型p值: {p_linear:.4f} {'(正态)' if p_linear > 0.05 else '(非正态)'}")
print(f"多项式模型p值: {p_poly:.4f} {'(正态)' if p_poly > 0.05 else '(非正态)'}")
# 添加随机误差分析
print("\n随机误差分析:")
print("="*50)
print(f"随机误差标准差: {df['随机误差'].std():.4f}")
print(f"随机误差对孔面积占比的贡献比例: {df['随机误差'].var() / df['孔面积占比'].var():.4f}")
# 输出多项式模型系数
print("\n多项式模型系数:")
print("="*50)
for i, coef in enumerate(poly_model.params):
if i == 0:
print(f"截距: {coef:.4f}")
else:
print(f"{poly_feature_names[i-1]}: {coef:.4f}")
# 解释p=0的结果
print("\n关于p=0结果的解释:")
print("="*50)
print("在ANOVA分析中,p值显示为0或极小值(如1.23e-15)表示:")
print("1. 统计显著性极高,几乎可以排除随机因素导致差异的可能性")
print("2. 在实际应用中应报告为p<0.001")
print("3. 但需注意:极小的p值可能源于过大的样本量或过小的随机误差")
print("4. 建议结合效应量(η²)评估实际意义:")
print(" - η² > 0.14: 大效应")
print(" - η² > 0.06: 中等效应")
print(" - η² > 0.01: 小效应")
为什么结果图上标题位置都是方框不显示结果
最新发布