# Summary statistics for our two regression variables
penguins[['flipper_length_mm', 'body_mass_g']].describe().round(2)

fig, axes = plt.subplots(1, 2, figsize=(12, 3.5))
for ax, col, label in zip(axes,
  ['flipper_length_mm', 'body_mass_g'],
  ['Flipper length (mm)', 'Body mass (g)']):
    for sp, grp in penguins.groupby('species'):
        ax.hist(grp[col], bins=18, alpha=0.55, label=sp, density=True)
    ax.set_xlabel(label); ax.set_ylabel('Density')
    ax.set_title(f'Distribution of {label}')
    ax.legend(fontsize=9)
plt.tight_layout(); plt.show()

sns.scatterplot(
  data=penguins, x='flipper_length_mm', y='body_mass_g',
  hue='species', alpha=0.6, s=40
)
plt.xlabel('Flipper length (mm)'); plt.ylabel('Body mass (g)')
plt.title('Scatter Plot: Body Mass vs. Flipper Length')
plt.show()

# Step 1: Specify and fit the model using the formula interface
model = smf.ols('body_mass_g ~ flipper_length_mm', data=penguins).fit()

# Step 2: Print the full regression summary
print(model.summary())

# Extract the coefficients table as a DataFrame
coef_table = model.summary2().tables[1]
print(coef_table.round(4))

# Intercept β₀
b0 = model.params['Intercept']
print(f"\nIntercept  β̂₀ = {b0:.2f} g")

# Slope β₁
b1 = model.params['flipper_length_mm']
print(f"Slope β̂₁ = {b1:.4f} g per mm")

# 95% confidence intervals
ci = model.conf_int()
print(f"\n95% CI for β̂₁: ({ci.loc['flipper_length_mm', 0]:.3f}, {ci.loc['flipper_length_mm', 1]:.3f})")

# Fitted values and residuals
y_hat = model.fittedvalues
resid = model.resid

# Model fit statistics
print(f"R² = {model.rsquared:.4f}")
print(f"Adj. R² = {model.rsquared_adj:.4f}")

# F-test
print(f"\nF-statistic = {model.fvalue:.2f}")
print(f"Prob(F) = {model.f_pvalue:.2e}")

# Residual std deviation (σ̂)
print(f"\nσ̂ (RMSE) = {np.sqrt(model.mse_resid):.2f} g")
print(f"n = {int(model.nobs)}, df_resid = {int(model.df_resid)}")

# Build a prediction grid
x_grid = pd.DataFrame({
    'flipper_length_mm': np.linspace(
        penguins['flipper_length_mm'].min(),
        penguins['flipper_length_mm'].max(), 200)
})

# Get predictions with confidence interval
pred_ci = model.get_prediction(x_grid).summary_frame(alpha=0.05)

fig, ax = plt.subplots(figsize=(10, 6))
ax.scatter(penguins['flipper_length_mm'], penguins['body_mass_g'],
           alpha=0.4, color='steelblue', s=35, label='Observed data', zorder=3)
ax.plot(x_grid['flipper_length_mm'], pred_ci['mean'],
        color='crimson', lw=2.5, label='OLS fit')
ax.fill_between(x_grid['flipper_length_mm'],
                pred_ci['mean_ci_lower'], pred_ci['mean_ci_upper'],
                alpha=0.25, color='crimson', label='95% CI for mean')
ax.set_xlabel('Flipper length (mm)'); ax.set_ylabel('Body mass (g)')
ax.set_title(f'SLR: Body Mass ~ Flipper Length   ($R^2$ = {model.rsquared:.3f})')
ax.legend(); plt.tight_layout(); plt.show()

pred_pi = model.get_prediction(x_grid).summary_frame(alpha=0.05)

fig, ax = plt.subplots(figsize=(10, 6))
ax.scatter(penguins['flipper_length_mm'], penguins['body_mass_g'],
           alpha=0.4, color='steelblue', s=35, label='Observed data', zorder=3)
