import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.formula.api import ols
from scipy import stats

# 数据
data = {
    'sales': [162, 120, 223, 131, 67, 169, 81, 192, 116, 55, 252, 232, 144, 103, 212],
    'population': [27.4, 18.0, 37.5, 20.5, 8.6, 26.5, 9.8, 33.0, 19.5, 5.3, 43.0, 37.2, 23.6, 15.7, 37.0],
    'income': [2450, 3254, 3802, 2838, 2347, 3782, 3008, 2450, 2137, 2560, 4020, 4427, 2660, 2088, 2605]
}

df = pd.DataFrame(data)

# 创建交互项
df['interaction'] = df['population'] * df['income']

# 模型1: 不包含交互项
model1 = ols('sales ~ population + income', data=df).fit()
sse1 = sum(model1.resid ** 2)
print("\n模型1(无交互项)摘要:")
print(model1.summary())

# 模型2: 包含交互项
model2 = ols('sales ~ population + income + interaction', data=df).fit()
sse2 = sum(model2.resid ** 2)
print("\n模型2(有交互项)摘要:")
print(model2.summary())

# 计算F统计量
df1 = 1  # 交互项增加的参数数量
df2 = 15 - 4  # n - p (样本量 - 扩展模型参数数量)
f_value = ((sse1 - sse2) / df1) / (sse2 / df2)
p_value = 1 - stats.f.cdf(f_value, df1, df2)

print("\n===== 交互项显著性检验 =====")
print(f"SSE1 (无交互项): {sse1:.4f}")
print(f"SSE2 (有交互项): {sse2:.4f}")
print(f"F值: {f_value:.4f}")
print(f"p值: {p_value:.4f}")
print(f"F临界值 (α=0.05): {stats.f.ppf(0.95, df1, df2):.4f}")

if p_value < 0.05:
    print("结论: 交互项显著 (p < 0.05), 应当保留在模型中")
else:
    print("结论: 交互项不显著 (p > 0.05), 可以从模型中移除")