cmkl/fall-2024/math/mat-206/00030/MAT-206-00030.ipynb

1362 lines
301 KiB
Plaintext
Raw Permalink Normal View History

2024-11-29 23:54:08 +07:00
{
"cells": [
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"import statsmodels.api as sm\n",
"import statsmodels.formula.api as smf\n",
"\n",
"import seaborn as sns\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from sklearn.preprocessing import MinMaxScaler\n",
"\n",
"from scipy.stats import norm, uniform\n",
"from statsmodels.stats import outliers_influence"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [],
"source": [
"#Significance level\n",
"ALPHA = 0.11"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [],
"source": [
"def PlotSimpleRegression(data, variable, ax):\n",
"\n",
" data = data.copy()\n",
" data = data.sort_values(variable).reset_index(drop=True)\n",
"\n",
" # Scatterplot of the observations\n",
" sns.scatterplot(\n",
" data = data,\n",
" x=variable,\n",
" y=\"Life Ladder\",\n",
" ax=ax,\n",
" label=\"Observations\"\n",
" )\n",
"\n",
" # Plot predicted mean\n",
" ax.plot(\n",
" data[variable],\n",
" data[\"mean\"],\n",
" color=\"k\",\n",
" label=\"Prediction\"\n",
" )\n",
"\n",
" # Plot prediction interval\n",
" ax.fill_between(\n",
" data[variable],\n",
" data[\"obs_ci_lower\"],\n",
" data[\"obs_ci_upper\"],\n",
" color=\"rebeccapurple\",\n",
" alpha=0.5,\n",
" label=\"Prediction interval\"\n",
" )\n",
"\n",
" # Plot confidence interval\n",
" ax.fill_between(\n",
" data[variable],\n",
" data[\"mean_ci_lower\"],\n",
" data[\"mean_ci_upper\"],\n",
" color=\"pink\",\n",
" alpha=0.5,\n",
" label=\"Confidence interval\"\n",
" )\n",
"\n",
" ax.legend(frameon=False)\n",
" ax.spines[['right', 'top']].set_visible(False)\n",
"\n",
" return ax"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [],
"source": [
"def PlotCompareYHatY(data, ax):\n",
" ax.scatter(data[\"Life Ladder\"], data[\"mean\"], color=\"k\")\n",
"\n",
" ax.errorbar(\n",
" data[\"Life Ladder\"],\n",
" data[\"mean\"],\n",
" yerr=data[\"obs_ci_upper\"] - data[\"mean\"],\n",
" fmt=\"o\",\n",
" color=\"k\"\n",
" )\n",
"\n",
" ax.plot(\n",
" [data[\"Life Ladder\"].min(), data[\"Life Ladder\"].max()]\n",
" , [data[\"Life Ladder\"].min(), data[\"Life Ladder\"].max()]\n",
" , color='r'\n",
" , linestyle='--'\n",
" )\n",
"\n",
" ax.set_xlabel(r\"$Y$\")\n",
" ax.set_ylabel(r\"$\\hat{Y}$\")\n",
" ax.spines[['right', 'top']].set_visible(False)\n",
"\n",
" return ax"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Reading and preprocessing data"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [],
"source": [
"DataWhr2024 = pd.read_csv(\"DataWhr2024.csv\")\n",
"UnM49 = pd.read_csv(\"UnM49.csv\", sep=';')"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [],
"source": [
"DataWhr2024.loc[DataWhr2024[\"Country name\"].str.startswith(\"Hong\"), \"Country name\"] = \"Hong Kong\"\n",
"DataWhr2024.loc[DataWhr2024[\"Country name\"].str.startswith(\"Somaliland\"), \"Country name\"] = \"Somaliland\"\n",
"DataWhr2024.loc[DataWhr2024[\"Country name\"].str.startswith(\"Taiwan\"), \"Country name\"] = \"Taiwan\""
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [],
"source": [
"UnM49 = UnM49[['Country or Area', 'Sub-region Name', 'Region Name']]\n",
"UnM49 = UnM49.rename({'Country or Area':'Country name', 'Sub-region Name':'Subregion', 'Region Name':'Continent'}, axis=1)"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [],
"source": [
"UnM49.loc[97, \"Country name\"] = \"Bolivia\"\n",
"UnM49.loc[33, \"Country name\"] = \"Congo (Brazzaville)\"\n",
"UnM49.loc[34, \"Country name\"] = \"Congo (Kinshasa)\"\n",
"UnM49.loc[124, \"Country name\"] = \"Hong Kong\"\n",
"UnM49.loc[125, \"Country name\"] = \"Macao\"\n",
"UnM49.loc[126, \"Country name\"] = \"North Korea\"\n",
"UnM49.loc[145, \"Country name\"] = \"Iran\"\n",
"UnM49.loc[46, \"Country name\"] = \"Ivory Coast\"\n",
"UnM49.loc[133, \"Country name\"] = \"Laos\"\n",
"UnM49.loc[129, \"Country name\"] = \"South Korea\"\n",
"UnM49.loc[173, \"Country name\"] = \"Moldova\"\n",
"UnM49.loc[217, \"Country name\"] = \"Netherlands\"\n",
"UnM49.loc[175, \"Country name\"] = \"Russia\"\n",
"UnM49.loc[164, \"Country name\"] = \"Syria\"\n",
"UnM49.loc[26, \"Country name\"] = \"Tanzania\"\n",
"UnM49.loc[116, \"Country name\"] = \"United States\"\n",
"UnM49.loc[193, \"Country name\"] = \"United Kingdom\"\n",
"UnM49.loc[111, \"Country name\"] = \"Venezuela\"\n",
"UnM49.loc[140, \"Country name\"] = \"Vietnam\""
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [],
"source": [
"_ = pd.DataFrame(\n",
" {\n",
" \"Country name\": [\"Kosovo\", \"Somaliland\", \"Taiwan\"],\n",
" \"Subregion\": [\"Southern Europe\", \"Sub-Saharan Africa\", \"Eastern Asia\"],\n",
" \"Continent\": [\"Europe\", \"Africa\", \"Asia\"],\n",
" }\n",
")\n",
"\n",
"UnM49 = pd.concat([UnM49, _], axis=0)\n",
"UnM49 = UnM49.reset_index(drop=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Merging the datasets"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [],
"source": [
"# Data\n",
"Dat = pd.merge(DataWhr2024, UnM49)\n",
"\n",
"# Data of 2023\n",
"Dat2023 = Dat[Dat['year'] == 2023]\n",
"Dat2023 = Dat2023.reset_index(drop=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In a previous analysis, I found that Afghanistan behaves as a leverage point, while Botswana and Sri Lanka bahave as outliers. Thus, we will not consider these countries in our analyses"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>Country name</th>\n",
" <th>year</th>\n",
" <th>Life Ladder</th>\n",
" <th>Log GDP per capita</th>\n",
" <th>Social support</th>\n",
" <th>Healthy life expectancy at birth</th>\n",
" <th>Freedom to make life choices</th>\n",
" <th>Generosity</th>\n",
" <th>Perceptions of corruption</th>\n",
" <th>Positive affect</th>\n",
" <th>Negative affect</th>\n",
" <th>Subregion</th>\n",
" <th>Continent</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>Afghanistan</td>\n",
" <td>2023</td>\n",
" <td>1.446</td>\n",
" <td>NaN</td>\n",
" <td>0.368</td>\n",
" <td>55.2</td>\n",
" <td>0.228</td>\n",
" <td>NaN</td>\n",
" <td>0.738</td>\n",
" <td>0.261</td>\n",
" <td>0.460</td>\n",
" <td>Southern Asia</td>\n",
" <td>Asia</td>\n",
" </tr>\n",
" <tr>\n",
" <th>13</th>\n",
" <td>Botswana</td>\n",
" <td>2023</td>\n",
" <td>3.332</td>\n",
" <td>9.673</td>\n",
" <td>0.701</td>\n",
" <td>55.0</td>\n",
" <td>0.741</td>\n",
" <td>-0.264</td>\n",
" <td>0.814</td>\n",
" <td>0.657</td>\n",
" <td>0.247</td>\n",
" <td>Sub-Saharan Africa</td>\n",
" <td>Africa</td>\n",
" </tr>\n",
" <tr>\n",
" <th>115</th>\n",
" <td>Sri Lanka</td>\n",
" <td>2023</td>\n",
" <td>3.602</td>\n",
" <td>9.364</td>\n",
" <td>0.790</td>\n",
" <td>67.4</td>\n",
" <td>0.754</td>\n",
" <td>0.050</td>\n",
" <td>0.922</td>\n",
" <td>0.709</td>\n",
" <td>0.353</td>\n",
" <td>Southern Asia</td>\n",
" <td>Asia</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" Country name year Life Ladder Log GDP per capita Social support \\\n",
"0 Afghanistan 2023 1.446 NaN 0.368 \n",
"13 Botswana 2023 3.332 9.673 0.701 \n",
"115 Sri Lanka 2023 3.602 9.364 0.790 \n",
"\n",
" Healthy life expectancy at birth Freedom to make life choices \\\n",
"0 55.2 0.228 \n",
"13 55.0 0.741 \n",
"115 67.4 0.754 \n",
"\n",
" Generosity Perceptions of corruption Positive affect Negative affect \\\n",
"0 NaN 0.738 0.261 0.460 \n",
"13 -0.264 0.814 0.657 0.247 \n",
"115 0.050 0.922 0.709 0.353 \n",
"\n",
" Subregion Continent \n",
"0 Southern Asia Asia \n",
"13 Sub-Saharan Africa Africa \n",
"115 Southern Asia Asia "
]
},
"execution_count": 37,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Dat2023.loc[[0, 13, 115]]"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [],
"source": [
"Dat2023 = Dat2023.drop([0, 13, 115])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"More preprocessing"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [],
"source": [
"Y = Dat2023[\"Life Ladder\"]\n",
"\n",
"X = Dat2023[[\n",
" 'Log GDP per capita',\n",
" 'Social support',\n",
" 'Healthy life expectancy at birth',\n",
" 'Freedom to make life choices',\n",
" 'Generosity',\n",
" 'Perceptions of corruption',\n",
" 'Positive affect',\n",
" 'Negative affect'\n",
"]]\n",
"\n",
"X = sm.add_constant(X)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Q1"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: Life Ladder R-squared: 0.001\n",
"Model: OLS Adj. R-squared: -0.007\n",
"Method: Least Squares F-statistic: 0.07861\n",
"Date: Fri, 29 Nov 2024 Prob (F-statistic): 0.780\n",
"Time: 17:57:08 Log-Likelihood: -187.60\n",
"No. Observations: 127 AIC: 379.2\n",
"Df Residuals: 125 BIC: 384.9\n",
"Df Model: 1 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"const 5.6752 0.097 58.421 0.000 5.483 5.867\n",
"Generosity 0.1657 0.591 0.280 0.780 -1.004 1.336\n",
"==============================================================================\n",
"Omnibus: 10.908 Durbin-Watson: 1.906\n",
"Prob(Omnibus): 0.004 Jarque-Bera (JB): 6.234\n",
"Skew: -0.367 Prob(JB): 0.0443\n",
"Kurtosis: 2.200 Cond. No. 6.24\n",
"==============================================================================\n",
"\n",
"Notes:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
]
}
],
"source": [
"# Extracting the design matrix and response variable\n",
"XGenerosity = X[[\"const\", \"Generosity\"]].dropna()\n",
"YGenerosity = Y[XGenerosity.index]\n",
"\n",
"# Fit the linear regression model\n",
"Model1 = sm.OLS(YGenerosity, XGenerosity).fit()\n",
"print(Model1.summary())"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA0UAAAIoCAYAAACmmkCFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAACt1klEQVR4nOzdd3xT5eIG8OdkN90LWlZZBUG2KEMFGTIUfqAIDlRwK+CEK4iooCKggjhRLwriVcAriFwQFJAlqOyhVJYUkF0KhdJmn98fzUlzMtq0TZr1fD+f0uSck5M3oU3f57xLEEVRBBERERERUZRSBLsAREREREREwcRQREREREREUY2hiIiIiIiIohpDERERERERRTWGIiIiIiIiimoMRUREREREFNUYioiIiIiIKKoxFBERERERUVRjKCIiIiIioqjGUEREROWqX78+6tevH+xihA1BEHDTTTe5bT948CBuu+02ZGZmQqFQICkpKSjl87eJEydCEASsW7fOp+Nzc3MhCAKGDx8e8LIREfmCoYiIosauXbvw+OOPo3nz5khISIBGo0FGRgZuvvlmTJ8+HefOnQt2EcNKRSvC4a5+/frQ6XSVfrzVasXAgQPxww8/4NZbb8XLL7+McePG+bWM69atgyAIePzxx/16XiKiSKcKdgGIiALNZrPh+eefx/Tp06FUKtGlSxf06tULsbGxOHv2LH799VeMGTMGr7zyCvbv34/atWsHu8ghZ82aNcEuQljJycmBXq+XbTty5Aj27duHRx55BJ9++mnQykZERO4Yiogo4r344ouYPn062rVrh4ULF6Jx48Zux+zYsQNjx45FcXFxUMoY6ho1ahTsIoSVq666ym3byZMnAQC1atUKQomIiKgs7D5HRBHtwIEDeOutt5Ceno6VK1d6DEQA0K5dO6xatcrjuJk9e/bgrrvuQmZmJjQaDbKysvDkk0/i/PnzsuOcx0kcOnQIt912G5KTkxEbG4uePXti9+7dHp/77NmzePbZZ9G4cWNotVqkpaVh0KBB+OOPP9yOlcb2XLx4EaNGjULdunWhUqkwd+5cxzH/+9//0K1bNyQmJiImJgatW7fGjBkzYLFY3M63du1a9O3bF7Vq1YJWq0XNmjVx4403urVkuI4puummmzBp0iQAQLdu3SAIAgRBQP369WGz2ZCVlYXU1FQYjUaPr7lLly5QqVT4559/PO4HgKNHj0KhUKB79+4e95vNZqSlpaFu3bqw2WwAgIKCArz88sto3rw54uLikJCQgMaNG2PYsGE4evSo1+fyN9cxRfXr10fXrl0BAJMmTXK8XxMnTnQcYzKZMGPGDLRr1w6xsbGIj4/HjTfeiKVLlwakjCdPnsQrr7yCjh07okaNGtBqtahfvz5GjBiBs2fPenzM8ePHcffddyMlJQVxcXHo2rUrNmzY4PU5rFYrpk2bhsaNG0On06Fx48aYMmWK4//LE3//PhAR+YItRUQU0b744gtYrVY89thjSE9PL/d4lUr+sbh06VIMGTIECoUCAwYMQN26dbFv3z588MEH+PHHH/H7778jOTlZ9pjc3Fx07NgRV199NR588EEcPnwY33//Pbp164acnBzUrFnTcezhw4dx00034Z9//kGvXr0wcOBAnD17FosWLcKPP/6INWvWoEOHDrLzG41GdO/eHYWFhfi///s/qFQqxzlnzJiB0aNHIyUlBffccw9iY2OxdOlSjB49Ghs3bsTixYshCAIAYPny5ejfvz+SkpIwYMAAZGZm4ty5c9i9eze+/PJLPProo17fJ2mA/Pr16zFs2DBHYEpKSoJCocDDDz+Ml19+GYsWLcI999wje+z+/fuxceNG3HrrrahTp47X58jKykKXLl2wfv16/PPPP27H/vDDDzh//jzGjh0LhUIBURTRu3dv/P7777j++uvRp08fKBQKHD16FEuXLsV9992HrKwsr88XSM888wx27dqFL774Al27dnUEJum70WhEnz59sG7dOrRp0wYPPfQQzGYzli9fjgEDBuD999/HqFGj/FqmDRs2YPr06ejRowc6dOgAtVqNnTt3YtasWfjxxx+xY8cOJCYmOo4/deoUOnXqhBMnTqB3795o164dcnJycPPNN6Nbt24en+PRRx/F559/jgYNGmDkyJEwGAyYMWMGNm/e7PF4f/8+EBH5TCQiimDdunUTAYhr1qyp8GPz8vLEhIQEsXbt2mJubq5s3/z580UA4qhRoxzbjhw5IgIQAYhTp06VHT9hwgQRgDhlyhTZ9s6dO4tKpVJcuXKlbPv+/fvF+Ph4sWXLlrLtWVlZIgCxd+/eYlFRkWzfoUOHRJVKJdaoUUM8duyYY7vBYBBvuOEGEYA4b948x/bbb79dBCDu2rXL42t3fd6srCzZtldeeUUEIK5du9bt8SdOnBBVKpV40003ue0bM2aMCEBcsmSJ2z5Xs2fPFgGI06ZNc9s3aNAgEYD4xx9/iKIoinv27BEBiAMHDnQ71mAwiJcvXy73+cqSlZUlarVan44FIHbt2lW2be3atSIA8ZVXXnE7fvz48SIA8aWXXhJtNptj+6VLl8T27duLGo1GPHHiRLnPKz3HY489Vu6xZ86c8fiefPHFFyIA8fXXX5dtHzZsmMftn3zyiePn3vlnQSpL69atxcLCQsf2f/75R0xLSxMBiMOGDZOdy5+/D0REFcHuc0QU0U6fPg14Gcexbt06TJw4UfblPJPavHnzcOnSJUyZMsWtheGuu+5Cu3btsGDBArfzNmjQAP/6179k2x566CEAwNatWx3bdu7cic2bN2PYsGHo3bu37PgmTZrgkUcewd69ez12G3rzzTcRExMj2/b111/DYrFg9OjRqFu3rmO7VqvFtGnTAMBjtyLX8wBAamqq27aKqFWrFvr374/169fj0KFDju1msxnz5s1DZmYmbr311nLPc8cdd0Cn0+E///mPbPvFixexbNkytGnTBldffXW5r0er1SIuLq5KrylQbDYbZs2ahUaNGjm61kni4+Px8ssvw2QyYfHixX593ho1anh8T+677z4kJCRg9erVjm0mkwkLFy5EjRo1MHr0aNnxDz/8MLKzs93OM2/ePADAyy+/jNjYWMf22rVr4+mnn3Y73t+/D0REFcHuc0QUtdatW+cYF+NM6tL022+/AQB+//13HD582O04g8GAvLw85OXlIS0tzbG9TZs2UCjk15ykrl8XL150bJPOf+bMGdnYEslff/3l+N6iRQvHdp1Oh5YtW7odv3PnTln5nXXq1Ak6nQ67du1ybLvrrruwePFidOzYEffccw969OiBG2+8UfZaquKxxx7Dd999h9mzZ2Pq1KmAvTvi2bNnMX78eLeuip4kJibi//7v//DNN99g9+7daN26NQDgv//9L4xGI+677z7Hsc2aNUOrVq0wf/58/PPPPxg4cCBuuukmj/8foWT//v24cOECatWq5fHnUZoqXvp58KfFixfjk08+wY4dO3DhwgVYrVbHPmliCKmMBoMB3bt3d5uWXKFQ4Prrr8fBgwdl26UxdDfeeKPb83ra5u/fByKiimAoIqKIVrNmTeTk5ODkyZNuM4JJrUMAsGDBAtx9992y/fn5+QCADz/8sMznuHLliixIJCQkuB0jBQDnSqd0/uXLl2P58uVlnt9ZjRo1ZK0JkkuXLjlesytBEFCzZk2cOHHCsW3w4MFYsmQJZsyYgY8//hgffvghBEFAt27dMH36dLRp06bM112eXr16oUGDBvjiiy/w+uuvQ6VSYfbs2RAEwdFy5ov77rsP33zzDf7zn/84QtGXX34JpVIpG6+kUqnw888/Y+LEiVi0aJGjRSM9PR2jRo3Ciy++CKVSWaXXFAjSz8Gff/6JP//80+txrj8HVTV9+nSMGTMG6enp6NWrF+rUqeNobZk5c6ZskoyCggLA/rPniaefuYKCAigUCo8h29Px/v59ICKqiNC9dEZE5AedO3cG7LOsVZQUbvbu3QtRFL1+VXbwvnT+999/v8zzDxs2TPY4bxVA6Xxnzpxx2yeKIs6cOeMW2AYMGID169fjwoULWLFiBR5++GGsW7c
"text/plain": [
"<Figure size 1000x600 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Ensure that PredictionTable only contains predictions for the same index as XGenerosity\n",
"PredictionTable1 = Model1.get_prediction(XGenerosity).summary_frame(alpha=0.11)\n",
"\n",
"# Create the plot\n",
"fig, ax = plt.subplots(figsize=(10, 6))\n",
"\n",
"# Scatterplot of observations\n",
"sns.scatterplot(\n",
" x=XGenerosity[\"Generosity\"], \n",
" y=YGenerosity, \n",
" ax=ax, \n",
" label=\"Observations\"\n",
")\n",
"\n",
"# Plot the predicted mean (regression line)\n",
"ax.plot(\n",
" XGenerosity[\"Generosity\"], \n",
" PredictionTable1[\"mean\"], \n",
" color=\"k\", \n",
" label=\"Prediction (Regression Line)\"\n",
")\n",
"\n",
"# Get the min and max of the x-axis for full range\n",
"x_min, x_max = XGenerosity[\"Generosity\"].min(), XGenerosity[\"Generosity\"].max()\n",
"\n",
"# Create a smoother x-range for the prediction lines and intervals\n",
"x_smooth = np.linspace(x_min, x_max, 300)\n",
"\n",
"# Get the predictions for the smooth x-range\n",
"PredictionSmooth = Model1.get_prediction(sm.add_constant(x_smooth)).summary_frame(alpha=0.11)\n",
"\n",
"# Plot prediction intervals across the full x-range\n",
"ax.fill_between(\n",
" x_smooth, \n",
" PredictionSmooth[\"obs_ci_lower\"], \n",
" PredictionSmooth[\"obs_ci_upper\"], \n",
" color=\"rebeccapurple\", \n",
" alpha=0.5, \n",
" label=\"Prediction Interval (89%)\"\n",
")\n",
"\n",
"# Plot confidence intervals across the full x-range\n",
"ax.fill_between(\n",
" x_smooth, \n",
" PredictionSmooth[\"mean_ci_lower\"], \n",
" PredictionSmooth[\"mean_ci_upper\"], \n",
" color=\"pink\", \n",
" alpha=0.5, \n",
" label=\"Confidence Interval (89%)\"\n",
")\n",
"\n",
"# Customize the plot\n",
"ax.set_title(\"Generosity vs. Life Ladder\", fontsize=14)\n",
"ax.set_xlabel(\"Generosity\", fontsize=12)\n",
"ax.set_ylabel(\"Life Ladder\", fontsize=12)\n",
"ax.legend()\n",
"ax.spines[['right', 'top']].set_visible(False)\n",
"\n",
"plt.show()\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Q2"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: Life Ladder R-squared: 0.271\n",
"Model: OLS Adj. R-squared: 0.266\n",
"Method: Least Squares F-statistic: 49.54\n",
"Date: Fri, 29 Nov 2024 Prob (F-statistic): 9.33e-11\n",
"Time: 17:57:08 Log-Likelihood: -177.57\n",
"No. Observations: 135 AIC: 359.1\n",
"Df Residuals: 133 BIC: 364.9\n",
"Df Model: 1 \n",
"Covariance Type: nonrobust \n",
"===================================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"-----------------------------------------------------------------------------------\n",
"const 2.2346 0.496 4.503 0.000 1.253 3.216\n",
"Positive affect 5.2694 0.749 7.039 0.000 3.789 6.750\n",
"==============================================================================\n",
"Omnibus: 6.132 Durbin-Watson: 1.815\n",
"Prob(Omnibus): 0.047 Jarque-Bera (JB): 4.574\n",
"Skew: -0.327 Prob(JB): 0.102\n",
"Kurtosis: 2.379 Cond. No. 13.7\n",
"==============================================================================\n",
"\n",
"Notes:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
]
}
],
"source": [
"# Create the design matrix XPossitive and the response vector YPossitive\n",
"XPossitive = X[['const', 'Positive affect']].dropna() # Drop missing values from XPossitive\n",
"YPossitive = Y.loc[XPossitive.index] # Align Y with XPossitive, matching indices\n",
"\n",
"# Fit the linear regression model (Model 2)\n",
"Model2 = sm.OLS(YPossitive, XPossitive).fit()\n",
"print(Model2.summary())"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA0UAAALCCAYAAAARRXhhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd3hU1dbA4d/0zGQmnRYMHREUEEW9iKIoAhY+UARRroDdC4qoqGBFEQUF7KBeFCyI3iuKXKyAUsSCjdBCCCENSO+TZPr5/kgyEjIhhZSZZL3Pk0c5+8w5e85MJmfNXnttlaIoCkIIIYQQQgjRRqlbugNCCCGEEEII0ZIkKBJCCCGEEEK0aRIUCSGEEEIIIdo0CYqEEEIIIYQQbZoERUIIIYQQQog2TYIiIYQQQgghRJsmQZEQQgghhBCiTZOgSAghhBBCCNGmSVAkhBBCCCGEaNMkKBJCiFokJyejUqmYNm1avR6nUqm49NJLm6xfTc3pdDJv3jx69+6NwWBApVKxbt26WtvEqZs3bx4qlYotW7ZUa3v11Vc588wzMZlMqFQqXn755RbpY2Or7+/LtGnTUKlUJCcnN2m/hBBtgwRFQgi/VhmQHP+j1+uJiYnhpptuYvfu3S3Wt0svvRSVStVi52+obdu2ea/lf//73xr3W7JkCU8//TTR0dHMnj2bp556ijPOOKPWtqZyskAhEKxatQqVSsXChQsbfIyPP/6Y++67D4PBwH333cdTTz3FP/7xj0btZ+X7OiMjo1GPK4QQ/kzb0h0QQoi66NmzJ//85z8BsFqt/PLLL6xZs4bPPvuMzZs3M3To0CY7d+fOnYmLiyM0NLRej4uLi8NkMjVZvxrqnXfegYpv5t99910mTJjgc78NGzZgNpvZuHEjer2+zm3i1N1zzz1MmjSJLl26VNm+YcMG73+jo6NbqHdCCNH6SFAkhAgIvXr1Yt68eVW2Pf744yxYsIDHHnusSUcPdDpdg0ZBmnrkpCGKior49NNPGTBgAB06dOC7774jLS2NmJiYavseO3aMyMhIn0HPydrEqYuKiiIqKqra9mPHjgFIQCSEEI1M0ueEEAHr3nvvBeC3337zbnO5XCxdupSBAwdiNBoJDQ1l+PDh/O9//6v2eI/Hw4oVKzj//POJiIjAaDRy2mmnMWbMmCpBlq85RSqViq1bt3r/v/LnxH2OnyNx2223oVKp2LZtm8/ns3TpUlQqFf/+97+rbN+9ezeTJk2iU6dO6PV6unbtyr333ktubm69r9maNWsoLS1lypQpTJkyBY/Hw6pVq6rsU5mmlpSUREpKive5devW7aRtx9u2bRtjxowhKioKg8FA7969efzxxyktLfXZr23btjFu3Dg6dOiAwWAgJiaG6667jh9//BEqUrqefvppAIYPH17jeU/UkGv+ww8/cOWVVxIdHY3BYKBDhw5cfPHFvP3223W8yqfuxFTBytS7H374AU54zx2vMd8rdfH5559z44030qtXL0wmE6GhoVx88cWsXbu2xsesWLGCs846i6CgIGJiYnj44Yex2Ww17r9v3z6uueYaLBYLoaGhXHXVVezdu/ek/friiy+4/PLLCQ8PJygoiLPOOovFixfjdrur7Fd5XVetWsX//vc/hg4disViqfV9JYRofWSkSAgR8CpvDBVF4frrr+eLL77g9NNPZ8aMGZSUlPDJJ5/wf//3fyxdupT777/f+7i5c+fywgsv0LNnT2666SYsFgtHjx7lxx9/ZNOmTSed9P3UU0+xatUqUlJSeOqpp7zbzz777Bofc/PNN/Puu+/y4YcfMmzYsGrtH3zwAQaDoUo62/r165k4cSJqtZqxY8cSExPD/v37ef311/n222/59ddfCQ8Pr/O1euedd9BoNEyePJmQkBD+9a9/sXLlSh5//HHvdax83pUT+GfNmgVAWFiY9/n5aqu0fPlyZsyYQVhYGGPGjKF9+/b8/vvvLFiwgB9++IEffvihygjTK6+8wv3334/RaOTaa6+lS5cu3tfh008/5aKLLvIGm1u3bmXq1Knem9bjz9sY1/zLL79kzJgxhIWFMXbsWDp16kR2djaxsbF88MEH3HnnnXW+1o3p7LPPrvE9V6mx3yt1MXfuXPR6PRdddJH3Wq1fv57rr7+eV1991fvFRaX58+fz5JNP0qFDB+644w50Oh2ffPIJcXFxPo+/d+9ehg4ditVq5brrrqN3797s3LmToUOHMnDgwBr7tHDhQjp37sx1111HaGgo27dv56GHHuLXX3/1OY/uv//9L9999x3XXHMN06dPp6ioqJGukBAiYChCCOHHkpKSFEAZNWpUtbYnn3xSAZThw4criqIo7733ngIol1xyiWK32737paSkKFFRUYpWq1USExO92yMiIpTo6GilpKSk2rFzc3Or9WHq1KlV9rnkkkuUk32MVvalksfjUbp06aKEh4crNputyr579uxRAOX666/3bsvJyVFCQkKUzp07K8nJyVX2X7NmjQIo99xzT43nP9Hu3burXcspU6YogLJp06Zq+3ft2lXp2rWrz2PV1LZv3z5Fq9UqAwcOVHJycqq0Pf/88wqgLF682Ltt165dilqtVqKjo5WkpKQq+3s8HuXo0aPefz/11FMKoPzwww91fs71vebXXXedAii7du2qdqwTn099rVy5UgGU559/vtZ9a3quNb3nGvO9UnmO9PT0Wvc9/vepUnFxsdK/f38lNDS0yu9WQkKCotVqlc6dOyuZmZne7YWFhUqfPn2q/b4c35cPP/ywyva5c+cqgAJUed9899133ve41Wr1bvd4PMrdd9+tAMqnn37q3V75mqjVamXjxo11uj5CiNZJ0ueEEAHh0KFDzJs3j3nz5vHQQw8xbNgwnnnmGYKCgliwYAEA7733HgAvvPBClZGILl26cP/99+NyuVi9enWV4+r1ejQaTbXzRURENPpzUKlUTJ48mfz8fL788ssqbR988AGAt5gEwPvvv09RURHPP/88Xbt2rbL/pEmTOOecc/j444/rfP7KAgtTpkzxbqv8/8q2U/XWW2/hcrl47bXXiIyMrNL28MMP065dO9asWVNlf4/Hw7PPPlstZUmlUp3y3Jn6XvNKRqOx2rYTn48/aez3Sl316NGj2jaz2cy0adMoLCysktr60Ucf4XK5eOCBB2jfvr13e0hICI8//ni146SmprJ161YGDBjA5MmTq7Q9+uijPkcJX3/9dQDefvttgoODvdsrq/6pVKoq779KY8eOZcSIEfV67kKI1kXS54QQASExMdE7p0Sn09GhQwduuukm5syZQ//+/QH466+/MJlMnH/++dUeP3z4cAB27drl3TZp0iSWLVvGWWedxaRJkxg+fDhDhgzxeUPcWG6++Waef/55PvjgA6677jqomNv00UcfERkZyVVXXeXd95dffgHg119/JTExsdqxbDYbOTk55OTk+JyUfzy73c6HH36IxWLh2muv9W4fPnw4MTExfP755+Tn559yelVln7/99ls2b95crV2n03HgwAHvv3fu3AnAyJEjT+m8J1Ofaz5p0iQ+++wz/vGPf3DTTTdx+eWXc/HFF9d6fVtaY75X6iMrK4uFCxfy9ddfk5KSQllZWZX2ysIQALGxsQBcfPHF1Y7ja1vl/hdddFG1NrPZzNlnn12twMovv/xCcHAw7777rs/+Go3GKu+/Sr4+M4QQbYsERUKIgDBq1Ci++eabk+5TVFTks4oaQKdOnbz7VHrllVfo3r07K1eu5Nlnn+XZZ58lKCiIiRMnsmTJkia5Ee7bty/nnnsuX331lTcI2bJlC0eOHGH69OnodDrvvnl5eQC88cYbJz1mSUlJrX1dt24dubm53HLLLVWCPrVazeTJk1m4cCEfffQRM2bMOKXnV9nnytG72hQWFqJSqbyvT1OozzWfMGEC69atY+nSpbz55pu88cYbqFQqhg8fzpIlS046Z6wlNeZ7pT7nPO+880hNTWXo0KGMGDGCsLAwNBoNu3bt4osvvsBut3v3LywsBKgySlSpQ4cO1badbP+aHpOXl4fL5fJ+geJLSUlJnY4lhGhbJH1
"text/plain": [
"<Figure size 1000x800 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Ensure that PredictionTable only contains predictions for the same index as XPossitive\n",
"PredictionTable2 = Model2.get_prediction(XPossitive).summary_frame(alpha=0.11)\n",
"\n",
"# Create the plot\n",
"fig, ax = plt.subplots(figsize=(10, 8))\n",
"\n",
"# Scatterplot of observations\n",
"sns.scatterplot(\n",
" x=XPossitive[\"Positive affect\"], \n",
" y=YPossitive, \n",
" ax=ax, \n",
" label=\"Observations\"\n",
")\n",
"\n",
"# Plot the predicted mean (regression line)\n",
"ax.plot(\n",
" XPossitive[\"Positive affect\"], \n",
" PredictionTable2[\"mean\"], \n",
" color=\"k\", \n",
" label=\"Prediction (Regression Line)\"\n",
")\n",
"\n",
"# Get the min and max of the x-axis for full range\n",
"x_min, x_max = XPossitive[\"Positive affect\"].min(), XPossitive[\"Positive affect\"].max()\n",
"\n",
"# Create a smoother x-range for the prediction lines and intervals\n",
"x_smooth = np.linspace(x_min, x_max, 300)\n",
"\n",
"# Get the predictions for the smooth x-range\n",
"PredictionSmooth = Model2.get_prediction(sm.add_constant(x_smooth)).summary_frame(alpha=0.11)\n",
"\n",
"# Plot prediction intervals across the full x-range\n",
"ax.fill_between(\n",
" x_smooth, \n",
" PredictionSmooth[\"obs_ci_lower\"], \n",
" PredictionSmooth[\"obs_ci_upper\"], \n",
" color=\"rebeccapurple\", \n",
" alpha=0.5, \n",
" label=\"Prediction Interval (89%)\"\n",
")\n",
"\n",
"# Plot confidence intervals across the full x-range\n",
"ax.fill_between(\n",
" x_smooth, \n",
" PredictionSmooth[\"mean_ci_lower\"], \n",
" PredictionSmooth[\"mean_ci_upper\"], \n",
" color=\"pink\", \n",
" alpha=0.5, \n",
" label=\"Confidence Interval (89%)\"\n",
")\n",
"\n",
"# Customize the plot\n",
"ax.set_title(\"Positive Affect vs. Life Ladder\", fontsize=14)\n",
"ax.set_xlabel(\"Positive Affect\", fontsize=12)\n",
"ax.set_ylabel(\"Life Ladder\", fontsize=12)\n",
"ax.legend()\n",
"ax.spines[['right', 'top']].set_visible(False)\n",
"\n",
"plt.show()\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Q3"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"['Log GDP per capita',\n",
" 'Social support',\n",
" 'Healthy life expectancy at birth',\n",
" 'Freedom to make life choices',\n",
" 'Generosity',\n",
" 'Perceptions of corruption',\n",
" 'Positive affect',\n",
" 'Negative affect']"
]
},
"execution_count": 44,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"['Log GDP per capita', 'Social support', 'Healthy life expectancy at birth', 'Freedom to make life choices', 'Generosity', 'Perceptions of corruption', 'Positive affect', 'Negative affect']"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: Life Ladder R-squared: 0.856\n",
"Model: OLS Adj. R-squared: 0.845\n",
"Method: Least Squares F-statistic: 80.73\n",
"Date: Fri, 29 Nov 2024 Prob (F-statistic): 2.99e-42\n",
"Time: 17:57:08 Log-Likelihood: -59.747\n",
"No. Observations: 118 AIC: 137.5\n",
"Df Residuals: 109 BIC: 162.4\n",
"Df Model: 8 \n",
"Covariance Type: nonrobust \n",
"====================================================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"----------------------------------------------------------------------------------------------------\n",
"const -2.5991 0.785 -3.309 0.001 -4.156 -1.042\n",
"Log GDP per capita 0.3141 0.086 3.638 0.000 0.143 0.485\n",
"Social support 3.2510 0.567 5.735 0.000 2.128 4.374\n",
"Healthy life expectancy at birth 0.0102 0.016 0.651 0.516 -0.021 0.041\n",
"Freedom to make life choices 1.3683 0.444 3.082 0.003 0.488 2.248\n",
"Generosity -0.4163 0.253 -1.646 0.103 -0.917 0.085\n",
"Perceptions of corruption -0.8887 0.269 -3.309 0.001 -1.421 -0.356\n",
"Positive affect 1.9932 0.461 4.322 0.000 1.079 2.907\n",
"Negative affect 1.0249 0.599 1.712 0.090 -0.162 2.212\n",
"==============================================================================\n",
"Omnibus: 4.969 Durbin-Watson: 2.170\n",
"Prob(Omnibus): 0.083 Jarque-Bera (JB): 4.505\n",
"Skew: -0.371 Prob(JB): 0.105\n",
"Kurtosis: 3.604 Cond. No. 1.51e+03\n",
"==============================================================================\n",
"\n",
"Notes:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
"[2] The condition number is large, 1.51e+03. This might indicate that there are\n",
"strong multicollinearity or other numerical problems.\n"
]
}
],
"source": [
"# Extract all covariates (independent variables) except the response variable\n",
"XAll = X[['const', 'Log GDP per capita', 'Social support', 'Healthy life expectancy at birth', 'Freedom to make life choices', \n",
" 'Generosity', 'Perceptions of corruption', 'Positive affect', 'Negative affect']].dropna()\n",
"\n",
"# Add a constant column to XAll for the intercept\n",
"XAll = sm.add_constant(XAll)\n",
"\n",
"# Ensure YAll is aligned with XAll\n",
"YAll = Y.loc[XAll.index]\n",
"\n",
"# Fit the linear regression model (Model 3)\n",
"Model3 = sm.OLS(YAll, XAll).fit()\n",
"\n",
"# Display the summary of Model 3\n",
"print(Model3.summary())\n"
]
},
{
"cell_type": "code",
"execution_count": 55,
"metadata": {},
"outputs": [],
"source": [
"PredictionTable3 = Model3.get_prediction(XAll).summary_frame(alpha=0.11)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Q4"
]
},
{
"cell_type": "code",
"execution_count": 47,
"metadata": {},
"outputs": [],
"source": [
"scaler = MinMaxScaler()\n",
"\n",
"X[\"Healthy life scaled\"] = scaler.fit_transform(X[[\"Healthy life expectancy at birth\"]])\n",
"\n",
"X[\"Log GDP scaled\"] = scaler.fit_transform(X[[\"Log GDP per capita\"]])"
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {},
"outputs": [],
"source": [
"XAllScaled = X.copy()\n",
"\n",
"XAllScaled = XAllScaled.drop([\"Healthy life expectancy at birth\", \"Log GDP per capita\"], axis=1)\n",
"\n",
"XAllScaled = XAllScaled.dropna()"
]
},
{
"cell_type": "code",
"execution_count": 49,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>const</th>\n",
" <th>Log GDP per capita</th>\n",
" <th>Social support</th>\n",
" <th>Healthy life expectancy at birth</th>\n",
" <th>Freedom to make life choices</th>\n",
" <th>Generosity</th>\n",
" <th>Perceptions of corruption</th>\n",
" <th>Positive affect</th>\n",
" <th>Negative affect</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>count</th>\n",
" <td>118.0</td>\n",
" <td>118.000000</td>\n",
" <td>118.000000</td>\n",
" <td>118.000000</td>\n",
" <td>118.000000</td>\n",
" <td>118.000000</td>\n",
" <td>118.000000</td>\n",
" <td>118.000000</td>\n",
" <td>118.000000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>mean</th>\n",
" <td>1.0</td>\n",
" <td>9.495669</td>\n",
" <td>0.790627</td>\n",
" <td>65.237288</td>\n",
" <td>0.796364</td>\n",
" <td>0.034373</td>\n",
" <td>0.722347</td>\n",
" <td>0.654653</td>\n",
" <td>0.293610</td>\n",
" </tr>\n",
" <tr>\n",
" <th>std</th>\n",
" <td>0.0</td>\n",
" <td>1.149838</td>\n",
" <td>0.131170</td>\n",
" <td>5.492634</td>\n",
" <td>0.113688</td>\n",
" <td>0.162590</td>\n",
" <td>0.173567</td>\n",
" <td>0.106431</td>\n",
" <td>0.088618</td>\n",
" </tr>\n",
" <tr>\n",
" <th>min</th>\n",
" <td>1.0</td>\n",
" <td>7.076000</td>\n",
" <td>0.398000</td>\n",
" <td>52.200000</td>\n",
" <td>0.452000</td>\n",
" <td>-0.268000</td>\n",
" <td>0.184000</td>\n",
" <td>0.344000</td>\n",
" <td>0.114000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>25%</th>\n",
" <td>1.0</td>\n",
" <td>8.612500</td>\n",
" <td>0.695750</td>\n",
" <td>60.700000</td>\n",
" <td>0.735250</td>\n",
" <td>-0.072500</td>\n",
" <td>0.663250</td>\n",
" <td>0.578250</td>\n",
" <td>0.229250</td>\n",
" </tr>\n",
" <tr>\n",
" <th>50%</th>\n",
" <td>1.0</td>\n",
" <td>9.636000</td>\n",
" <td>0.837500</td>\n",
" <td>66.100000</td>\n",
" <td>0.817500</td>\n",
" <td>0.022000</td>\n",
" <td>0.767500</td>\n",
" <td>0.667000</td>\n",
" <td>0.283000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>75%</th>\n",
" <td>1.0</td>\n",
" <td>10.470250</td>\n",
" <td>0.894250</td>\n",
" <td>69.650000</td>\n",
" <td>0.877000</td>\n",
" <td>0.134250</td>\n",
" <td>0.844500</td>\n",
" <td>0.738750</td>\n",
" <td>0.357500</td>\n",
" </tr>\n",
" <tr>\n",
" <th>max</th>\n",
" <td>1.0</td>\n",
" <td>11.676000</td>\n",
" <td>0.979000</td>\n",
" <td>74.600000</td>\n",
" <td>0.965000</td>\n",
" <td>0.590000</td>\n",
" <td>0.948000</td>\n",
" <td>0.843000</td>\n",
" <td>0.516000</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" const Log GDP per capita Social support \\\n",
"count 118.0 118.000000 118.000000 \n",
"mean 1.0 9.495669 0.790627 \n",
"std 0.0 1.149838 0.131170 \n",
"min 1.0 7.076000 0.398000 \n",
"25% 1.0 8.612500 0.695750 \n",
"50% 1.0 9.636000 0.837500 \n",
"75% 1.0 10.470250 0.894250 \n",
"max 1.0 11.676000 0.979000 \n",
"\n",
" Healthy life expectancy at birth Freedom to make life choices \\\n",
"count 118.000000 118.000000 \n",
"mean 65.237288 0.796364 \n",
"std 5.492634 0.113688 \n",
"min 52.200000 0.452000 \n",
"25% 60.700000 0.735250 \n",
"50% 66.100000 0.817500 \n",
"75% 69.650000 0.877000 \n",
"max 74.600000 0.965000 \n",
"\n",
" Generosity Perceptions of corruption Positive affect Negative affect \n",
"count 118.000000 118.000000 118.000000 118.000000 \n",
"mean 0.034373 0.722347 0.654653 0.293610 \n",
"std 0.162590 0.173567 0.106431 0.088618 \n",
"min -0.268000 0.184000 0.344000 0.114000 \n",
"25% -0.072500 0.663250 0.578250 0.229250 \n",
"50% 0.022000 0.767500 0.667000 0.283000 \n",
"75% 0.134250 0.844500 0.738750 0.357500 \n",
"max 0.590000 0.948000 0.843000 0.516000 "
]
},
"execution_count": 49,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"XAll.describe()"
]
},
{
"cell_type": "code",
"execution_count": 50,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>const</th>\n",
" <th>Social support</th>\n",
" <th>Freedom to make life choices</th>\n",
" <th>Generosity</th>\n",
" <th>Perceptions of corruption</th>\n",
" <th>Positive affect</th>\n",
" <th>Negative affect</th>\n",
" <th>Healthy life scaled</th>\n",
" <th>Log GDP scaled</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>count</th>\n",
" <td>118.0</td>\n",
" <td>118.000000</td>\n",
" <td>118.000000</td>\n",
" <td>118.000000</td>\n",
" <td>118.000000</td>\n",
" <td>118.000000</td>\n",
" <td>118.000000</td>\n",
" <td>118.000000</td>\n",
" <td>118.000000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>mean</th>\n",
" <td>1.0</td>\n",
" <td>0.790627</td>\n",
" <td>0.796364</td>\n",
" <td>0.034373</td>\n",
" <td>0.722347</td>\n",
" <td>0.654653</td>\n",
" <td>0.293610</td>\n",
" <td>0.582022</td>\n",
" <td>0.526015</td>\n",
" </tr>\n",
" <tr>\n",
" <th>std</th>\n",
" <td>0.0</td>\n",
" <td>0.131170</td>\n",
" <td>0.113688</td>\n",
" <td>0.162590</td>\n",
" <td>0.173567</td>\n",
" <td>0.106431</td>\n",
" <td>0.088618</td>\n",
" <td>0.245207</td>\n",
" <td>0.249965</td>\n",
" </tr>\n",
" <tr>\n",
" <th>min</th>\n",
" <td>1.0</td>\n",
" <td>0.398000</td>\n",
" <td>0.452000</td>\n",
" <td>-0.268000</td>\n",
" <td>0.184000</td>\n",
" <td>0.344000</td>\n",
" <td>0.114000</td>\n",
" <td>0.000000</td>\n",
" <td>0.000000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>25%</th>\n",
" <td>1.0</td>\n",
" <td>0.695750</td>\n",
" <td>0.735250</td>\n",
" <td>-0.072500</td>\n",
" <td>0.663250</td>\n",
" <td>0.578250</td>\n",
" <td>0.229250</td>\n",
" <td>0.379464</td>\n",
" <td>0.334022</td>\n",
" </tr>\n",
" <tr>\n",
" <th>50%</th>\n",
" <td>1.0</td>\n",
" <td>0.837500</td>\n",
" <td>0.817500</td>\n",
" <td>0.022000</td>\n",
" <td>0.767500</td>\n",
" <td>0.667000</td>\n",
" <td>0.283000</td>\n",
" <td>0.620536</td>\n",
" <td>0.556522</td>\n",
" </tr>\n",
" <tr>\n",
" <th>75%</th>\n",
" <td>1.0</td>\n",
" <td>0.894250</td>\n",
" <td>0.877000</td>\n",
" <td>0.134250</td>\n",
" <td>0.844500</td>\n",
" <td>0.738750</td>\n",
" <td>0.357500</td>\n",
" <td>0.779018</td>\n",
" <td>0.737880</td>\n",
" </tr>\n",
" <tr>\n",
" <th>max</th>\n",
" <td>1.0</td>\n",
" <td>0.979000</td>\n",
" <td>0.965000</td>\n",
" <td>0.590000</td>\n",
" <td>0.948000</td>\n",
" <td>0.843000</td>\n",
" <td>0.516000</td>\n",
" <td>1.000000</td>\n",
" <td>1.000000</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" const Social support Freedom to make life choices Generosity \\\n",
"count 118.0 118.000000 118.000000 118.000000 \n",
"mean 1.0 0.790627 0.796364 0.034373 \n",
"std 0.0 0.131170 0.113688 0.162590 \n",
"min 1.0 0.398000 0.452000 -0.268000 \n",
"25% 1.0 0.695750 0.735250 -0.072500 \n",
"50% 1.0 0.837500 0.817500 0.022000 \n",
"75% 1.0 0.894250 0.877000 0.134250 \n",
"max 1.0 0.979000 0.965000 0.590000 \n",
"\n",
" Perceptions of corruption Positive affect Negative affect \\\n",
"count 118.000000 118.000000 118.000000 \n",
"mean 0.722347 0.654653 0.293610 \n",
"std 0.173567 0.106431 0.088618 \n",
"min 0.184000 0.344000 0.114000 \n",
"25% 0.663250 0.578250 0.229250 \n",
"50% 0.767500 0.667000 0.283000 \n",
"75% 0.844500 0.738750 0.357500 \n",
"max 0.948000 0.843000 0.516000 \n",
"\n",
" Healthy life scaled Log GDP scaled \n",
"count 118.000000 118.000000 \n",
"mean 0.582022 0.526015 \n",
"std 0.245207 0.249965 \n",
"min 0.000000 0.000000 \n",
"25% 0.379464 0.334022 \n",
"50% 0.620536 0.556522 \n",
"75% 0.779018 0.737880 \n",
"max 1.000000 1.000000 "
]
},
"execution_count": 50,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"XAllScaled.describe()"
]
},
{
"cell_type": "code",
"execution_count": 51,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(np.float64(1508.9221), np.float64(43.9715))"
]
},
"execution_count": 51,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from numpy.linalg import cond\n",
"\n",
"condition_XAll = cond(XAll.values)\n",
"condition_XAllScaled = cond(XAllScaled.values)\n",
"\n",
"condition_XAll.round(4), condition_XAllScaled.round(4)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Q5"
]
},
{
"cell_type": "code",
"execution_count": 52,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: Life Ladder R-squared: 0.841\n",
"Model: OLS Adj. R-squared: 0.834\n",
"Method: Least Squares F-statistic: 119.1\n",
"Date: Fri, 29 Nov 2024 Prob (F-statistic): 2.26e-43\n",
"Time: 17:57:08 Log-Likelihood: -66.306\n",
"No. Observations: 119 AIC: 144.6\n",
"Df Residuals: 113 BIC: 161.3\n",
"Df Model: 5 \n",
"Covariance Type: nonrobust \n",
"================================================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"------------------------------------------------------------------------------------------------\n",
"const 5.6699 0.040 142.681 0.000 5.591 5.749\n",
"Social support 2.8341 0.538 5.266 0.000 1.768 3.900\n",
"Freedom to make life choices 1.3681 0.452 3.027 0.003 0.473 2.264\n",
"Perceptions of corruption -0.7368 0.272 -2.711 0.008 -1.275 -0.198\n",
"Positive affect 1.7803 0.472 3.773 0.000 0.845 2.715\n",
"Log GDP scaled 1.7166 0.289 5.940 0.000 1.144 2.289\n",
"==============================================================================\n",
"Omnibus: 1.443 Durbin-Watson: 2.108\n",
"Prob(Omnibus): 0.486 Jarque-Bera (JB): 0.998\n",
"Skew: -0.194 Prob(JB): 0.607\n",
"Kurtosis: 3.226 Cond. No. 15.6\n",
"==============================================================================\n",
"\n",
"Notes:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
]
}
],
"source": [
"XScaleFewVariables = X.drop([\"Healthy life expectancy at birth\", \"Log GDP per capita\"], axis=1)\n",
"\n",
"XScaleFewVariables = XScaleFewVariables.drop([\"const\", \"Healthy life scaled\", \"Generosity\", \"Negative affect\"], axis=1).dropna()\n",
"\n",
"XScaleFewVariables = XScaleFewVariables - XScaleFewVariables.mean()\n",
"\n",
"XScaleFewVariables = sm.add_constant(XScaleFewVariables)\n",
"\n",
"YScaleFewVariables = Y[XScaleFewVariables.index]\n",
"\n",
"Model4 = sm.OLS(YScaleFewVariables, XScaleFewVariables).fit()\n",
"\n",
"print(Model4.summary())\n"
]
},
{
"cell_type": "code",
"execution_count": 53,
"metadata": {},
"outputs": [],
"source": [
"PredictionTable4 = Model4.get_prediction(XScaleFewVariables).summary_frame(alpha=0.11)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Q6"
]
},
{
"cell_type": "code",
"execution_count": 54,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABW0AAAPeCAYAAAB3GThSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd3xT9f7H8Vea0gKtogxx/JDlYogTlKvI9ep1oV73FbeooGVYqy2glVmlFoUKWLUIuPe67nlduABxgOBgiHhVhijQUtqS5vdH2pKmGSfJSc5J8n4+Hnm0OTk555PkJDl5n+/5fh1ut9uNiIiIiIiIiIiIiNhCmtUFiIiIiIiIiIiIiMhOCm1FREREREREREREbEShrYiIiIiIiIiIiIiNKLQVERERERERERERsRGFtiIiIiIiIiIiIiI2otBWRERERERERERExEYU2oqIiIiIiIiIiIjYiEJbERERERERERERERtRaCsiIiIiIiIiIiJiIwptRUQSnMPhYMKECWHf76effsLhcPDggw/GpC4RERERESO0Pysi0pxCWxEREzz44IM4HA4cDgfz589vdrvb7aZTp044HA5OP/10S2qMxm233caZZ55Jx44dI96pFhERERH7Sub92e+++46CggIOPfRQdtllF/baay8GDRrEokWLrC5NRCQghbYiIiZq2bIljz/+eLPpH3zwAb/88guZmZmW1BWtwsJCFi5cyGGHHWZ1KSIiIiISQ8m4P/vAAw8we/ZsjjzySO666y7y8vL4/vvvOfroo3nnnXesLk9ExC+FtiIiJjrttNN45pln2LFjR5Ppjz/+OEcccQR77rmnZbVFY/Xq1fz22288+uijVpciIiIiIjGUjPuzgwcPZu3atTzwwAMMHTqU/Px8Pv/8c9q2baszyETEthTaioiYaPDgwfzxxx+8/fbbjdNqamp49tlnueiii/zep7KykhtvvJFOnTqRmZnJgQceyJ133onb7W4yX3V1NTfccAMdOnRgl1124cwzz+SXX37xu8z//e9/DBkyhI4dO5KZmUmvXr2YO3duxI+rS5cuEd9XRERERBJHMu7PHnHEEWRnZzeZ1q5dOwYMGMDy5csjWqaISKwptBURMVGXLl3o378/TzzxROO0119/nc2bN3PhhRc2m9/tdnPmmWcyffp0TjnlFKZNm8aBBx5Ifn4+eXl5Tea9+uqrKS0t5aSTTqK4uJgWLVowaNCgZstct25d46leI0aM4O6772a//fbjqquuorS0NEaPXERERESSQSrtz/7++++0b9/etOWJiJhJoa2IiMkuuugiXnzxRaqqqgB47LHHGDhwIHvvvXezeV966SX++9//MnnyZGbPns3w4cN56aWXOO+887j77rtZuXIlAF9//TWPPvooOTk5PPbYYwwfPpznnnuO3r17N1vmLbfcgsvl4ssvv+TWW2/l2muv5T//+Q8XXnghEyZMaKxLRERERMSfVNif/eijj/j000/597//HfWyRERiQaGtiIjJLrjgAqqqqnjllVfYunUrr7zySsBTyV577TWcTiejRo1qMv3GG2/E7Xbz+uuvN84HNJsvNze3yXW3281zzz3HGWecgdvtZuPGjY2Xk08+mc2bN7N48WKTH7GIiIiIJJNk359dv349F110EV27dqWgoCCqZYmIxEq61QWIiCSbDh06cOKJJ/L444+zbds2XC4X5513nt9516xZw957780uu+zSZHqPHj0ab2/4m5aWRvfu3ZvMd+CBBza5vmHDBv766y/Ky8spLy/3u87169dH9fhEREREJLkl8/5sZWUlp59+Olu3bmX+/PnN+roVEbELhbYiIjFw0UUXcc011/D7779z6qmnsttuu8VlvXV1dQBccsklXH755X7n6dOnT1xqEREREZHElYz7szU1NZxzzjl88803vPnmm367ZhARsQuFtiIiMXD22WczbNgwPvvsM5566qmA83Xu3Jl33nmHrVu3Nmmd8N133zXe3vC3rq6OlStXNmmN8P333zdZXsNIvC6XixNPPDEGj0xEREREUkGy7c/W1dVx2WWX8e677/L0008zcOBA05YtIhIL6tNWRCQGsrOzuffee5kwYQJnnHFGwPlOO+00XC4Xs2bNajJ9+vTpOBwOTj31VIDGvzNmzGgyn+/ouU6nk3PPPZfnnnuOpUuXNlvfhg0bonpcIiIiIpIakm1/duTIkTz11FOUlZVxzjnnRLQMEZF4UktbEZEYCXQ6l7czzjiD448/nltuuYWffvqJQw45hLfeeov//Oc/5ObmNvb5deihhzJ48GDKysrYvHkzf/vb33j33XdZsWJFs2UWFxfz3nvvcdRRR3HNNdfQs2dPNm3axOLFi3nnnXfYtGlT2I/lkUceYc2aNWzbtg2ADz/8kKKiIgAuvfTSxhYUIiIiIpI8kmV/trS0lLKyMvr370/r1q159NFHm9x+9tlnk5WVFdYyRURiTaGtiIiF0tLSeOmllxg3bhxPPfUU8+bNo0uXLkydOpUbb7yxybxz586lQ4cOPPbYY7z44ov84x//4NVXX6VTp05N5uvYsSMLFixg0qRJPP/885SVldGuXTt69erFHXfcEVGdc+bM4YMPPmi8/t577/Hee+8BcOyxxyq0FREREUlRibA/+9VXXwHw6aef8umnnza7ffXq1QptRcR2HG632211ESIiIiIiIiIiIiLioT5tRURERERERERERGxEoa2IiIiIiIiIiIiIjSi0FREREREREREREbERhbYiIiIiIiIiIiIiNqLQVkRERERERERERMRGFNqKiIiIiIiIiIiI2EjShbZut5stW7bgdrutLkVEREREJCLapxURERFJbUkX2m7dupU2bdqwdetWq0sREREREYmI9mlFREREUlvShbYiIiIiIiIiIiIiiUyhrYiIiIiIiIiIiIiNKLQVERERERERERERsRGFtiIiIiIiIiIiIiI2otBWRERERERERERExEYU2oqIiIiIiIiIiIjYiEJbERERERERERERERtRaCsiIiIiIiIiIiJiIwptRURERERERERERGxEoa2IiIiIiIiIiIiIjSi0FREREREREREREbERhbYiIiIiIiIiIiIiNqLQVkRERERERERERMRGFNqKiIiIiIiIiIiI2IhCWxEREREREREREREbUWgrIiIiIiIiIiIiYiMKbUVERERERERERERsRKGtiIiIiISnrg5mzYKXXrK6EhERERGRpJRudQEiIiIikkA2bYIzz4SPP4Y994TjjoPddrO6KhEREbGxyspKsrOzAaioqCArK8vqkkRsTy1tRURERMS43XaDFi0gOxtuvRV23dXqikREREREko5a2oqIiIhIcF9+CfvtB7vsAmlp8OCD4HDAvvtaXZmIiIiISFJSS1sRERER8a+qCsaMgb59YezYndM7d1ZgKyIiIiISQwptRURERKS5jz6CQw+FO+4Alws2bvT8FRERkaRQWVmJw+HA4XBQWVlpdTki4kOhrYiIiIjstHUrDB/uGWDshx9gr73ghRfgySfB6bS6OhEREUlBCpglFSm0FRERERGPRYugVy8oK/Ncv+oqWLYMzjrL6soSisvl4tZbb6Vr1660atWK7t27M3nyZNxut9WliYiIiEiC0EBkIiIiIuLRqRNUVEDXrjB7NpxwgtUVJaQ77riDe++9l4ceeohevXqxaNEirrzyStq0acOoUaOsLk9EkkBlZSXZ2dkAVFRUkJWVZXVJIiJiMoW2IiIiIqnK7YZPP4W//c1zvWNHeOMNT2tbBQAR++STT/jXv/7FoEGDAOjSpQtPPPEECxYssLo0EREREUkQ6h5BREREJBX99huccw4ccwy8+OLO6f36KbCN0t/+9jfeffddfvjhBwC+/vpr5s+fz6mnnmp1aSIiIiKSINTSVkREDNFpeCJJwu2GefMgLw82b4b0dFizxuqqksqYMWPYsmULBx10EE6nE5fLxW233cbFF18c8D7V1dVUV1c3Xt+yZUucqhURERERO1JLWxEREZFUsXo1nHSSZ4CxzZvhyCNh8WK4/nqrK0sqTz/9NI8
"text/plain": [
"<Figure size 1400x1000 with 4 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import matplotlib.pyplot as plt\n",
"\n",
"# Assuming PredictionTables and Y variables for each model are already calculated:\n",
"# Replace these with actual data (e.g., YModel1, PredictionTableModel1)\n",
"models = [\n",
" (\"Model 1\", YGenerosity, PredictionTable1),\n",
" (\"Model 2\", YPossitive, PredictionTable2),\n",
" (\"Model 3\", YAll, PredictionTable3),\n",
" (\"Model 4\", YScaleFewVariables, PredictionTable4),\n",
"]\n",
"\n",
"fig, axes = plt.subplots(2, 2, figsize=(14, 10))\n",
"axes = axes.ravel()\n",
"\n",
"for ax, (model_name, Y, PredictionTable) in zip(axes, models):\n",
" # Scatter plot of actual vs predicted values\n",
" ax.scatter(Y, PredictionTable[\"mean\"], color=\"k\")\n",
" \n",
" # Add error bars for prediction interval\n",
" yerr = PredictionTable[\"obs_ci_upper\"] - PredictionTable[\"mean\"]\n",
" ax.errorbar(Y, PredictionTable[\"mean\"], yerr=yerr, fmt=\"o\", color=\"k\")\n",
" \n",
" # Add identity line\n",
" ax.plot(\n",
" [Y.min(), Y.max()],\n",
" [Y.min(), Y.max()],\n",
" color=\"r\",\n",
" linestyle=\"--\",\n",
" )\n",
" \n",
" # Set labels and title\n",
" ax.set_xlabel(r\"$Y$\")\n",
" ax.set_ylabel(r\"$\\hat{Y}$\")\n",
" ax.set_title(model_name)\n",
" ax.spines[[\"right\", \"top\"]].set_visible(False)\n",
"\n",
"# Adjust layout\n",
"plt.tight_layout()\n",
"plt.show()\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.6"
}
},
"nbformat": 4,
"nbformat_minor": 2
}