ax.plot(x_grid['flipper_length_mm'], pred_pi['mean'],
        color='crimson', lw=2.5, label='OLS fit')
ax.fill_between(x_grid['flipper_length_mm'],
                pred_pi['mean_ci_lower'], pred_pi['mean_ci_upper'],
                alpha=0.35, color='crimson', label='95% CI for mean response')
ax.fill_between(x_grid['flipper_length_mm'],
                pred_pi['obs_ci_lower'], pred_pi['obs_ci_upper'],
                alpha=0.12, color='steelblue', label='95% Prediction interval')
ax.set_xlabel('Flipper length (mm)'); ax.set_ylabel('Body mass (g)')
ax.set_title('Confidence Band vs. Prediction Band')
ax.legend(); plt.tight_layout(); plt.show()

# A new penguin with flipper length 195 mm
new_penguin = pd.DataFrame({'flipper_length_mm': [195]})

# Point prediction
y_hat_new = model.predict(new_penguin)
print(f"Predicted mass: {y_hat_new.values[0]:.1f} g")

# Full prediction with intervals
pred_new = (
  model.get_prediction(new_penguin)
  .summary_frame(alpha=0.05)
  )
print(f"\n95% CI for mean:  ({pred_new['mean_ci_lower'].values[0]:.1f},"
      f" {pred_new['mean_ci_upper'].values[0]:.1f}) g")
print(f"95% PI for obs:   ({pred_new['obs_ci_lower'].values[0]:.1f},"
      f" {pred_new['obs_ci_upper'].values[0]:.1f}) g")

# Extract fitted values and residuals
fitted   = model.fittedvalues
residual = model.resid

fig, ax = plt.subplots(figsize=(9, 5))
ax.scatter(fitted, residual, alpha=0.5, color='steelblue', s=35, zorder=3)
ax.axhline(0, color='crimson', lw=1.5, linestyle='--')
lowess = sm.nonparametric.lowess(residual, fitted, frac=0.5)
ax.plot(lowess[:, 0], lowess[:, 1], color='crimson', lw=2, label='LOWESS smoother')
ax.set_xlabel('Fitted values $\\hat{y}$'); ax.set_ylabel('Residuals $e$')
ax.set_title('Residuals vs. Fitted'); ax.legend()
plt.tight_layout(); plt.show()

from scipy.stats import probplot

qq_result = probplot(residual, dist='norm')

fig, ax = plt.subplots(figsize=(7, 6))
ax.scatter(qq_result[0][0], qq_result[0][1],
           alpha=0.55, color='steelblue', s=30, zorder=3)
ax.plot(qq_result[0][0],
        qq_result[1][1] + qq_result[1][0] * qq_result[0][0],
        color='crimson', lw=2, label='Reference line')
ax.set_xlabel('Theoretical quantiles'); ax.set_ylabel('Sample quantiles')
ax.set_title('Normal Q-Q Plot of Residuals'); ax.legend()
plt.tight_layout(); plt.show()

# Standardised residuals
resid_std = residual / np.sqrt(np.sum(residual**2) / (len(residual) - 2))
sqrt_abs  = np.sqrt(np.abs(resid_std))

fig, ax = plt.subplots(figsize=(9, 5))
ax.scatter(fitted, sqrt_abs, alpha=0.5, color='steelblue', s=35, zorder=3)
lowess = sm.nonparametric.lowess(sqrt_abs, fitted, frac=0.5)
ax.plot(lowess[:, 0], lowess[:, 1], color='crimson', lw=2, label='LOWESS smoother')
ax.set_xlabel('Fitted values $\\hat{y}$')
ax.set_ylabel('$\\sqrt{|\\text{Standardised residuals}|}$')
ax.set_title('Scale-Location Plot'); ax.legend()
plt.tight_layout(); plt.show()

# Compute influence measures via statsmodels
influence = model.get_influence()
cooks_d   = influence.cooks_distance[0]
threshold = 4 / len(penguins)

