Week 8: Correlation, Linear and Logistic Regression

import numpy as np
import pandas as pd
from scipy import stats
import matplotlib.pyplot as plt

print("="*50)
print("CORRELATION ANALYSIS")
print("="*50)

print("""
WHAT IS CORRELATION?
- Measures LINEAR relationship between two variables
- Ranges from -1 to +1
- Does NOT imply causation

CORRELATION INTERPRETATION:
- r = 1: Perfect positive relationship
- r = 0.5 to 0.8: Strong positive
- r = 0.3 to 0.5: Moderate positive
- r = 0 to 0.3: Weak positive
- r = 0: No relationship
- r = -0.3 to 0: Weak negative
- r = -0.5 to -0.3: Moderate negative
- r = -0.8 to -0.5: Strong negative
- r = -1: Perfect negative relationship
""")

# Create sample data
np.random.seed(42)
x = np.linspace(0, 10, 50)
y_positive = 2*x + np.random.normal(0, 3, 50)
y_negative = -1.5*x + np.random.normal(50, 3, 50)
y_no_corr = np.random.normal(20, 5, 50)

# Calculate correlations
corr_pos = np.corrcoef(x, y_positive)[0, 1]
corr_neg = np.corrcoef(x, y_negative)[0, 1]
corr_none = np.corrcoef(x, y_no_corr)[0, 1]

print(f"\nExample Correlations:")
print(f"Positive relationship: r = {corr_pos:.3f}")
print(f"Negative relationship: r = {corr_neg:.3f}")
print(f"No relationship: r = {corr_none:.3f}")

# Pearson Correlation
print("\n" + "="*50)
print("PEARSON CORRELATION COEFFICIENT")
print("="*50)

print("""
WHEN TO USE:
- Two continuous variables
- Linear relationship
- Data approximately normally distributed
- No extreme outliers

FORMULA:
r = Σ[(x - x̄)(y - ȳ)] / √[Σ(x - x̄)² × Σ(y - ȳ)²]
""")

# Real data example
data = pd.DataFrame({
'Age': [25, 30, 35, 40, 45, 50, 55, 60],
'Income': [35000, 45000, 55000, 65000, 75000, 85000, 95000, 105000],
'Experience': [2, 5, 8, 12, 15, 18, 20, 22]
})

print(f"\nData:")
print(data)

# Pearson correlation
pearson_r, p_value = stats.pearsonr(data['Age'], data['Income'])
print(f"\nPearson correlation (Age vs Income):")
print(f"r = {pearson_r:.4f}")
print(f"p-value = {p_value:.6f}")

if p_value < 0.05:
print("Correlation is STATISTICALLY SIGNIFICANT")
else:
print("Correlation is NOT statistically significant")

# Spearman Correlation
print("\n" + "="*50)
print("SPEARMAN CORRELATION COEFFICIENT")
print("="*50)

print("""
WHEN TO USE:
- Ordinal (ranked) data
- Non-linear monotonic relationships
- Data may not be normal
- More robust to outliers

Uses ranks instead of actual values
""")

# Example with ordinal data
rankings = pd.DataFrame({
'Movie_Rank': [1, 2, 3, 4, 5, 6, 7, 8],
'Audience_Score': [1, 3, 2, 5, 4, 8, 6, 7] # Different ranking
})

spearman_rho, p_spearman = stats.spearmanr(rankings['Movie_Rank'], rankings['Audience_Score'])
print(f"\nSpearman correlation (Movie rank vs audience score):")
print(f"ρ = {spearman_rho:.4f}")
print(f"p-value = {p_spearman:.6f}")

# Kendall Correlation
print("\n" + "="*50)
print("KENDALL CORRELATION COEFFICIENT")
print("="*50)

print("""
WHEN TO USE:
- Very small sample sizes
- Ordinal data
- Similar robustness to Spearman but different interpretation
""")

kendall_tau, p_kendall = stats.kendalltau(rankings['Movie_Rank'], rankings['Audience_Score'])
print(f"\nKendall correlation (Movie rank vs audience score):")
print(f"τ = {kendall_tau:.4f}")
print(f"p-value = {p_kendall:.6f}")

# Comparison
print("\n" + "="*50)
print("CORRELATION METHODS COMPARISON")
print("="*50)

print(f"Pearson r: {pearson_r:.4f} (parametric, assumes linearity)")
print(f"Spearman ρ: {spearman_rho:.4f} (non-parametric, ranks)")
print(f"Kendall τ: {kendall_tau:.4f} (non-parametric, ordinal associations)")

# Partial Correlation
print("\n" + "="*50)
print("PARTIAL CORRELATION")
print("="*50)

print("""
WHEN TO USE:
- Controlling for confounding variables
- Finding relationship between X and Y after removing effect of Z
- Understanding unique contribution of each variable
""")

