{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a478b36b",
   "metadata": {},
   "source": [
    "# Week 8: Correlation, Linear and Logistic Regression\n",
    "\n",
    "*Kmex Consult — Data Science Fundamentals Course*\n",
    "\n",
    "Work through each code cell in order. Edit and re-run to experiment."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39702252",
   "metadata": {},
   "source": [
    "*Data Science Fundamentals Course*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5df558a2",
   "metadata": {},
   "source": [
    "## Week 8 Overview\n",
    "Week 8 focuses on understanding relationships between variables. These techniques are fundamental to predictive modeling and understanding data patterns. You'll learn how to:\n",
    "\n",
    "- Quantify relationships between variables\n",
    "- Build predictive models using regression\n",
    "- Understand the critical distinction between correlation and causation\n",
    "- Diagnose and fix regression problems\n",
    "- Extend regression to binary outcomes\n",
    "\n",
    "Regression is one of the most widely used techniques in data science and business:\n",
    "- Sales forecasting based on marketing spend\n",
    "- House price prediction from features\n",
    "- Customer lifetime value prediction\n",
    "- Risk assessment and credit scoring\n",
    "- Treatment effect estimation in experiments\n",
    "\n",
    "By the end of Week 8, you will be able to:\n",
    "- Calculate and interpret various correlation coefficients\n",
    "- Build simple and multiple linear regression models\n",
    "- Diagnose regression assumptions and address violations\n",
    "- Evaluate regression model performance\n",
    "- Understand polynomial and non-linear relationships\n",
    "- Build logistic regression for classification\n",
    "- Understand causality and confounding bias\n",
    "- Make causal inferences appropriately\n",
    "\n",
    "Week 8 is divided into three 2-hour sessions:\n",
    "- Session 1: Correlation Analysis and Relationships\n",
    "- Session 2: Linear Regression: Simple and Multiple\n",
    "- Session 3: Logistic Regression, Model Diagnostics, and Causality"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c11cc832",
   "metadata": {},
   "source": [
    "## SESSION 1: Correlation Analysis and Relationships"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28e68b24",
   "metadata": {},
   "source": [
    "### Duration: 2 hours"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e66d0fc",
   "metadata": {},
   "source": [
    "### 1.1 Understanding Correlation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37d0dcc0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy import stats\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "print(\"=\"*50)\n",
    "print(\"CORRELATION ANALYSIS\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "print(\"\"\"\n",
    "WHAT IS CORRELATION?\n",
    "- Measures LINEAR relationship between two variables\n",
    "- Ranges from -1 to +1\n",
    "- Does NOT imply causation\n",
    "\n",
    "CORRELATION INTERPRETATION:\n",
    "- r = 1: Perfect positive relationship\n",
    "- r = 0.5 to 0.8: Strong positive\n",
    "- r = 0.3 to 0.5: Moderate positive\n",
    "- r = 0 to 0.3: Weak positive\n",
    "- r = 0: No relationship\n",
    "- r = -0.3 to 0: Weak negative\n",
    "- r = -0.5 to -0.3: Moderate negative\n",
    "- r = -0.8 to -0.5: Strong negative\n",
    "- r = -1: Perfect negative relationship\n",
    "\"\"\")\n",
    "\n",
    "# Create sample data\n",
    "np.random.seed(42)\n",
    "x = np.linspace(0, 10, 50)\n",
    "y_positive = 2*x + np.random.normal(0, 3, 50)\n",
    "y_negative = -1.5*x + np.random.normal(50, 3, 50)\n",
    "y_no_corr = np.random.normal(20, 5, 50)\n",
    "\n",
    "# Calculate correlations\n",
    "corr_pos = np.corrcoef(x, y_positive)[0, 1]\n",
    "corr_neg = np.corrcoef(x, y_negative)[0, 1]\n",
    "corr_none = np.corrcoef(x, y_no_corr)[0, 1]\n",
    "\n",
    "print(f\"\\nExample Correlations:\")\n",
    "print(f\"Positive relationship: r = {corr_pos:.3f}\")\n",
    "print(f\"Negative relationship: r = {corr_neg:.3f}\")\n",
    "print(f\"No relationship: r = {corr_none:.3f}\")\n",
    "\n",
    "# Pearson Correlation\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"PEARSON CORRELATION COEFFICIENT\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "print(\"\"\"\n",
    "WHEN TO USE:\n",
    "- Two continuous variables\n",
    "- Linear relationship\n",
    "- Data approximately normally distributed\n",
    "- No extreme outliers\n",
    "\n",
    "FORMULA:\n",
    "r = Σ[(x - x̄)(y - ȳ)] / √[Σ(x - x̄)² × Σ(y - ȳ)²]\n",
    "\"\"\")\n",
    "\n",
    "# Real data example\n",
    "data = pd.DataFrame({\n",
    "'Age': [25, 30, 35, 40, 45, 50, 55, 60],\n",
    "'Income': [35000, 45000, 55000, 65000, 75000, 85000, 95000, 105000],\n",
    "'Experience': [2, 5, 8, 12, 15, 18, 20, 22]\n",
    "})\n",
    "\n",
    "print(f\"\\nData:\")\n",
    "print(data)\n",
    "\n",
    "# Pearson correlation\n",
    "pearson_r, p_value = stats.pearsonr(data['Age'], data['Income'])\n",
    "print(f\"\\nPearson correlation (Age vs Income):\")\n",
    "print(f\"r = {pearson_r:.4f}\")\n",
    "print(f\"p-value = {p_value:.6f}\")\n",
    "\n",
    "if p_value < 0.05:\n",
    "print(\"Correlation is STATISTICALLY SIGNIFICANT\")\n",
    "else:\n",
    "print(\"Correlation is NOT statistically significant\")\n",
    "\n",
    "# Spearman Correlation\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"SPEARMAN CORRELATION COEFFICIENT\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "print(\"\"\"\n",
    "WHEN TO USE:\n",
    "- Ordinal (ranked) data\n",
    "- Non-linear monotonic relationships\n",
    "- Data may not be normal\n",
    "- More robust to outliers\n",
    "\n",
    "Uses ranks instead of actual values\n",
    "\"\"\")\n",
    "\n",
    "# Example with ordinal data\n",
    "rankings = pd.DataFrame({\n",
    "'Movie_Rank': [1, 2, 3, 4, 5, 6, 7, 8],\n",
    "'Audience_Score': [1, 3, 2, 5, 4, 8, 6, 7] # Different ranking\n",
    "})\n",
    "\n",
    "spearman_rho, p_spearman = stats.spearmanr(rankings['Movie_Rank'], rankings['Audience_Score'])\n",
    "print(f\"\\nSpearman correlation (Movie rank vs audience score):\")\n",
    "print(f\"ρ = {spearman_rho:.4f}\")\n",
    "print(f\"p-value = {p_spearman:.6f}\")\n",
    "\n",
    "# Kendall Correlation\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"KENDALL CORRELATION COEFFICIENT\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "print(\"\"\"\n",
    "WHEN TO USE:\n",
    "- Very small sample sizes\n",
    "- Ordinal data\n",
    "- Similar robustness to Spearman but different interpretation\n",
    "\"\"\")\n",
    "\n",
    "kendall_tau, p_kendall = stats.kendalltau(rankings['Movie_Rank'], rankings['Audience_Score'])\n",
    "print(f\"\\nKendall correlation (Movie rank vs audience score):\")\n",
    "print(f\"τ = {kendall_tau:.4f}\")\n",
    "print(f\"p-value = {p_kendall:.6f}\")\n",
    "\n",
    "# Comparison\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"CORRELATION METHODS COMPARISON\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "print(f\"Pearson r: {pearson_r:.4f} (parametric, assumes linearity)\")\n",
    "print(f\"Spearman ρ: {spearman_rho:.4f} (non-parametric, ranks)\")\n",
    "print(f\"Kendall τ: {kendall_tau:.4f} (non-parametric, ordinal associations)\")\n",
    "\n",
    "# Partial Correlation\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"PARTIAL CORRELATION\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "print(\"\"\"\n",
    "WHEN TO USE:\n",
    "- Controlling for confounding variables\n",
    "- Finding relationship between X and Y after removing effect of Z\n",
    "- Understanding unique contribution of each variable\n",
    "\"\"\")\n",
    "\n",
    "# Create data with confounding\n",
    "np.random.seed(42)\n",
    "n = 100\n",
    "Z = np.random.normal(100, 15, n) # Confounding variable\n",
    "X = 0.5*Z + np.random.normal(0, 5, n)\n",
    "Y = 0.6*Z + np.random.normal(0, 5, n)\n",
    "\n",
    "# Simple correlations\n",
    "r_xy = np.corrcoef(X, Y)[0, 1]\n",
    "r_xz = np.corrcoef(X, Z)[0, 1]\n",
    "r_yz = np.corrcoef(Y, Z)[0, 1]\n",
    "\n",
    "print(f\"\\nCorrelation X-Y: {r_xy:.4f}\")\n",
    "print(f\"Correlation X-Z: {r_xz:.4f}\")\n",
    "print(f\"Correlation Y-Z: {r_yz:.4f}\")\n",
    "\n",
    "# Partial correlation: correlation of X and Y controlling for Z\n",
    "# Using residuals approach\n",
    "residuals_x = X - (r_xz * Z / np.std(Z)) * np.std(X)\n",
    "residuals_y = Y - (r_yz * Z / np.std(Z)) * np.std(Y)\n",
    "partial_r = np.corrcoef(residuals_x, residuals_y)[0, 1]\n",
    "\n",
    "print(f\"\\nPartial correlation (X-Y | Z): {partial_r:.4f}\")\n",
    "print(\"The confounding Z variable reduced apparent correlation!\")\n",
    "\n",
    "# Visualization\n",
    "fig, axes = plt.subplots(2, 3, figsize=(15, 10))\n",
    "\n",
    "# Scatter plots with different correlations\n",
    "datasets = [\n",
    "(np.random.randn(100), np.random.randn(100), \"No correlation (r≈0)\"),\n",
    "(x, y_positive, \"Positive correlation (r≈0.87)\"),\n",
    "(x, y_negative, \"Negative correlation (r≈-0.89)\"),\n",
    "(np.random.randn(100), np.random.randn(100), \"Weak (r≈0)\"),\n",
    "(np.random.randn(100)*10, np.cumsum(np.random.randn(100)), \"Non-linear\"),\n",
    "(data['Age'], data['Income']/1000, \"Age vs Income (r≈1.00)\")\n",
    "]\n",
    "\n",
    "for idx, (x_data, y_data, title) in enumerate(datasets):\n",
    "row = idx // 3\n",
    "col = idx % 3\n",
    "\n",
    "r = np.corrcoef(x_data, y_data)[0, 1]\n",
    "axes[row, col].scatter(x_data, y_data, alpha=0.6, s=50)\n",
    "axes[row, col].set_title(f\"{title}\\nr = {r:.3f}\")\n",
    "axes[row, col].grid(True, alpha=0.3)\n",
    "\n",
    "# Add trend line if sufficient correlation\n",
    "if abs(r) > 0.3:\n",
    "z = np.polyfit(x_data, y_data, 1)\n",
    "p = np.poly1d(z)\n",
    "x_line = np.linspace(x_data.min(), x_data.max(), 100)\n",
    "axes[row, col].plot(x_line, p(x_line), \"r--\", alpha=0.8)\n",
    "\n",
    "plt.suptitle('Correlation Examples', fontsize=14, fontweight='bold')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "611ac39e",
   "metadata": {},
   "source": [
    "### 1.2 Correlation Matrices and Multivariate Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4493a192",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "print(\"=\"*50)\n",
    "print(\"CORRELATION MATRIX ANALYSIS\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "# Create dataset with multiple variables\n",
    "np.random.seed(42)\n",
    "data = pd.DataFrame({\n",
    "'Age': np.random.randint(20, 60, 100),\n",
    "'Salary': np.random.randint(30000, 150000, 100),\n",
    "'Experience': np.random.randint(0, 30, 100),\n",
    "'Performance': np.random.randint(1, 5, 100),\n",
    "'Education': np.random.randint(10, 16, 100)\n",
    "})\n",
    "\n",
    "# Add some correlations to make it realistic\n",
    "data['Salary'] = data['Experience'] * 3000 + data['Age'] * 500 + np.random.randint(-20000, 20000, 100)\n",
    "data['Performance'] = 4 - (data['Age'] / 20) + np.random.normal(0, 0.3, 100)\n",
    "\n",
    "print(\"Data:\")\n",
    "print(data.head())\n",
    "\n",
    "# Calculate correlation matrix\n",
    "corr_matrix = data.corr()\n",
    "\n",
    "print(\"\\nCorrelation Matrix:\")\n",
    "print(corr_matrix)\n",
    "\n",
    "# Find strongest correlations\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"STRONGEST CORRELATIONS\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "# Get upper triangle of correlation matrix\n",
    "mask = np.triu(np.ones_like(corr_matrix, dtype=bool), k=1)\n",
    "upper_corr = corr_matrix.where(mask)\n",
    "\n",
    "# Flatten and sort\n",
    "corr_pairs = []\n",
    "for i in range(len(corr_matrix.columns)):\n",
    "for j in range(i+1, len(corr_matrix.columns)):\n",
    "corr_pairs.append({\n",
    "'Var1': corr_matrix.columns[i],\n",
    "'Var2': corr_matrix.columns[j],\n",
    "'Correlation': corr_matrix.iloc[i, j]\n",
    "})\n",
    "\n",
    "corr_df = pd.DataFrame(corr_pairs).sort_values('Correlation', key=abs, ascending=False)\n",
    "print(\"\\nTop 5 correlations:\")\n",
    "print(corr_df.head())\n",
    "\n",
    "# Visualization: Heatmap\n",
    "plt.figure(figsize=(10, 8))\n",
    "sns.heatmap(corr_matrix, annot=True, fmt='.2f', cmap='coolwarm', center=0,\n",
    "square=True, linewidths=1, cbar_kws={\"shrink\": 0.8})\n",
    "plt.title('Correlation Matrix Heatmap', fontsize=14, fontweight='bold')\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "# Pairplot (all relationships)\n",
    "print(\"\\nGenerating pairplot...\")\n",
    "sns.pairplot(data)\n",
    "plt.suptitle('Pairplot: All Relationships', fontsize=14, fontweight='bold', y=0.995)\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "# Multicollinearity check\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"MULTICOLLINEARITY CHECK\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "print(\"\"\"\n",
    "MULTICOLLINEARITY:\n",
    "- When independent variables are highly correlated\n",
    "- Causes instability in regression coefficients\n",
    "- Rule of thumb: r > 0.7 suggests multicollinearity\n",
    "- Can be detected with VIF (Variance Inflation Factor)\n",
    "\"\"\")\n",
    "\n",
    "# Check for high correlations (excluding diagonal)\n",
    "high_corr = []\n",
    "for i in range(len(corr_matrix.columns)):\n",
    "for j in range(i+1, len(corr_matrix.columns)):\n",
    "if abs(corr_matrix.iloc[i, j]) > 0.7:\n",
    "high_corr.append((corr_matrix.columns[i], corr_matrix.columns[j],\n",
    "corr_matrix.iloc[i, j]))\n",
    "\n",
    "if high_corr:\n",
    "print(\"\\nHigh correlations (|r| > 0.7):\")\n",
    "for var1, var2, r in high_corr:\n",
    "print(f\" {var1} <-> {var2}: r = {r:.3f}\")\n",
    "else:\n",
    "print(\"\\nNo high correlations detected (multicollinearity not severe)\")\n",
    "\n",
    "# Calculate VIF manually\n",
    "print(\"\\n\" + \"-\"*50)\n",
    "print(\"VARIANCE INFLATION FACTOR (VIF)\")\n",
    "print(\"-\"*50)\n",
    "\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from numpy.linalg import inv\n",
    "\n",
    "X = data[['Experience', 'Age', 'Education']].values\n",
    "X_scaled = StandardScaler().fit_transform(X)\n",
    "corr_X = np.corrcoef(X_scaled.T)\n",
    "vif_matrix = inv(corr_X)\n",
    "\n",
    "print(f\"\\nVIF for each variable:\")\n",
    "for i, var in enumerate(['Experience', 'Age', 'Education']):\n",
    "vif = vif_matrix[i, i]\n",
    "print(f\" {var}: VIF = {vif:.3f}\")\n",
    "interpretation = \"PROBLEMATIC\" if vif > 10 else \"MODERATE CONCERN\" if vif > 5 else \"OK\"\n",
    "print(f\" → {interpretation}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fd4d31e",
   "metadata": {},
   "source": [
    "## SESSION 2: Linear Regression - Simple and Multiple"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c251ffd",
   "metadata": {},
   "source": [
    "### Duration: 2 hours"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25128846",
   "metadata": {},
   "source": [
    "### 2.1 Simple Linear Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b0b9fe1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy import stats\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.metrics import mean_squared_error, r2_score\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "print(\"=\"*50)\n",
    "print(\"SIMPLE LINEAR REGRESSION\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "print(\"\"\"\n",
    "GOAL:\n",
    "Predict Y from X using linear relationship: Y = β₀ + β₁X + ε\n",
    "\n",
    "COMPONENTS:\n",
    "- Y: Dependent variable (what we predict)\n",
    "- X: Independent variable (predictor)\n",
    "- β₀: Intercept (Y when X=0)\n",
    "- β₁: Slope (change in Y per unit X)\n",
    "- ε: Error/residual (unexplained variation)\n",
    "\n",
    "ASSUMPTIONS:\n",
    "1. Linearity: Relationship is linear\n",
    "2. Independence: Observations are independent\n",
    "3. Homoscedasticity: Constant variance of errors\n",
    "4. Normality: Errors are normally distributed\n",
    "\"\"\")\n",
    "\n",
    "# Real example: House price prediction\n",
    "np.random.seed(42)\n",
    "n = 100\n",
    "square_feet = np.random.uniform(1000, 5000, n)\n",
    "price = 50000 + 200*square_feet + np.random.normal(0, 50000, n)\n",
    "\n",
    "data = pd.DataFrame({\n",
    "'Square_Feet': square_feet,\n",
    "'Price': price\n",
    "})\n",
    "\n",
    "print(f\"\\nExample: House Price Prediction\")\n",
    "print(f\"Sample size: {n}\")\n",
    "print(data.head())\n",
    "\n",
    "# Fit regression using sklearn\n",
    "model = LinearRegression()\n",
    "X = data[['Square_Feet']].values\n",
    "y = data['Price'].values\n",
    "model.fit(X, y)\n",
    "\n",
    "# Get coefficients\n",
    "intercept = model.intercept_\n",
    "slope = model.coef_[0]\n",
    "\n",
    "print(f\"\\nRegression Equation:\")\n",
    "print(f\"Price = {intercept:.2f} + {slope:.2f} × Square_Feet\")\n",
    "\n",
    "# Predictions\n",
    "y_pred = model.predict(X)\n",
    "\n",
    "# Model evaluation\n",
    "r2 = r2_score(y, y_pred)\n",
    "rmse = np.sqrt(mean_squared_error(y, y_pred))\n",
    "mae = np.mean(np.abs(y - y_pred))\n",
    "\n",
    "print(f\"\\nModel Performance:\")\n",
    "print(f\"R² (coefficient of determination): {r2:.4f}\")\n",
    "print(f\" → {r2*100:.1f}% of variance in price explained by square feet\")\n",
    "print(f\"RMSE (root mean squared error): ${rmse:,.2f}\")\n",
    "print(f\"MAE (mean absolute error): ${mae:,.2f}\")\n",
    "\n",
    "# Statistical significance\n",
    "# Calculate t-statistics manually\n",
    "residuals = y - y_pred\n",
    "mse = np.sum(residuals**2) / (len(y) - 2)\n",
    "se_slope = np.sqrt(mse / np.sum((X - X.mean())**2))\n",
    "t_stat = slope / se_slope\n",
    "p_value = 2 * (1 - stats.t.cdf(abs(t_stat), len(y) - 2))\n",
    "\n",
    "print(f\"\\nStatistical Significance:\")\n",
    "print(f\"Slope coefficient: {slope:.2f}\")\n",
    "print(f\"Standard error: {se_slope:.4f}\")\n",
    "print(f\"t-statistic: {t_stat:.4f}\")\n",
    "print(f\"p-value: {p_value:.2e}\")\n",
    "\n",
    "if p_value < 0.05:\n",
    "print(\"→ Slope is SIGNIFICANTLY different from zero\")\n",
    "else:\n",
    "print(\"→ Slope is NOT significantly different from zero\")\n",
    "\n",
    "# Confidence interval for slope\n",
    "t_crit = stats.t.ppf(0.975, len(y) - 2)\n",
    "ci_lower = slope - t_crit * se_slope\n",
    "ci_upper = slope + t_crit * se_slope\n",
    "print(f\"95% CI for slope: [{ci_lower:.2f}, {ci_upper:.2f}]\")\n",
    "\n",
    "# Predict for new values\n",
    "new_sqft = np.array([[2500], [3500], [4000]])\n",
    "new_pred = model.predict(new_sqft)\n",
    "\n",
    "print(f\"\\nPredictions for new houses:\")\n",
    "for sqft, pred in zip(new_sqft.flatten(), new_pred):\n",
    "print(f\" {sqft:.0f} sq ft → ${pred:,.0f}\")\n",
    "\n",
    "# Visualization\n",
    "fig, axes = plt.subplots(1, 2, figsize=(14, 5))\n",
    "\n",
    "# Scatter plot with regression line\n",
    "axes[0].scatter(X, y, alpha=0.6, s=50, label='Actual data')\n",
    "axes[0].plot(X, y_pred, 'r-', linewidth=2, label='Regression line')\n",
    "axes[0].set_xlabel('Square Feet', fontsize=12)\n",
    "axes[0].set_ylabel('Price ($)', fontsize=12)\n",
    "axes[0].set_title(f'Linear Regression\\nR² = {r2:.3f}', fontsize=12, fontweight='bold')\n",
    "axes[0].legend()\n",
    "axes[0].grid(True, alpha=0.3)\n",
    "\n",
    "# Residuals plot\n",
    "axes[1].scatter(y_pred, residuals, alpha=0.6, s=50)\n",
    "axes[1].axhline(y=0, color='r', linestyle='--', linewidth=2)\n",
    "axes[1].set_xlabel('Predicted Values', fontsize=12)\n",
    "axes[1].set_ylabel('Residuals', fontsize=12)\n",
    "axes[1].set_title('Residual Plot', fontsize=12, fontweight='bold')\n",
    "axes[1].grid(True, alpha=0.3)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b1f0954",
   "metadata": {},
   "source": [
    "### 2.2 Multiple Linear Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ddff94a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.metrics import r2_score, mean_squared_error\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "print(\"=\"*50)\n",
    "print(\"MULTIPLE LINEAR REGRESSION\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "print(\"\"\"\n",
    "GOAL:\n",
    "Predict Y from multiple X variables\n",
    "Y = β₀ + β₁X₁ + β₂X₂ + β₃X₃ + ... + ε\n",
    "\n",
    "ADVANTAGES:\n",
    "- Use multiple predictors\n",
    "- Control for confounding variables\n",
    "- Improve predictions\n",
    "- Understand individual variable effects\n",
    "\n",
    "CHALLENGES:\n",
    "- Multicollinearity: predictors correlated\n",
    "- More parameters to estimate\n",
    "- Risk of overfitting with too many variables\n",
    "\"\"\")\n",
    "\n",
    "# Example: Salary prediction with multiple variables\n",
    "np.random.seed(42)\n",
    "n = 100\n",
    "data = pd.DataFrame({\n",
    "'Age': np.random.randint(25, 60, n),\n",
    "'Experience': np.random.randint(0, 35, n),\n",
    "'Education_Years': np.random.randint(12, 20, n),\n",
    "'Performance_Score': np.random.uniform(2, 5, n)\n",
    "})\n",
    "\n",
    "# Create salary with known relationships\n",
    "data['Salary'] = (\n",
    "20000 +\n",
    "1500 * data['Age'] +\n",
    "2000 * data['Experience'] +\n",
    "3000 * data['Education_Years'] +\n",
    "10000 * data['Performance_Score'] +\n",
    "np.random.normal(0, 10000, n)\n",
    ")\n",
    "\n",
    "print(f\"\\nData (Salary Prediction):\")\n",
    "print(data.head())\n",
    "\n",
    "# Fit multiple regression\n",
    "X = data[['Age', 'Experience', 'Education_Years', 'Performance_Score']]\n",
    "y = data['Salary']\n",
    "\n",
    "model = LinearRegression()\n",
    "model.fit(X, y)\n",
    "\n",
    "# Coefficients\n",
    "intercept = model.intercept_\n",
    "coefs = model.coef_\n",
    "\n",
    "print(f\"\\nRegression Equation:\")\n",
    "print(f\"Salary = {intercept:,.0f}\")\n",
    "for var, coef in zip(X.columns, coefs):\n",
    "print(f\" + {coef:.2f} × {var}\")\n",
    "\n",
    "# Model performance\n",
    "y_pred = model.predict(X)\n",
    "r2 = r2_score(y, y_pred)\n",
    "rmse = np.sqrt(mean_squared_error(y, y_pred))\n",
    "\n",
    "print(f\"\\nModel Performance:\")\n",
    "print(f\"R²: {r2:.4f} ({r2*100:.1f}% variance explained)\")\n",
    "print(f\"RMSE: ${rmse:,.0f}\")\n",
    "\n",
    "# Feature importance (standardized coefficients)\n",
    "print(f\"\\nFeature Importance (Standardized Coefficients):\")\n",
    "X_std = (X - X.mean()) / X.std()\n",
    "model_std = LinearRegression()\n",
    "model_std.fit(X_std, y)\n",
    "\n",
    "for var, coef in zip(X.columns, model_std.coef_):\n",
    "print(f\" {var}: {abs(coef):,.0f} (absolute importance)\")\n",
    "\n",
    "# Prediction with new data\n",
    "new_employee = pd.DataFrame({\n",
    "'Age': [35],\n",
    "'Experience': [10],\n",
    "'Education_Years': [16],\n",
    "'Performance_Score': [4.5]\n",
    "})\n",
    "\n",
    "predicted_salary = model.predict(new_employee)\n",
    "print(f\"\\nPredicted salary for new employee: ${predicted_salary[0]:,.0f}\")\n",
    "\n",
    "# Adjusted R²\n",
    "adj_r2 = 1 - (1 - r2) * (n - 1) / (n - X.shape[1] - 1)\n",
    "print(f\"\\nAdjusted R²: {adj_r2:.4f}\")\n",
    "print(f\"(Adjusted for number of predictors)\")\n",
    "\n",
    "# Visualization: Actual vs Predicted\n",
    "fig, ax = plt.subplots(figsize=(10, 6))\n",
    "ax.scatter(y, y_pred, alpha=0.6, s=50)\n",
    "ax.plot([y.min(), y.max()], [y.min(), y.max()], 'r--', linewidth=2)\n",
    "ax.set_xlabel('Actual Salary', fontsize=12)\n",
    "ax.set_ylabel('Predicted Salary', fontsize=12)\n",
    "ax.set_title(f'Multiple Regression: Actual vs Predicted\\nR² = {r2:.3f}',\n",
    "fontsize=12, fontweight='bold')\n",
    "ax.grid(True, alpha=0.3)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9001c2d",
   "metadata": {},
   "source": [
    "### 2.3 Regression Assumptions and Diagnostics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5f4c8b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy import stats\n",
    "from sklearn.linear_model import LinearRegression\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "print(\"=\"*50)\n",
    "print(\"REGRESSION DIAGNOSTICS\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "print(\"\"\"\n",
    "REGRESSION ASSUMPTIONS:\n",
    "1. Linearity: Y = f(X) is linear\n",
    "2. Independence: Errors are independent\n",
    "3. Homoscedasticity: Constant error variance\n",
    "4. Normality: Errors are normally distributed\n",
    "5. No multicollinearity: Predictors not highly correlated\n",
    "\n",
    "CHECKING ASSUMPTIONS:\n",
    "- Residual plots\n",
    "- Q-Q plot for normality\n",
    "- Breusch-Pagan test for heteroscedasticity\n",
    "- Durbin-Watson test for autocorrelation\n",
    "\"\"\")\n",
    "\n",
    "# Create dataset\n",
    "np.random.seed(42)\n",
    "X = np.linspace(0, 10, 100)\n",
    "y = 5 + 2*X + np.random.normal(0, 2, 100)\n",
    "\n",
    "# Fit model\n",
    "model = LinearRegression()\n",
    "model.fit(X.reshape(-1, 1), y)\n",
    "y_pred = model.predict(X.reshape(-1, 1))\n",
    "residuals = y - y_pred\n",
    "\n",
    "print(f\"\\nResidual Statistics:\")\n",
    "print(f\"Mean: {residuals.mean():.6f} (should be ~0)\")\n",
    "print(f\"Std Dev: {residuals.std():.4f}\")\n",
    "\n",
    "# Normality test (Shapiro-Wilk)\n",
    "stat_shapiro, p_shapiro = stats.shapiro(residuals)\n",
    "print(f\"\\nShapiro-Wilk Test for Normality:\")\n",
    "print(f\"p-value: {p_shapiro:.4f}\")\n",
    "print(f\"Residuals are {'NORMAL' if p_shapiro > 0.05 else 'NOT NORMAL'}\")\n",
    "\n",
    "# Homoscedasticity test (Breusch-Pagan)\n",
    "from scipy.stats import f as f_dist\n",
    "residuals_sq = residuals ** 2\n",
    "X_const = np.column_stack([np.ones(len(X)), X])\n",
    "model_aux = LinearRegression()\n",
    "model_aux.fit(X.reshape(-1, 1), residuals_sq)\n",
    "ss_res_aux = np.sum((residuals_sq - model_aux.predict(X.reshape(-1, 1)))**2)\n",
    "ss_tot_aux = np.sum((residuals_sq - residuals_sq.mean())**2)\n",
    "r2_aux = 1 - ss_res_aux / ss_tot_aux\n",
    "\n",
    "bp_stat = len(X) * r2_aux\n",
    "bp_pvalue = 1 - f_dist.cdf(bp_stat, 1, len(X) - 2)\n",
    "\n",
    "print(f\"\\nBreusch-Pagan Test for Homoscedasticity:\")\n",
    "print(f\"BP statistic: {bp_stat:.4f}\")\n",
    "print(f\"p-value: {bp_pvalue:.4f}\")\n",
    "print(f\"Variance is {'CONSTANT' if bp_pvalue > 0.05 else 'NOT CONSTANT'}\")\n",
    "\n",
    "# Diagnostic plots\n",
    "fig, axes = plt.subplots(2, 2, figsize=(12, 10))\n",
    "\n",
    "# Plot 1: Residuals vs Fitted\n",
    "axes[0, 0].scatter(y_pred, residuals, alpha=0.6)\n",
    "axes[0, 0].axhline(y=0, color='r', linestyle='--')\n",
    "axes[0, 0].set_xlabel('Fitted Values')\n",
    "axes[0, 0].set_ylabel('Residuals')\n",
    "axes[0, 0].set_title('Residuals vs Fitted')\n",
    "axes[0, 0].grid(True, alpha=0.3)\n",
    "\n",
    "# Plot 2: Q-Q plot\n",
    "stats.probplot(residuals, dist=\"norm\", plot=axes[0, 1])\n",
    "axes[0, 1].set_title('Normal Q-Q Plot')\n",
    "axes[0, 1].grid(True, alpha=0.3)\n",
    "\n",
    "# Plot 3: Scale-Location\n",
    "standardized_residuals = residuals / residuals.std()\n",
    "axes[1, 0].scatter(y_pred, np.sqrt(np.abs(standardized_residuals)), alpha=0.6)\n",
    "axes[1, 0].set_xlabel('Fitted Values')\n",
    "axes[1, 0].set_ylabel('√|Standardized Residuals|')\n",
    "axes[1, 0].set_title('Scale-Location Plot')\n",
    "axes[1, 0].grid(True, alpha=0.3)\n",
    "\n",
    "# Plot 4: Residuals histogram\n",
    "axes[1, 1].hist(residuals, bins=20, edgecolor='black', alpha=0.7)\n",
    "axes[1, 1].set_xlabel('Residuals')\n",
    "axes[1, 1].set_ylabel('Frequency')\n",
    "axes[1, 1].set_title('Distribution of Residuals')\n",
    "axes[1, 1].grid(True, alpha=0.3, axis='y')\n",
    "\n",
    "plt.suptitle('Regression Diagnostic Plots', fontsize=14, fontweight='bold')\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "# Detecting outliers\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"DETECTING OUTLIERS AND INFLUENTIAL POINTS\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "# Cook's distance\n",
    "n = len(X)\n",
    "p = 2 # number of parameters\n",
    "mse = np.sum(residuals**2) / (n - p)\n",
    "leverage = 1/n + (X - X.mean())**2 / np.sum((X - X.mean())**2)\n",
    "cooks_d = (residuals**2 / (p * mse)) * (leverage / (1 - leverage))\n",
    "\n",
    "print(f\"\\nCook's distance (measures influence of each point):\")\n",
    "outlier_threshold = 4 / (n - p)\n",
    "outliers = np.where(cooks_d > outlier_threshold)[0]\n",
    "\n",
    "if len(outliers) > 0:\n",
    "print(f\"Found {len(outliers)} influential points\")\n",
    "for idx in outliers:\n",
    "print(f\" Point {idx}: Cook's D = {cooks_d[idx]:.4f}\")\n",
    "else:\n",
    "print(\"No extreme outliers detected\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc52cb8d",
   "metadata": {},
   "source": [
    "## SESSION 3: Logistic Regression, Causality, and Model Selection"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18026d9e",
   "metadata": {},
   "source": [
    "### Duration: 2 hours"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9beeab64",
   "metadata": {},
   "source": [
    "### 3.1 Logistic Regression for Binary Classification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18a8dad0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score\n",
    "from sklearn.metrics import confusion_matrix, roc_auc_score, roc_curve\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "print(\"=\"*50)\n",
    "print(\"LOGISTIC REGRESSION\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "print(\"\"\"\n",
    "WHEN TO USE:\n",
    "- Binary classification problems (yes/no, success/failure)\n",
    "- Predicting probabilities of belonging to class 1\n",
    "- When you need interpretable model\n",
    "\n",
    "LINEAR REGRESSION LIMITATIONS:\n",
    "- Predicts continuous values (-∞ to +∞)\n",
    "- Cannot directly give probabilities (0 to 1)\n",
    "- Violates assumptions\n",
    "\n",
    "LOGISTIC REGRESSION SOLUTION:\n",
    "- Uses sigmoid function: P(Y=1) = 1 / (1 + e^(-z))\n",
    "- Outputs probability (0 to 1)\n",
    "- Uses threshold (usually 0.5) to classify\n",
    "\n",
    "ODDS AND LOG-ODDS:\n",
    "- Odds = P(Y=1) / P(Y=0)\n",
    "- Log-odds = ln(Odds)\n",
    "- Logistic regression models log-odds linearly\n",
    "\"\"\")\n",
    "\n",
    "# Example: Customer churn prediction\n",
    "np.random.seed(42)\n",
    "n = 200\n",
    "\n",
    "data = pd.DataFrame({\n",
    "'Months_Customer': np.random.randint(1, 60, n),\n",
    "'Monthly_Bill': np.random.uniform(20, 150, n),\n",
    "'Tech_Support': np.random.choice([0, 1], n),\n",
    "})\n",
    "\n",
    "# Create churn probability based on features\n",
    "prob_churn = (\n",
    "0.9 / (1 + np.exp(-(0.1*data['Months_Customer'] - 2))) +\n",
    "0.3 * (data['Monthly_Bill'] / 150) -\n",
    "0.2 * data['Tech_Support']\n",
    ")\n",
    "prob_churn = np.clip(prob_churn, 0, 1)\n",
    "data['Churn'] = (np.random.rand(n) < prob_churn).astype(int)\n",
    "\n",
    "print(f\"\\nExample: Customer Churn Prediction\")\n",
    "print(data.head())\n",
    "print(f\"\\nChurn rate: {data['Churn'].mean():.1%}\")\n",
    "\n",
    "# Fit logistic regression\n",
    "X = data[['Months_Customer', 'Monthly_Bill', 'Tech_Support']]\n",
    "y = data['Churn']\n",
    "\n",
    "model = LogisticRegression()\n",
    "model.fit(X, y)\n",
    "\n",
    "# Coefficients\n",
    "print(f\"\\nModel Coefficients:\")\n",
    "print(f\"Intercept: {model.intercept_[0]:.4f}\")\n",
    "for feature, coef in zip(X.columns, model.coef_[0]):\n",
    "print(f\" {feature}: {coef:.4f}\")\n",
    "\n",
    "# Predictions\n",
    "y_pred_prob = model.predict_proba(X)[:, 1]\n",
    "y_pred = model.predict(X)\n",
    "\n",
    "# Evaluation metrics\n",
    "accuracy = accuracy_score(y, y_pred)\n",
    "precision = precision_score(y, y_pred)\n",
    "recall = recall_score(y, y_pred)\n",
    "f1 = f1_score(y, y_pred)\n",
    "auc = roc_auc_score(y, y_pred_prob)\n",
    "\n",
    "print(f\"\\nModel Performance:\")\n",
    "print(f\"Accuracy: {accuracy:.4f}\")\n",
    "print(f\"Precision: {precision:.4f}\")\n",
    "print(f\"Recall: {recall:.4f}\")\n",
    "print(f\"F1-Score: {f1:.4f}\")\n",
    "print(f\"AUC-ROC: {auc:.4f}\")\n",
    "\n",
    "# Confusion matrix\n",
    "cm = confusion_matrix(y, y_pred)\n",
    "print(f\"\\nConfusion Matrix:\")\n",
    "print(cm)\n",
    "print(f\"True Negatives: {cm[0,0]}\")\n",
    "print(f\"False Positives: {cm[0,1]}\")\n",
    "print(f\"False Negatives: {cm[1,0]}\")\n",
    "print(f\"True Positives: {cm[1,1]}\")\n",
    "\n",
    "# Predictions for new data\n",
    "new_customer = pd.DataFrame({\n",
    "'Months_Customer': [6],\n",
    "'Monthly_Bill': [100],\n",
    "'Tech_Support': [0]\n",
    "})\n",
    "\n",
    "prob = model.predict_proba(new_customer)[0, 1]\n",
    "prediction = model.predict(new_customer)[0]\n",
    "\n",
    "print(f\"\\nPrediction for new customer:\")\n",
    "print(f\" Probability of churn: {prob:.2%}\")\n",
    "print(f\" Prediction: {'CHURN' if prediction == 1 else 'RETAIN'}\")\n",
    "\n",
    "# ROC Curve\n",
    "fpr, tpr, thresholds = roc_curve(y, y_pred_prob)\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(14, 5))\n",
    "\n",
    "# ROC Curve\n",
    "axes[0].plot(fpr, tpr, linewidth=2, label=f'AUC = {auc:.3f}')\n",
    "axes[0].plot([0, 1], [0, 1], 'r--', label='Random classifier')\n",
    "axes[0].set_xlabel('False Positive Rate')\n",
    "axes[0].set_ylabel('True Positive Rate')\n",
    "axes[0].set_title('ROC Curve')\n",
    "axes[0].legend()\n",
    "axes[0].grid(True, alpha=0.3)\n",
    "\n",
    "# Confusion Matrix Heatmap\n",
    "import seaborn as sns\n",
    "sns.heatmap(cm, annot=True, fmt='d', cmap='Blues', ax=axes[1],\n",
    "xticklabels=['No Churn', 'Churn'],\n",
    "yticklabels=['No Churn', 'Churn'])\n",
    "axes[1].set_ylabel('Actual')\n",
    "axes[1].set_xlabel('Predicted')\n",
    "axes[1].set_title('Confusion Matrix')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee33695c",
   "metadata": {},
   "source": [
    "### 3.2 Correlation vs Causation and Confounding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1390696e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "print(\"=\"*50)\n",
    "print(\"CORRELATION VS CAUSATION\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "print(\"\"\"\n",
    "CRITICAL DISTINCTION:\n",
    "\n",
    "CORRELATION:\n",
    "- Two variables move together\n",
    "- Can be observed in data\n",
    "- Does NOT imply causation\n",
    "\n",
    "CAUSATION:\n",
    "- Change in X causes change in Y\n",
    "- Requires causal mechanism\n",
    "- Cannot be determined from observational data alone\n",
    "\n",
    "CONFOUNDING:\n",
    "- Third variable affects both X and Y\n",
    "- Creates spurious correlation\n",
    "- Makes X appear to cause Y when it doesn't\n",
    "\"\"\")\n",
    "\n",
    "# Classic examples\n",
    "examples = [\n",
    "(\"Shoe size & Reading ability\", \"AGE (confound)\"),\n",
    "(\"Ice cream sales & Drowning deaths\", \"TEMPERATURE (confound)\"),\n",
    "(\"Number of firefighters & Damage from fire\", \"FIRE SEVERITY (confound)\"),\n",
    "]\n",
    "\n",
    "print(\"\\nClassic confounding examples:\")\n",
    "for corr, confound in examples:\n",
    "print(f\" Apparent: {corr}\")\n",
    "print(f\" Real confounder: {confound}\\n\")\n",
    "\n",
    "# Simulate confounding\n",
    "print(\"=\"*50)\n",
    "print(\"EXAMPLE: CONFOUNDING IN ACTION\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "np.random.seed(42)\n",
    "n = 200\n",
    "\n",
    "# Confounding variable\n",
    "age = np.random.normal(40, 10, n)\n",
    "\n",
    "# Exposure (influenced by confounder)\n",
    "exercise = age * 0.5 + np.random.normal(0, 5, n)\n",
    "\n",
    "# Outcome (influenced by confounder, NOT by exposure)\n",
    "heart_disease = (age * 0.8 + np.random.normal(0, 10, n)) > 40\n",
    "\n",
    "# But appears correlated\n",
    "correlation = np.corrcoef(exercise, heart_disease.astype(int))[0, 1]\n",
    "\n",
    "print(f\"\\nScenario: Does exercise prevent heart disease?\")\n",
    "print(f\"Observed correlation: {correlation:.4f}\")\n",
    "print(f\"Interpretation: Appears exercise {'prevents' if correlation < 0 else 'causes'} disease\")\n",
    "\n",
    "print(f\"\\nBUT the true story:\")\n",
    "print(f\"- Age affects both exercise AND heart disease\")\n",
    "print(f\"- Older people exercise less AND have more disease\")\n",
    "print(f\"- It's AGE, not exercise!\")\n",
    "\n",
    "# Simpson's Paradox\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"SIMPSON'S PARADOX\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "print(\"\"\"\n",
    "A trend appears in groups separately but REVERSES when groups are combined!\n",
    "\"\"\")\n",
    "\n",
    "# Example: Treatment effectiveness\n",
    "data = pd.DataFrame({\n",
    "'Group': ['Mild', 'Mild', 'Severe', 'Severe', 'Overall', 'Overall'],\n",
    "'Treatment': ['A', 'B', 'A', 'B', 'A', 'B'],\n",
    "'Success_Rate': [0.95, 0.85, 0.65, 0.55, 0.80, 0.70]\n",
    "})\n",
    "\n",
    "print(\"\\nTreatment success rates by severity:\")\n",
    "print(data.to_string(index=False))\n",
    "\n",
    "print(\"\\nWithin each group: Treatment A is better\")\n",
    "print(\"BUT: Treatment B looks better overall (because different sample sizes!)\")\n",
    "\n",
    "# Causal inference methods\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"ADDRESSING CONFOUNDING\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "print(\"\"\"\n",
    "RANDOMIZED EXPERIMENTS:\n",
    "- Gold standard for causal inference\n",
    "- Randomization balances confounders\n",
    "- Can establish causation\n",
    "\n",
    "OBSERVATIONAL DATA METHODS:\n",
    "1. Matching: Compare similar subjects\n",
    "2. Stratification: Analyze within strata\n",
    "3. Regression adjustment: Control statistically\n",
    "4. Instrumental variables: Find natural experiments\n",
    "5. Propensity score: Match on probability of exposure\n",
    "\"\"\")\n",
    "\n",
    "# Visualization\n",
    "fig, axes = plt.subplots(2, 2, figsize=(14, 10))\n",
    "\n",
    "# Scenario 1: Correlation with causation\n",
    "np.random.seed(42)\n",
    "x1 = np.random.randn(100)\n",
    "y1 = 2*x1 + np.random.randn(100)*0.3\n",
    "axes[0, 0].scatter(x1, y1, alpha=0.6)\n",
    "axes[0, 0].set_title(\"Scenario 1: True Causation\n",
    "X causes Y\", fontweight='bold')\n",
    "axes[0, 0].set_xlabel(\"X\")\n",
    "axes[0, 0].set_ylabel(\"Y\")\n",
    "axes[0, 0].grid(True, alpha=0.3)\n",
    "\n",
    "# Scenario 2: Confounding\n",
    "z = np.random.randn(100)\n",
    "x2 = 0.5*z + np.random.randn(100)*0.3\n",
    "y2 = 0.5*z + np.random.randn(100)*0.3\n",
    "axes[0, 1].scatter(x2, y2, alpha=0.6)\n",
    "axes[0, 1].set_title(\"Scenario 2: Confounding\n",
    "Z causes both X and Y\", fontweight='bold')\n",
    "axes[0, 1].set_xlabel(\"X\")\n",
    "axes[0, 1].set_ylabel(\"Y\")\n",
    "axes[0, 1].grid(True, alpha=0.3)\n",
    "\n",
    "# Scenario 3: No relationship\n",
    "x3 = np.random.randn(100)\n",
    "y3 = np.random.randn(100)\n",
    "axes[1, 0].scatter(x3, y3, alpha=0.6)\n",
    "axes[1, 0].set_title(\"Scenario 3: No Relationship\", fontweight='bold')\n",
    "axes[1, 0].set_xlabel(\"X\")\n",
    "axes[1, 0].set_ylabel(\"Y\")\n",
    "axes[1, 0].grid(True, alpha=0.3)\n",
    "\n",
    "# Scenario 4: Reverse causation\n",
    "y4 = np.random.randn(100)\n",
    "x4 = 1.5*y4 + np.random.randn(100)*0.3\n",
    "axes[1, 1].scatter(x4, y4, alpha=0.6)\n",
    "axes[1, 1].set_title(\"Scenario 4: Reverse Causation\n",
    "Y causes X\", fontweight='bold')\n",
    "axes[1, 1].set_xlabel(\"X\")\n",
    "axes[1, 1].set_ylabel(\"Y\")\n",
    "axes[1, 1].grid(True, alpha=0.3)\n",
    "\n",
    "plt.suptitle('Causal Scenarios', fontsize=14, fontweight='bold')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59adc7d1",
   "metadata": {},
   "source": [
    "### 3.3 Model Selection and Validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de3880e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from sklearn.model_selection import train_test_split, cross_val_score\n",
    "from sklearn.linear_model import LinearRegression, Ridge, Lasso\n",
    "from sklearn.metrics import mean_squared_error, r2_score\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "print(\"=\"*50)\n",
    "print(\"MODEL SELECTION AND VALIDATION\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "print(\"\"\"\n",
    "TRAIN-TEST SPLIT:\n",
    "- Divide data: typically 70-80% train, 20-30% test\n",
    "- Train model on training set\n",
    "- Evaluate on test set (never seen by model)\n",
    "- Prevents overfitting\n",
    "\n",
    "CROSS-VALIDATION:\n",
    "- k-fold: Divide into k chunks, train k times\n",
    "- Leave-one-out: n-fold cross-validation\n",
    "- Better use of data, more stable estimates\n",
    "\n",
    "MODEL COMPLEXITY:\n",
    "- Underfitting: Model too simple, high bias\n",
    "- Overfitting: Model too complex, high variance\n",
    "- Sweet spot: Good generalization\n",
    "\"\"\")\n",
    "\n",
    "# Create dataset\n",
    "np.random.seed(42)\n",
    "n = 200\n",
    "X = np.linspace(0, 10, n).reshape(-1, 1)\n",
    "y_true = np.sin(X.ravel())\n",
    "y = y_true + np.random.normal(0, 0.15, n)\n",
    "\n",
    "# Train-test split\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)\n",
    "\n",
    "print(f\"\\nData split:\")\n",
    "print(f\"Training set: {len(X_train)} samples\")\n",
    "print(f\"Test set: {len(X_test)} samples\")\n",
    "\n",
    "# Test different polynomial degrees\n",
    "degrees = [1, 2, 3, 5, 10, 20]\n",
    "train_errors = []\n",
    "test_errors = []\n",
    "\n",
    "for degree in degrees:\n",
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "\n",
    "# Create polynomial features\n",
    "poly = PolynomialFeatures(degree=degree)\n",
    "X_train_poly = poly.fit_transform(X_train)\n",
    "X_test_poly = poly.transform(X_test)\n",
    "\n",
    "# Train model\n",
    "model = LinearRegression()\n",
    "model.fit(X_train_poly, y_train)\n",
    "\n",
    "# Evaluate\n",
    "y_train_pred = model.predict(X_train_poly)\n",
    "y_test_pred = model.predict(X_test_poly)\n",
    "\n",
    "train_mse = mean_squared_error(y_train, y_train_pred)\n",
    "test_mse = mean_squared_error(y_test, y_test_pred)\n",
    "\n",
    "train_errors.append(train_mse)\n",
    "test_errors.append(test_mse)\n",
    "\n",
    "print(f\"\\nDegree {degree}:\")\n",
    "print(f\" Train MSE: {train_mse:.6f}\")\n",
    "print(f\" Test MSE: {test_mse:.6f}\")\n",
    "print(f\" Status: {'UNDERFITTING' if degree < 3 else 'GOOD FIT' if degree == 3 else 'OVERFITTING'}\")\n",
    "\n",
    "# Regularization: Ridge vs Lasso\n",
    "print(\"\\n\" + \"=\"*50)\n",
    "print(\"REGULARIZATION (Ridge and Lasso)\")\n",
    "print(\"=\"*50)\n",
    "\n",
    "print(\"\"\"\n",
    "Regularization adds penalty for complex models:\n",
    "- Ridge (L2): Reduces large coefficients\n",
    "- Lasso (L1): Shrinks some coefficients to zero (feature selection)\n",
    "- Elastic Net: Combination of Ridge and Lasso\n",
    "\"\"\")\n",
    "\n",
    "from sklearn.linear_model import Ridge, Lasso\n",
    "\n",
    "# Create polynomial features\n",
    "poly = PolynomialFeatures(degree=10)\n",
    "X_train_poly = poly.fit_transform(X_train)\n",
    "X_test_poly = poly.transform(X_test)\n",
    "\n",
    "models = {\n",
    "'Linear': LinearRegression(),\n",
    "'Ridge (α=1)': Ridge(alpha=1),\n",
    "'Ridge (α=10)': Ridge(alpha=10),\n",
    "'Lasso (α=0.1)': Lasso(alpha=0.1)\n",
    "}\n",
    "\n",
    "for name, model in models.items():\n",
    "model.fit(X_train_poly, y_train)\n",
    "train_score = model.score(X_train_poly, y_train)\n",
    "test_score = model.score(X_test_poly, y_test)\n",
    "\n",
    "print(f\"\\n{name}:\")\n",
    "print(f\" Train R²: {train_score:.4f}\")\n",
    "print(f\" Test R²: {test_score:.4f}\")\n",
    "\n",
    "# Visualization\n",
    "fig, axes = plt.subplots(1, 2, figsize=(14, 5))\n",
    "\n",
    "# Plot 1: Bias-Variance Tradeoff\n",
    "axes[0].plot(degrees, train_errors, 'b-o', linewidth=2, label='Training error')\n",
    "axes[0].plot(degrees, test_errors, 'r-o', linewidth=2, label='Test error')\n",
    "axes[0].set_xlabel('Model Complexity (Polynomial Degree)')\n",
    "axes[0].set_ylabel('Mean Squared Error')\n",
    "axes[0].set_title('Bias-Variance Tradeoff', fontweight='bold')\n",
    "axes[0].legend()\n",
    "axes[0].grid(True, alpha=0.3)\n",
    "\n",
    "# Plot 2: Predictions for optimal model\n",
    "poly3 = PolynomialFeatures(degree=3)\n",
    "X_plot_poly = poly3.fit_transform(np.linspace(0, 10, 100).reshape(-1, 1))\n",
    "model_opt = LinearRegression()\n",
    "X_train_poly3 = poly3.fit_transform(X_train)\n",
    "model_opt.fit(X_train_poly3, y_train)\n",
    "y_plot_pred = model_opt.predict(X_plot_poly)\n",
    "\n",
    "axes[1].scatter(X_train, y_train, alpha=0.5, label='Training data')\n",
    "axes[1].scatter(X_test, y_test, alpha=0.5, label='Test data')\n",
    "axes[1].plot(np.linspace(0, 10, 100), y_true, 'g-', linewidth=2, label='True function')\n",
    "axes[1].plot(np.linspace(0, 10, 100), y_plot_pred, 'r--', linewidth=2, label='Degree 3 model')\n",
    "axes[1].set_xlabel('X')\n",
    "axes[1].set_ylabel('Y')\n",
    "axes[1].set_title('Optimal Model Fit', fontweight='bold')\n",
    "axes[1].legend()\n",
    "axes[1].grid(True, alpha=0.3)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9cddeebb",
   "metadata": {},
   "source": [
    "## Week 8 Summary\n",
    "By completing Week 8, you have learned:\n",
    "- Correlation: measures linear relationships between variables\n",
    "- Pearson correlation: parametric, assumes normality\n",
    "- Spearman correlation: non-parametric, rank-based\n",
    "- Kendall correlation: for small samples and ordinal data\n",
    "- Partial correlation: controlling for confounding variables\n",
    "- Correlation matrices and multicollinearity detection\n",
    "- VIF (Variance Inflation Factor) for detecting multicollinearity\n",
    "- Simple linear regression: predicting Y from single X\n",
    "- Multiple linear regression: predicting Y from multiple X\n",
    "- R² and adjusted R²: model fit measures\n",
    "- RMSE and MAE: prediction error metrics\n",
    "- Regression assumptions: linearity, independence, homoscedasticity, normality\n",
    "- Diagnostic plots: residuals, Q-Q, scale-location\n",
    "- Detecting outliers: Cook&#x27;s distance, leverage\n",
    "- Logistic regression: binary classification\n",
    "- Probabilities and threshold selection\n",
    "- Confusion matrix and performance metrics\n",
    "- ROC curve and AUC\n",
    "- Correlation vs causation: critical distinction\n",
    "- Confounding bias and Simpson&#x27;s Paradox\n",
    "- Methods for causal inference from observational data\n",
    "- Train-test split and cross-validation\n",
    "- Bias-variance tradeoff and overfitting\n",
    "- Regularization: Ridge and Lasso regression"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a04a130",
   "metadata": {},
   "source": [
    "## Week 8 Assignments"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6625b0a2",
   "metadata": {},
   "source": [
    "### Assignment 1: Correlation and Regression Analysis\n",
    "Analyze relationships in a dataset:\n",
    "- Calculate all correlation coefficients (Pearson, Spearman, Kendall)\n",
    "- Create correlation matrix and heatmap\n",
    "- Check for multicollinearity (VIF analysis)\n",
    "- Build simple linear regression model\n",
    "- Build multiple linear regression model\n",
    "- Compare model performance (R², RMSE)\n",
    "- Create diagnostic plots and interpret\n",
    "- Report findings professionally"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5524f9eb",
   "metadata": {},
   "source": [
    "### Assignment 2: Classification Model Development\n",
    "Build and evaluate logistic regression:\n",
    "- Prepare binary classification dataset\n",
    "- Split into train and test sets\n",
    "- Train logistic regression model\n",
    "- Calculate all performance metrics (accuracy, precision, recall, F1, AUC)\n",
    "- Create confusion matrix\n",
    "- Plot ROC curve\n",
    "- Interpret coefficients and make predictions\n",
    "- Compare with baseline model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e99b9e18",
   "metadata": {},
   "source": [
    "### Assignment 3: Model Selection and Validation\n",
    "Perform comprehensive model evaluation:\n",
    "- Implement train-test split with different ratios\n",
    "- Perform k-fold cross-validation\n",
    "- Test different model complexities (polynomial degrees)\n",
    "- Compare regularized models (Ridge, Lasso)\n",
    "- Visualize bias-variance tradeoff\n",
    "- Identify optimal model\n",
    "- Evaluate on held-out test set\n",
    "- Discuss overfitting and underfitting"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14a42283",
   "metadata": {},
   "source": [
    "## Recommended Practice Problems\n",
    "- Calculate and interpret all correlation types on same data\n",
    "- Detect and address multicollinearity in multiple regression\n",
    "- Build models with different feature combinations\n",
    "- Practice threshold selection for classification\n",
    "- Implement cross-validation from scratch\n",
    "- Compare regularized vs unregularized models\n",
    "- Identify confounding variables in real datasets\n",
    "- Design controlled experiments vs observational analysis\n",
    "- Create publication-quality diagnostic plots"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f742ee22",
   "metadata": {},
   "source": [
    "## Additional Resources"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3d09315",
   "metadata": {},
   "source": [
    "### Books\n",
    "- Chapter 8: Regression - &quot;Python for Data Analysis&quot; by Wes McKinney\n",
    "- The Elements of Statistical Learning by Hastie, Tibshirani, Friedman\n",
    "- Causal Inference: The Mixtape by Scott Cunningham (free online)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0b4fbfb",
   "metadata": {},
   "source": [
    "### Online Resources\n",
    "- scikit-learn documentation: https://scikit-learn.org/\n",
    "- Correlation explained: https://en.wikipedia.org/wiki/Correlation\n",
    "- Causality guide: https://www.britannica.com/science/causality\n",
    "- Statistical Learning tutorials"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}