import numpy as np, pandas as pd, matplotlib.pyplot as plt
from scipy import stats
import seaborn as sns   # seaborn = high-level statistical plots (heatmaps, regplots, pairplots)
import statsmodels.api as sm
from statsmodels.formula.api import ols
from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.nonparametric.smoothers_lowess import lowess
from sklearn.linear_model import LinearRegression, Ridge, Lasso, LogisticRegression
from sklearn.preprocessing import StandardScaler, PolynomialFeatures
from sklearn.pipeline import make_pipeline
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.metrics import r2_score, mean_squared_error
ORG="#ea580c"; DEEP="#c2410c"; LIGHT="#fdba74"; INK="#1a2138"; GRID="#e6e9f2"; GREEN="#059669"; RED="#ef4444"; AMBER="#d97706"; BLUE="#2563eb"; PUR="#9333ea"; GREY="#94a3b8"; SLATE="#475569"
plt.rcParams.update({"figure.facecolor":"white","axes.facecolor":"white","figure.dpi":110,"font.size":11,
   "axes.edgecolor":GRID,"axes.grid":True,"grid.color":GRID,"axes.axisbelow":True,"axes.spines.top":False,
   "axes.spines.right":False,"axes.titlesize":12,"axes.titleweight":"bold","legend.frameon":False})
sns.set_style("whitegrid")
BASE="https://raw.githubusercontent.com/johnfisher-ai/Statistics-Data-Science-AI-Visual-Book/main/data/"

try:
    cars = pd.read_excel('../../data/multiple-linear-regression--cars.xlsx', sheet_name='Cars')
except FileNotFoundError:
    cars = pd.read_excel(BASE + 'multiple-linear-regression--cars.xlsx', sheet_name='Cars')

m = ols('mpg ~ weight + horsepower + C(origin)', cars).fit()
b = m.params['weight']
print(f'Holding horsepower and origin constant, each extra pound changes mpg by {b:.5f}')
print(f'-> about {b*500:.2f} mpg per 500 lbs.')

Holding horsepower and origin constant, each extra pound changes mpg by -0.00750
-> about -3.75 mpg per 500 lbs.

rng=np.random.default_rng(2); c2=cars.copy(); c2['noise']=rng.normal(size=len(cars))
a=ols('mpg ~ weight + C(origin)', cars).fit(); b=ols('mpg ~ weight + C(origin) + noise', c2).fit()
print(f'base   R2={a.rsquared:.4f} adjR2={a.rsquared_adj:.4f}')
print(f'+noise R2={b.rsquared:.4f} adjR2={b.rsquared_adj:.4f}  -> R2 up, adjR2 down')

base   R2=0.9020 adjR2=0.9009
+noise R2=0.9022 adjR2=0.9007  -> R2 up, adjR2 down

m1=ols('mpg ~ weight + C(origin)', cars).fit()
print('reference Europe:', {k:round(v,2) for k,v in m1.params.items() if 'origin' in k})
m2=ols('mpg ~ weight + C(origin, Treatment(reference="Japan"))', cars).fit()
print('reference Japan :', {k:round(v,2) for k,v in m2.params.items() if 'origin' in k})
print(f'Same fit either way: R2 {m1.rsquared:.4f} vs {m2.rsquared:.4f}')

reference Europe: {'C(origin)[T.Japan]': 1.41, 'C(origin)[T.USA]': -1.72}
reference Japan : {'C(origin, Treatment(reference="Japan"))[T.Europe]': -1.41, 'C(origin, Treatment(reference="Japan"))[T.USA]': -3.14}
Same fit either way: R2 0.9020 vs 0.9020

from statsmodels.stats.outliers_influence import variance_inflation_factor
Xv=sm.add_constant(cars[['weight','horsepower','displacement']].astype(float))
print({c:round(variance_inflation_factor(Xv.values,i),1) for i,c in enumerate(Xv.columns)})
full=ols('mpg ~ weight + horsepower + displacement + C(origin)', cars).fit()
red =ols('mpg ~ weight + C(origin)', cars).fit()
print(f'full adjR2={full.rsquared_adj:.3f}  reduced adjR2={red.rsquared_adj:.3f} (drop displacement+hp)')

{'const': np.float64(18.9), 'weight': np.float64(8.3), 'horsepower': np.float64(6.6), 'displacement': np.float64(13.3)}
full adjR2=0.902  reduced adjR2=0.901 (drop displacement+hp)

inter=ols('mpg ~ weight * C(origin)', cars).fit()
print(inter.pvalues.filter(like='weight:').round(3).to_string())
print('p-values > 0.05 -> no significant interaction; parallel-line (main-effects) model is preferred.')

weight:C(origin)[T.Japan]    0.196
weight:C(origin)[T.USA]      0.731
p-values > 0.05 -> no significant interaction; parallel-line (main-effects) model is preferred.