# Create data with confounding
np.random.seed(42)
n = 100
Z = np.random.normal(100, 15, n) # Confounding variable
X = 0.5*Z + np.random.normal(0, 5, n)
Y = 0.6*Z + np.random.normal(0, 5, n)

# Simple correlations
r_xy = np.corrcoef(X, Y)[0, 1]
r_xz = np.corrcoef(X, Z)[0, 1]
r_yz = np.corrcoef(Y, Z)[0, 1]

print(f"\nCorrelation X-Y: {r_xy:.4f}")
print(f"Correlation X-Z: {r_xz:.4f}")
print(f"Correlation Y-Z: {r_yz:.4f}")

# Partial correlation: correlation of X and Y controlling for Z
# Using residuals approach
residuals_x = X - (r_xz * Z / np.std(Z)) * np.std(X)
residuals_y = Y - (r_yz * Z / np.std(Z)) * np.std(Y)
partial_r = np.corrcoef(residuals_x, residuals_y)[0, 1]

print(f"\nPartial correlation (X-Y | Z): {partial_r:.4f}")
print("The confounding Z variable reduced apparent correlation!")

# Visualization
fig, axes = plt.subplots(2, 3, figsize=(15, 10))

# Scatter plots with different correlations
datasets = [
(np.random.randn(100), np.random.randn(100), "No correlation (r≈0)"),
(x, y_positive, "Positive correlation (r≈0.87)"),
(x, y_negative, "Negative correlation (r≈-0.89)"),
(np.random.randn(100), np.random.randn(100), "Weak (r≈0)"),
(np.random.randn(100)*10, np.cumsum(np.random.randn(100)), "Non-linear"),
(data['Age'], data['Income']/1000, "Age vs Income (r≈1.00)")
]

for idx, (x_data, y_data, title) in enumerate(datasets):
row = idx // 3
col = idx % 3

r = np.corrcoef(x_data, y_data)[0, 1]
axes[row, col].scatter(x_data, y_data, alpha=0.6, s=50)
axes[row, col].set_title(f"{title}\nr = {r:.3f}")
axes[row, col].grid(True, alpha=0.3)

# Add trend line if sufficient correlation
if abs(r) > 0.3:
z = np.polyfit(x_data, y_data, 1)
p = np.poly1d(z)
x_line = np.linspace(x_data.min(), x_data.max(), 100)
axes[row, col].plot(x_line, p(x_line), "r--", alpha=0.8)