fig, ax = plt.subplots(figsize=(11, 4))
ax.stem(range(len(cooks_d)), cooks_d,
        linefmt='steelblue', markerfmt='o', basefmt=' ')
ax.axhline(threshold, color='crimson', linestyle='--',
           label=f'Threshold 4/n = {threshold:.4f}')

# Flag influential points
flagged = np.where(cooks_d > threshold)[0]
print(f"Observations above threshold: {len(flagged)}")
ax.set_xlabel('Observation index'); ax.set_ylabel("Cook's distance $D_i$")
ax.set_title("Cook's Distance"); ax.legend()
plt.tight_layout(); plt.show()

fig, axes = plt.subplots(2, 2, figsize=(13, 6))

# 1. Residuals vs Fitted
axes[0,0].scatter(fitted, residual, alpha=0.45, color='steelblue', s=25, zorder=3)
axes[0,0].axhline(0, color='crimson', lw=1.5, linestyle='--')
lw1 = sm.nonparametric.lowess(residual, fitted, frac=0.5)
axes[0,0].plot(lw1[:,0], lw1[:,1], 'crimson', lw=2)
axes[0,0].set_xlabel('Fitted values'); axes[0,0].set_ylabel('Residuals')
axes[0,0].set_title('Residuals vs. Fitted')

# 2. Normal Q-Q
qq = probplot(residual, dist='norm')
axes[0,1].scatter(qq[0][0], qq[0][1], alpha=0.5, color='steelblue', s=25, zorder=3)
axes[0,1].plot(qq[0][0], qq[1][1] + qq[1][0]*qq[0][0], 'crimson', lw=2, label='Reference')
axes[0,1].set_xlabel('Theoretical quantiles'); axes[0,1].set_ylabel('Sample quantiles')
axes[0,1].set_title('Normal Q-Q'); axes[0,1].legend(fontsize=9)

# 3. Scale-Location
axes[1,0].scatter(fitted, sqrt_abs, alpha=0.45, color='steelblue', s=25, zorder=3)
lw2 = sm.nonparametric.lowess(sqrt_abs, fitted, frac=0.5)
axes[1,0].plot(lw2[:,0], lw2[:,1], 'crimson', lw=2)
axes[1,0].set_xlabel('Fitted values')
axes[1,0].set_ylabel('$\\sqrt{|\\mathrm{Std.\\ Residuals}|}$')
axes[1,0].set_title('Scale-Location')

# 4. Cook's Distance
axes[1,1].stem(range(len(cooks_d)), cooks_d,
               linefmt='steelblue', markerfmt='o', basefmt=' ')
axes[1,1].axhline(threshold, color='crimson', linestyle='--', label=f'4/n = {threshold:.3f}')
axes[1,1].set_xlabel('Observation index')
axes[1,1].set_ylabel("Cook's $D_i$")
axes[1,1].set_title("Cook's Distance"); axes[1,1].legend(fontsize=9)

plt.suptitle("Regression Diagnostic Panel", fontsize=14, y=1.01)
plt.tight_layout(); plt.show()

fig, axes = plt.subplots(2, 2, figsize=(12, 6), sharex=False, sharey=False)

for ax, (ds, grp) in zip(axes.flat, anscombe.groupby('dataset')):
    model_a = smf.ols('y ~ x', data=grp).fit()
    x_a     = np.linspace(grp['x'].min(), grp['x'].max(), 100)
    y_a     = model_a.params['Intercept'] + model_a.params['x'] * x_a

    ax.scatter(grp['x'], grp['y'], color='steelblue', s=60, zorder=3, alpha=0.8)
    ax.plot(x_a, y_a, color='crimson', lw=2)
    ax.set_title(f"Dataset {ds}  (R² = {model_a.rsquared:.3f})", fontsize=11)
    ax.set_xlabel('x'); ax.set_ylabel('y')

plt.suptitle("Anscombe's Quartet: Same Statistics, Different Structure", fontsize=13, y=1.01)
plt.tight_layout(); plt.show()