plt.suptitle('Correlation Examples', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

print("="*50)
print("CORRELATION MATRIX ANALYSIS")
print("="*50)

# Create dataset with multiple variables
np.random.seed(42)
data = pd.DataFrame({
'Age': np.random.randint(20, 60, 100),
'Salary': np.random.randint(30000, 150000, 100),
'Experience': np.random.randint(0, 30, 100),
'Performance': np.random.randint(1, 5, 100),
'Education': np.random.randint(10, 16, 100)
})

# Add some correlations to make it realistic
data['Salary'] = data['Experience'] * 3000 + data['Age'] * 500 + np.random.randint(-20000, 20000, 100)
data['Performance'] = 4 - (data['Age'] / 20) + np.random.normal(0, 0.3, 100)

print("Data:")
print(data.head())

# Calculate correlation matrix
corr_matrix = data.corr()

print("\nCorrelation Matrix:")
print(corr_matrix)

# Find strongest correlations
print("\n" + "="*50)
print("STRONGEST CORRELATIONS")
print("="*50)

# Get upper triangle of correlation matrix
mask = np.triu(np.ones_like(corr_matrix, dtype=bool), k=1)
upper_corr = corr_matrix.where(mask)

# Flatten and sort
corr_pairs = []
for i in range(len(corr_matrix.columns)):
for j in range(i+1, len(corr_matrix.columns)):
corr_pairs.append({
'Var1': corr_matrix.columns[i],
'Var2': corr_matrix.columns[j],
'Correlation': corr_matrix.iloc[i, j]
})

corr_df = pd.DataFrame(corr_pairs).sort_values('Correlation', key=abs, ascending=False)
print("\nTop 5 correlations:")
print(corr_df.head())

# Visualization: Heatmap
plt.figure(figsize=(10, 8))
sns.heatmap(corr_matrix, annot=True, fmt='.2f', cmap='coolwarm', center=0,
square=True, linewidths=1, cbar_kws={"shrink": 0.8})
plt.title('Correlation Matrix Heatmap', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()

# Pairplot (all relationships)
print("\nGenerating pairplot...")
sns.pairplot(data)
plt.suptitle('Pairplot: All Relationships', fontsize=14, fontweight='bold', y=0.995)
plt.tight_layout()
plt.show()

# Multicollinearity check
print("\n" + "="*50)
print("MULTICOLLINEARITY CHECK")
print("="*50)

print("""
MULTICOLLINEARITY:
- When independent variables are highly correlated
- Causes instability in regression coefficients
- Rule of thumb: r > 0.7 suggests multicollinearity
- Can be detected with VIF (Variance Inflation Factor)
""")

# Check for high correlations (excluding diagonal)
high_corr = []
for i in range(len(corr_matrix.columns)):
for j in range(i+1, len(corr_matrix.columns)):
if abs(corr_matrix.iloc[i, j]) > 0.7:
high_corr.append((corr_matrix.columns[i], corr_matrix.columns[j],
corr_matrix.iloc[i, j]))

if high_corr:
print("\nHigh correlations (|r| > 0.7):")
for var1, var2, r in high_corr:
print(f" {var1} <-> {var2}: r = {r:.3f}")
else:
print("\nNo high correlations detected (multicollinearity not severe)")

# Calculate VIF manually
print("\n" + "-"*50)
print("VARIANCE INFLATION FACTOR (VIF)")
print("-"*50)

from sklearn.preprocessing import StandardScaler
from numpy.linalg import inv

X = data[['Experience', 'Age', 'Education']].values
X_scaled = StandardScaler().fit_transform(X)
corr_X = np.corrcoef(X_scaled.T)
vif_matrix = inv(corr_X)

print(f"\nVIF for each variable:")
for i, var in enumerate(['Experience', 'Age', 'Education']):
vif = vif_matrix[i, i]
print(f" {var}: VIF = {vif:.3f}")
interpretation = "PROBLEMATIC" if vif > 10 else "MODERATE CONCERN" if vif > 5 else "OK"
print(f" → {interpretation}")

import numpy as np
import pandas as pd
from scipy import stats
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
import matplotlib.pyplot as plt

print("="*50)
print("SIMPLE LINEAR REGRESSION")
print("="*50)

print("""
GOAL:
Predict Y from X using linear relationship: Y = β₀ + β₁X + ε

COMPONENTS:
- Y: Dependent variable (what we predict)
- X: Independent variable (predictor)
- β₀: Intercept (Y when X=0)
- β₁: Slope (change in Y per unit X)
- ε: Error/residual (unexplained variation)

ASSUMPTIONS:
1. Linearity: Relationship is linear
2. Independence: Observations are independent
3. Homoscedasticity: Constant variance of errors
4. Normality: Errors are normally distributed
""")

# Real example: House price prediction
np.random.seed(42)
n = 100
square_feet = np.random.uniform(1000, 5000, n)
price = 50000 + 200*square_feet + np.random.normal(0, 50000, n)

data = pd.DataFrame({
'Square_Feet': square_feet,
'Price': price
})

print(f"\nExample: House Price Prediction")
print(f"Sample size: {n}")
print(data.head())

# Fit regression using sklearn
model = LinearRegression()
X = data[['Square_Feet']].values
y = data['Price'].values
model.fit(X, y)

# Get coefficients
intercept = model.intercept_
slope = model.coef_[0]

print(f"\nRegression Equation:")
print(f"Price = {intercept:.2f} + {slope:.2f} × Square_Feet")

# Predictions
y_pred = model.predict(X)

# Model evaluation
r2 = r2_score(y, y_pred)
rmse = np.sqrt(mean_squared_error(y, y_pred))
mae = np.mean(np.abs(y - y_pred))

print(f"\nModel Performance:")
print(f"R² (coefficient of determination): {r2:.4f}")
print(f" → {r2*100:.1f}% of variance in price explained by square feet")
print(f"RMSE (root mean squared error): ${rmse:,.2f}")
print(f"MAE (mean absolute error): ${mae:,.2f}")

# Statistical significance
# Calculate t-statistics manually
residuals = y - y_pred
mse = np.sum(residuals**2) / (len(y) - 2)
se_slope = np.sqrt(mse / np.sum((X - X.mean())**2))
t_stat = slope / se_slope
p_value = 2 * (1 - stats.t.cdf(abs(t_stat), len(y) - 2))

print(f"\nStatistical Significance:")
print(f"Slope coefficient: {slope:.2f}")
print(f"Standard error: {se_slope:.4f}")
print(f"t-statistic: {t_stat:.4f}")
print(f"p-value: {p_value:.2e}")

if p_value < 0.05:
print("→ Slope is SIGNIFICANTLY different from zero")
else:
print("→ Slope is NOT significantly different from zero")

# Confidence interval for slope
t_crit = stats.t.ppf(0.975, len(y) - 2)
ci_lower = slope - t_crit * se_slope
ci_upper = slope + t_crit * se_slope
print(f"95% CI for slope: [{ci_lower:.2f}, {ci_upper:.2f}]")

# Predict for new values
new_sqft = np.array([[2500], [3500], [4000]])
new_pred = model.predict(new_sqft)

print(f"\nPredictions for new houses:")
for sqft, pred in zip(new_sqft.flatten(), new_pred):
print(f" {sqft:.0f} sq ft → ${pred:,.0f}")

# Visualization
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Scatter plot with regression line
axes[0].scatter(X, y, alpha=0.6, s=50, label='Actual data')
axes[0].plot(X, y_pred, 'r-', linewidth=2, label='Regression line')
axes[0].set_xlabel('Square Feet', fontsize=12)
axes[0].set_ylabel('Price ($)', fontsize=12)
axes[0].set_title(f'Linear Regression\nR² = {r2:.3f}', fontsize=12, fontweight='bold')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Residuals plot
axes[1].scatter(y_pred, residuals, alpha=0.6, s=50)
axes[1].axhline(y=0, color='r', linestyle='--', linewidth=2)
axes[1].set_xlabel('Predicted Values', fontsize=12)
axes[1].set_ylabel('Residuals', fontsize=12)
axes[1].set_title('Residual Plot', fontsize=12, fontweight='bold')
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score, mean_squared_error
import matplotlib.pyplot as plt

print("="*50)
print("MULTIPLE LINEAR REGRESSION")
print("="*50)

print("""
GOAL:
Predict Y from multiple X variables
Y = β₀ + β₁X₁ + β₂X₂ + β₃X₃ + ... + ε

ADVANTAGES:
- Use multiple predictors
- Control for confounding variables
- Improve predictions
- Understand individual variable effects

CHALLENGES:
- Multicollinearity: predictors correlated
- More parameters to estimate
- Risk of overfitting with too many variables
""")

# Example: Salary prediction with multiple variables
np.random.seed(42)
n = 100
data = pd.DataFrame({
'Age': np.random.randint(25, 60, n),
'Experience': np.random.randint(0, 35, n),
'Education_Years': np.random.randint(12, 20, n),
'Performance_Score': np.random.uniform(2, 5, n)
})

# Create salary with known relationships
data['Salary'] = (
20000 +
1500 * data['Age'] +
2000 * data['Experience'] +
3000 * data['Education_Years'] +
10000 * data['Performance_Score'] +
np.random.normal(0, 10000, n)
)

print(f"\nData (Salary Prediction):")
print(data.head())

# Fit multiple regression
X = data[['Age', 'Experience', 'Education_Years', 'Performance_Score']]
y = data['Salary']

model = LinearRegression()
model.fit(X, y)

# Coefficients
intercept = model.intercept_
coefs = model.coef_

print(f"\nRegression Equation:")
print(f"Salary = {intercept:,.0f}")
for var, coef in zip(X.columns, coefs):
print(f" + {coef:.2f} × {var}")

# Model performance
y_pred = model.predict(X)
r2 = r2_score(y, y_pred)
rmse = np.sqrt(mean_squared_error(y, y_pred))

print(f"\nModel Performance:")
print(f"R²: {r2:.4f} ({r2*100:.1f}% variance explained)")
print(f"RMSE: ${rmse:,.0f}")

# Feature importance (standardized coefficients)
print(f"\nFeature Importance (Standardized Coefficients):")
X_std = (X - X.mean()) / X.std()
model_std = LinearRegression()
model_std.fit(X_std, y)

for var, coef in zip(X.columns, model_std.coef_):
print(f" {var}: {abs(coef):,.0f} (absolute importance)")

# Prediction with new data
new_employee = pd.DataFrame({
'Age': [35],
'Experience': [10],
'Education_Years': [16],
'Performance_Score': [4.5]
})

predicted_salary = model.predict(new_employee)
print(f"\nPredicted salary for new employee: ${predicted_salary[0]:,.0f}")

# Adjusted R²
adj_r2 = 1 - (1 - r2) * (n - 1) / (n - X.shape[1] - 1)
print(f"\nAdjusted R²: {adj_r2:.4f}")
print(f"(Adjusted for number of predictors)")

# Visualization: Actual vs Predicted
fig, ax = plt.subplots(figsize=(10, 6))
ax.scatter(y, y_pred, alpha=0.6, s=50)
ax.plot([y.min(), y.max()], [y.min(), y.max()], 'r--', linewidth=2)
ax.set_xlabel('Actual Salary', fontsize=12)
ax.set_ylabel('Predicted Salary', fontsize=12)
ax.set_title(f'Multiple Regression: Actual vs Predicted\nR² = {r2:.3f}',
fontsize=12, fontweight='bold')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

import numpy as np
import pandas as pd
from scipy import stats
from sklearn.linear_model import LinearRegression
import matplotlib.pyplot as plt

print("="*50)
print("REGRESSION DIAGNOSTICS")
print("="*50)

print("""
REGRESSION ASSUMPTIONS:
1. Linearity: Y = f(X) is linear
2. Independence: Errors are independent
3. Homoscedasticity: Constant error variance
4. Normality: Errors are normally distributed
5. No multicollinearity: Predictors not highly correlated

CHECKING ASSUMPTIONS:
- Residual plots
- Q-Q plot for normality
- Breusch-Pagan test for heteroscedasticity
- Durbin-Watson test for autocorrelation
""")

# Create dataset
np.random.seed(42)
X = np.linspace(0, 10, 100)
y = 5 + 2*X + np.random.normal(0, 2, 100)

# Fit model
model = LinearRegression()
model.fit(X.reshape(-1, 1), y)
y_pred = model.predict(X.reshape(-1, 1))
residuals = y - y_pred

print(f"\nResidual Statistics:")
print(f"Mean: {residuals.mean():.6f} (should be ~0)")
print(f"Std Dev: {residuals.std():.4f}")

# Normality test (Shapiro-Wilk)
stat_shapiro, p_shapiro = stats.shapiro(residuals)
print(f"\nShapiro-Wilk Test for Normality:")
print(f"p-value: {p_shapiro:.4f}")
print(f"Residuals are {'NORMAL' if p_shapiro > 0.05 else 'NOT NORMAL'}")

# Homoscedasticity test (Breusch-Pagan)
from scipy.stats import f as f_dist
residuals_sq = residuals ** 2
X_const = np.column_stack([np.ones(len(X)), X])
model_aux = LinearRegression()
model_aux.fit(X.reshape(-1, 1), residuals_sq)
ss_res_aux = np.sum((residuals_sq - model_aux.predict(X.reshape(-1, 1)))**2)
ss_tot_aux = np.sum((residuals_sq - residuals_sq.mean())**2)
r2_aux = 1 - ss_res_aux / ss_tot_aux

bp_stat = len(X) * r2_aux
bp_pvalue = 1 - f_dist.cdf(bp_stat, 1, len(X) - 2)

print(f"\nBreusch-Pagan Test for Homoscedasticity:")
print(f"BP statistic: {bp_stat:.4f}")
print(f"p-value: {bp_pvalue:.4f}")
print(f"Variance is {'CONSTANT' if bp_pvalue > 0.05 else 'NOT CONSTANT'}")

# Diagnostic plots
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# Plot 1: Residuals vs Fitted
axes[0, 0].scatter(y_pred, residuals, alpha=0.6)
axes[0, 0].axhline(y=0, color='r', linestyle='--')
axes[0, 0].set_xlabel('Fitted Values')
axes[0, 0].set_ylabel('Residuals')
axes[0, 0].set_title('Residuals vs Fitted')
axes[0, 0].grid(True, alpha=0.3)

# Plot 2: Q-Q plot
stats.probplot(residuals, dist="norm", plot=axes[0, 1])
axes[0, 1].set_title('Normal Q-Q Plot')
axes[0, 1].grid(True, alpha=0.3)

# Plot 3: Scale-Location
standardized_residuals = residuals / residuals.std()
axes[1, 0].scatter(y_pred, np.sqrt(np.abs(standardized_residuals)), alpha=0.6)
axes[1, 0].set_xlabel('Fitted Values')
axes[1, 0].set_ylabel('√|Standardized Residuals|')
axes[1, 0].set_title('Scale-Location Plot')
axes[1, 0].grid(True, alpha=0.3)

# Plot 4: Residuals histogram
axes[1, 1].hist(residuals, bins=20, edgecolor='black', alpha=0.7)
axes[1, 1].set_xlabel('Residuals')
axes[1, 1].set_ylabel('Frequency')
axes[1, 1].set_title('Distribution of Residuals')
axes[1, 1].grid(True, alpha=0.3, axis='y')

plt.suptitle('Regression Diagnostic Plots', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()

# Detecting outliers
print("\n" + "="*50)
print("DETECTING OUTLIERS AND INFLUENTIAL POINTS")
print("="*50)

# Cook's distance
n = len(X)
p = 2 # number of parameters
mse = np.sum(residuals**2) / (n - p)
leverage = 1/n + (X - X.mean())**2 / np.sum((X - X.mean())**2)
cooks_d = (residuals**2 / (p * mse)) * (leverage / (1 - leverage))

print(f"\nCook's distance (measures influence of each point):")
outlier_threshold = 4 / (n - p)
outliers = np.where(cooks_d > outlier_threshold)[0]

if len(outliers) > 0:
print(f"Found {len(outliers)} influential points")
for idx in outliers:
print(f" Point {idx}: Cook's D = {cooks_d[idx]:.4f}")
else:
print("No extreme outliers detected")

import numpy as np
import pandas as pd
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score
from sklearn.metrics import confusion_matrix, roc_auc_score, roc_curve
import matplotlib.pyplot as plt

print("="*50)
print("LOGISTIC REGRESSION")
print("="*50)

print("""
WHEN TO USE:
- Binary classification problems (yes/no, success/failure)
- Predicting probabilities of belonging to class 1
- When you need interpretable model

LINEAR REGRESSION LIMITATIONS:
- Predicts continuous values (-∞ to +∞)
- Cannot directly give probabilities (0 to 1)
- Violates assumptions

LOGISTIC REGRESSION SOLUTION:
- Uses sigmoid function: P(Y=1) = 1 / (1 + e^(-z))
- Outputs probability (0 to 1)
- Uses threshold (usually 0.5) to classify

ODDS AND LOG-ODDS:
- Odds = P(Y=1) / P(Y=0)
- Log-odds = ln(Odds)
- Logistic regression models log-odds linearly
""")

# Example: Customer churn prediction
np.random.seed(42)
n = 200

data = pd.DataFrame({
'Months_Customer': np.random.randint(1, 60, n),
'Monthly_Bill': np.random.uniform(20, 150, n),
'Tech_Support': np.random.choice([0, 1], n),
})

# Create churn probability based on features
prob_churn = (
0.9 / (1 + np.exp(-(0.1*data['Months_Customer'] - 2))) +
0.3 * (data['Monthly_Bill'] / 150) -
0.2 * data['Tech_Support']
)
prob_churn = np.clip(prob_churn, 0, 1)
data['Churn'] = (np.random.rand(n) < prob_churn).astype(int)

print(f"\nExample: Customer Churn Prediction")
print(data.head())
print(f"\nChurn rate: {data['Churn'].mean():.1%}")

# Fit logistic regression
X = data[['Months_Customer', 'Monthly_Bill', 'Tech_Support']]
y = data['Churn']

model = LogisticRegression()
model.fit(X, y)

# Coefficients
print(f"\nModel Coefficients:")
print(f"Intercept: {model.intercept_[0]:.4f}")
for feature, coef in zip(X.columns, model.coef_[0]):
print(f" {feature}: {coef:.4f}")

# Predictions
y_pred_prob = model.predict_proba(X)[:, 1]
y_pred = model.predict(X)

# Evaluation metrics
accuracy = accuracy_score(y, y_pred)
precision = precision_score(y, y_pred)
recall = recall_score(y, y_pred)
f1 = f1_score(y, y_pred)
auc = roc_auc_score(y, y_pred_prob)

print(f"\nModel Performance:")
print(f"Accuracy: {accuracy:.4f}")
print(f"Precision: {precision:.4f}")
print(f"Recall: {recall:.4f}")
print(f"F1-Score: {f1:.4f}")
print(f"AUC-ROC: {auc:.4f}")

# Confusion matrix
cm = confusion_matrix(y, y_pred)
print(f"\nConfusion Matrix:")
print(cm)
print(f"True Negatives: {cm[0,0]}")
print(f"False Positives: {cm[0,1]}")
print(f"False Negatives: {cm[1,0]}")
print(f"True Positives: {cm[1,1]}")

# Predictions for new data
new_customer = pd.DataFrame({
'Months_Customer': [6],
'Monthly_Bill': [100],
'Tech_Support': [0]
})

prob = model.predict_proba(new_customer)[0, 1]
prediction = model.predict(new_customer)[0]

print(f"\nPrediction for new customer:")
print(f" Probability of churn: {prob:.2%}")
print(f" Prediction: {'CHURN' if prediction == 1 else 'RETAIN'}")

# ROC Curve
fpr, tpr, thresholds = roc_curve(y, y_pred_prob)

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# ROC Curve
axes[0].plot(fpr, tpr, linewidth=2, label=f'AUC = {auc:.3f}')
axes[0].plot([0, 1], [0, 1], 'r--', label='Random classifier')
axes[0].set_xlabel('False Positive Rate')
axes[0].set_ylabel('True Positive Rate')
axes[0].set_title('ROC Curve')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Confusion Matrix Heatmap
import seaborn as sns
sns.heatmap(cm, annot=True, fmt='d', cmap='Blues', ax=axes[1],
xticklabels=['No Churn', 'Churn'],
yticklabels=['No Churn', 'Churn'])
axes[1].set_ylabel('Actual')
axes[1].set_xlabel('Predicted')
axes[1].set_title('Confusion Matrix')

plt.tight_layout()
plt.show()

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

print("="*50)
print("CORRELATION VS CAUSATION")
print("="*50)

print("""
CRITICAL DISTINCTION:

CORRELATION:
- Two variables move together
- Can be observed in data
- Does NOT imply causation

CAUSATION:
- Change in X causes change in Y
- Requires causal mechanism
- Cannot be determined from observational data alone

CONFOUNDING:
- Third variable affects both X and Y
- Creates spurious correlation
- Makes X appear to cause Y when it doesn't
""")

# Classic examples
examples = [
("Shoe size & Reading ability", "AGE (confound)"),
("Ice cream sales & Drowning deaths", "TEMPERATURE (confound)"),
("Number of firefighters & Damage from fire", "FIRE SEVERITY (confound)"),
]

print("\nClassic confounding examples:")
for corr, confound in examples:
print(f" Apparent: {corr}")
print(f" Real confounder: {confound}\n")

# Simulate confounding
print("="*50)
print("EXAMPLE: CONFOUNDING IN ACTION")
print("="*50)

np.random.seed(42)
n = 200

# Confounding variable
age = np.random.normal(40, 10, n)

# Exposure (influenced by confounder)
exercise = age * 0.5 + np.random.normal(0, 5, n)

# Outcome (influenced by confounder, NOT by exposure)
heart_disease = (age * 0.8 + np.random.normal(0, 10, n)) > 40

# But appears correlated
correlation = np.corrcoef(exercise, heart_disease.astype(int))[0, 1]

print(f"\nScenario: Does exercise prevent heart disease?")
print(f"Observed correlation: {correlation:.4f}")
print(f"Interpretation: Appears exercise {'prevents' if correlation < 0 else 'causes'} disease")

print(f"\nBUT the true story:")
print(f"- Age affects both exercise AND heart disease")
print(f"- Older people exercise less AND have more disease")
print(f"- It's AGE, not exercise!")

# Simpson's Paradox
print("\n" + "="*50)
print("SIMPSON'S PARADOX")
print("="*50)

print("""
A trend appears in groups separately but REVERSES when groups are combined!
""")

# Example: Treatment effectiveness
data = pd.DataFrame({
'Group': ['Mild', 'Mild', 'Severe', 'Severe', 'Overall', 'Overall'],
'Treatment': ['A', 'B', 'A', 'B', 'A', 'B'],
'Success_Rate': [0.95, 0.85, 0.65, 0.55, 0.80, 0.70]
})

print("\nTreatment success rates by severity:")
print(data.to_string(index=False))

print("\nWithin each group: Treatment A is better")
print("BUT: Treatment B looks better overall (because different sample sizes!)")

# Causal inference methods
print("\n" + "="*50)
print("ADDRESSING CONFOUNDING")
print("="*50)

print("""
RANDOMIZED EXPERIMENTS:
- Gold standard for causal inference
- Randomization balances confounders
- Can establish causation

OBSERVATIONAL DATA METHODS:
1. Matching: Compare similar subjects
2. Stratification: Analyze within strata
3. Regression adjustment: Control statistically
4. Instrumental variables: Find natural experiments
5. Propensity score: Match on probability of exposure
""")

# Visualization
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Scenario 1: Correlation with causation
np.random.seed(42)
x1 = np.random.randn(100)
y1 = 2*x1 + np.random.randn(100)*0.3
axes[0, 0].scatter(x1, y1, alpha=0.6)
axes[0, 0].set_title("Scenario 1: True Causation
X causes Y", fontweight='bold')
axes[0, 0].set_xlabel("X")
axes[0, 0].set_ylabel("Y")
axes[0, 0].grid(True, alpha=0.3)

# Scenario 2: Confounding
z = np.random.randn(100)
x2 = 0.5*z + np.random.randn(100)*0.3
y2 = 0.5*z + np.random.randn(100)*0.3
axes[0, 1].scatter(x2, y2, alpha=0.6)
axes[0, 1].set_title("Scenario 2: Confounding
Z causes both X and Y", fontweight='bold')
axes[0, 1].set_xlabel("X")
axes[0, 1].set_ylabel("Y")
axes[0, 1].grid(True, alpha=0.3)

# Scenario 3: No relationship
x3 = np.random.randn(100)
y3 = np.random.randn(100)
axes[1, 0].scatter(x3, y3, alpha=0.6)
axes[1, 0].set_title("Scenario 3: No Relationship", fontweight='bold')
axes[1, 0].set_xlabel("X")
axes[1, 0].set_ylabel("Y")
axes[1, 0].grid(True, alpha=0.3)

# Scenario 4: Reverse causation
y4 = np.random.randn(100)
x4 = 1.5*y4 + np.random.randn(100)*0.3
axes[1, 1].scatter(x4, y4, alpha=0.6)
axes[1, 1].set_title("Scenario 4: Reverse Causation
Y causes X", fontweight='bold')
axes[1, 1].set_xlabel("X")
axes[1, 1].set_ylabel("Y")
axes[1, 1].grid(True, alpha=0.3)

plt.suptitle('Causal Scenarios', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()

import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.linear_model import LinearRegression, Ridge, Lasso
from sklearn.metrics import mean_squared_error, r2_score
import matplotlib.pyplot as plt

print("="*50)
print("MODEL SELECTION AND VALIDATION")
print("="*50)

print("""
TRAIN-TEST SPLIT:
- Divide data: typically 70-80% train, 20-30% test
- Train model on training set
- Evaluate on test set (never seen by model)
- Prevents overfitting

CROSS-VALIDATION:
- k-fold: Divide into k chunks, train k times
- Leave-one-out: n-fold cross-validation
- Better use of data, more stable estimates

MODEL COMPLEXITY:
- Underfitting: Model too simple, high bias
- Overfitting: Model too complex, high variance
- Sweet spot: Good generalization
""")

# Create dataset
np.random.seed(42)
n = 200
X = np.linspace(0, 10, n).reshape(-1, 1)
y_true = np.sin(X.ravel())
y = y_true + np.random.normal(0, 0.15, n)

# Train-test split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

print(f"\nData split:")
print(f"Training set: {len(X_train)} samples")
print(f"Test set: {len(X_test)} samples")

# Test different polynomial degrees
degrees = [1, 2, 3, 5, 10, 20]
train_errors = []
test_errors = []

for degree in degrees:
from sklearn.preprocessing import PolynomialFeatures

# Create polynomial features
poly = PolynomialFeatures(degree=degree)
X_train_poly = poly.fit_transform(X_train)
X_test_poly = poly.transform(X_test)

# Train model
model = LinearRegression()
model.fit(X_train_poly, y_train)

# Evaluate
y_train_pred = model.predict(X_train_poly)
y_test_pred = model.predict(X_test_poly)

train_mse = mean_squared_error(y_train, y_train_pred)
test_mse = mean_squared_error(y_test, y_test_pred)

train_errors.append(train_mse)
test_errors.append(test_mse)

print(f"\nDegree {degree}:")
print(f" Train MSE: {train_mse:.6f}")
print(f" Test MSE: {test_mse:.6f}")
print(f" Status: {'UNDERFITTING' if degree < 3 else 'GOOD FIT' if degree == 3 else 'OVERFITTING'}")

# Regularization: Ridge vs Lasso
print("\n" + "="*50)
print("REGULARIZATION (Ridge and Lasso)")
print("="*50)

print("""
Regularization adds penalty for complex models:
- Ridge (L2): Reduces large coefficients
- Lasso (L1): Shrinks some coefficients to zero (feature selection)
- Elastic Net: Combination of Ridge and Lasso
""")

from sklearn.linear_model import Ridge, Lasso

# Create polynomial features
poly = PolynomialFeatures(degree=10)
X_train_poly = poly.fit_transform(X_train)
X_test_poly = poly.transform(X_test)

models = {
'Linear': LinearRegression(),
'Ridge (α=1)': Ridge(alpha=1),
'Ridge (α=10)': Ridge(alpha=10),
'Lasso (α=0.1)': Lasso(alpha=0.1)
}

for name, model in models.items():
model.fit(X_train_poly, y_train)
train_score = model.score(X_train_poly, y_train)
test_score = model.score(X_test_poly, y_test)

print(f"\n{name}:")
print(f" Train R²: {train_score:.4f}")
print(f" Test R²: {test_score:.4f}")

# Visualization
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Plot 1: Bias-Variance Tradeoff
axes[0].plot(degrees, train_errors, 'b-o', linewidth=2, label='Training error')
axes[0].plot(degrees, test_errors, 'r-o', linewidth=2, label='Test error')
axes[0].set_xlabel('Model Complexity (Polynomial Degree)')
axes[0].set_ylabel('Mean Squared Error')
axes[0].set_title('Bias-Variance Tradeoff', fontweight='bold')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Plot 2: Predictions for optimal model
poly3 = PolynomialFeatures(degree=3)
X_plot_poly = poly3.fit_transform(np.linspace(0, 10, 100).reshape(-1, 1))
model_opt = LinearRegression()
X_train_poly3 = poly3.fit_transform(X_train)
model_opt.fit(X_train_poly3, y_train)
y_plot_pred = model_opt.predict(X_plot_poly)

axes[1].scatter(X_train, y_train, alpha=0.5, label='Training data')
axes[1].scatter(X_test, y_test, alpha=0.5, label='Test data')
axes[1].plot(np.linspace(0, 10, 100), y_true, 'g-', linewidth=2, label='True function')
axes[1].plot(np.linspace(0, 10, 100), y_plot_pred, 'r--', linewidth=2, label='Degree 3 model')
axes[1].set_xlabel('X')
axes[1].set_ylabel('Y')
axes[1].set_title('Optimal Model Fit', fontweight='bold')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()