{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 2.1 - Intermediate - Introduction to Regression\n",
"\n",
"COMET Team
*Emrul Hasan, Jonah Heyl, Shiming Wu, William Co, Avi\n",
"Woodward-Kelen, Anneke Dresselhuis, Rathin Dharani, Devan Rawlings,\n",
"Jasmine Arora, Jonathan Graves* \n",
"2024-07-25\n",
"\n",
"## Outline\n",
"\n",
"### Prerequisites\n",
"\n",
"- Basic R and Jupyter skills\n",
"- A theoretical understanding of simple linear relationship\n",
"- An understanding of hypothesis testing\n",
"- Types of variables (qualitative, quantitative)\n",
"\n",
"### Outcomes\n",
"\n",
"By the end of this notebook, you will be able to:\n",
"\n",
"- Learn how to run a simple linear regression using R\n",
"- Create and understand regression outputs in R\n",
"- Understand how to interpret coefficient estimates from simple linear\n",
" regressions in terms of an econometric model\n",
"- Examine the various elements of regression objects in R (including\n",
" fitted values, residuals and coefficients)\n",
"- Understand the relationship between $t$-tests and the estimates from\n",
" simple linear regressions\n",
"- Understand the role of qualitative variables in regression analysis,\n",
" with a particular emphasis on dummies\n",
"- Explain how adding variables to a model changes the results\n",
"\n",
"> Note that the data in this exercise is provided under the Statistics\n",
"> Canada Open License: \\> [1](#fn1s)Statistics\n",
"> Canada, Survey of Financial Security, 2019, 2021. Reproduced and\n",
"> distributed on an “as is” basis with the permission of Statistics\n",
"> Canada.Adapted from Statistics Canada, Survey of Financial Security,\n",
"> 2019, 2021. This does not constitute an endorsement by Statistics\n",
"> Canada of this product."
],
"id": "7ed055d1-a55e-4e2f-9f9c-d0a28a1cf85d"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"library(tidyverse)\n",
"library(haven)\n",
"library(dplyr)\n",
"library(scales)\n",
"source(\"intermediate_intro_to_regression_tests.r\")"
],
"id": "c987ad96-b99f-4db7-a415-f50a791d197c"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"SFS_data <- read_dta(\"../datasets_intermediate/SFS_2019_Eng.dta\") #this code is discussed in module 1\n",
"\n",
"SFS_data <- filter(SFS_data, !is.na(SFS_data$pefmtinc))\n",
"SFS_data <- rename(SFS_data, income_before_tax = pefmtinc)\n",
"SFS_data <- rename(SFS_data, income_after_tax = pefatinc)\n",
"SFS_data <- rename(SFS_data, wealth = pwnetwpg)\n",
"SFS_data <- rename(SFS_data, gender = pgdrmie)\n",
"SFS_data <- rename(SFS_data, education = peducmie)\n",
"\n",
"SFS_data <- SFS_data[!(SFS_data$education==\"9\"),]\n",
"SFS_data$education <- as.numeric(SFS_data$education)\n",
"SFS_data <- SFS_data[order(SFS_data$education),]\n",
"SFS_data$education <- as.character(SFS_data$education)\n",
"SFS_data$education[SFS_data$education == \"1\"] <- \"Less than high school\"\n",
"SFS_data$education[SFS_data$education == \"2\"] <- \"High school\"\n",
"SFS_data$education[SFS_data$education == \"3\"] <- \"Non-university post-secondary\"\n",
"SFS_data$education[SFS_data$education == \"4\"] <- \"University\"\n",
"\n",
"SFS_data$gender <- as_factor(SFS_data$gender)\n",
"SFS_data$education <- as_factor(SFS_data$education)"
],
"id": "3b9157aa-b3a9-4f9b-8a97-dc99a101b307"
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Part 1: Learning About Regressions\n",
"\n",
"What is a regression? What is the relationship of a regression to other\n",
"statistical concepts? How do we use regressions to answer economic\n",
"questions?\n",
"\n",
"In this notebook, we will explore these questions using data from the\n",
"2019 Survey of Financial Security (SFS). If you did module *1.0.1\n",
"Beginner Introduction to Statistics using R* you might recall we used it\n",
"to learn more about the gender wealth gap between male and female lead\n",
"households. You don’t need to work through the module in order to follow\n",
"along, but if you didn’t do it and you feel you might benefit from a\n",
"statistical refresher - now would be a great time for a little\n",
"refresher.\n",
"\n",
"Either way, we’ll begin our analysis by exploring the relationship\n",
"between wealth and income. Let’s start off with a visualization:"
],
"id": "18d6de00-2438-4a2a-9f6f-7722cdfba495"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"options(repr.plot.width=8,repr.plot.height=8) #controls the image size\n",
"f <- ggplot(data = SFS_data, xlim=c(0,2.4*10^6), ylim=c(0,3.4*10^7), aes(x = income_after_tax, y = wealth)) + \n",
" xlab(\"Income After Tax\") + \n",
" ylab(\"Wealth\") + scale_x_continuous()\n",
"\n",
"f + geom_point()"
],
"id": "0b9cb554-210a-4b55-9a47-a50249cea5e6"
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"> *Think Deeper*: What do you see here? Is there anything about this\n",
"> relationship that sticks out to you? Why does it have the shape it\n",
"> does?\n",
"\n",
"You can probably tell that there is definitely some relationship between\n",
"wealth and after-tax income - but it can be difficult to visualize using\n",
"a scatterplot alone. There are far too many points to make out a\n",
"discernable pattern or relationship here.\n",
"\n",
"### Regression Models\n",
"\n",
"This is where a **regression model** comes in. A regression model\n",
"specifies the relationship between two variables. For example, a linear\n",
"relationship would be:\n",
"\n",
"$$ W_i = \\beta_0 + \\beta_1I_i$$\n",
"\n",
"Where $W_i$ is wealth of family $i$, and $I_i$ is their after-tax\n",
"income. In econometrics, we typically refer to $W_i$ as the **outcome**\n",
"variable, and $I_i$ as the **explanatory** variable; you may have also\n",
"heard the terms *dependent* and *independent* variables respectively,\n",
"but these aren’t actually very good descriptions of what these variables\n",
"are in econometrics which is why we won’t use them here.\n",
"\n",
"A model like this is our description of what this relationship is - but\n",
"it depends on two unknowns: $\\beta_0$, $\\beta_1$.\n",
"\n",
"- $\\beta_0$ and $\\beta_1$ are **parameters** of the model: they are\n",
" numbers that determine the relationship (intercept and slope,\n",
" respectively) between $W_i$ and $I_i$\n",
"- This is a *linear* relationship because the model we have specified\n",
" uses coefficients that are characteristic of linear model formulas -\n",
" note that there are many other kinds of models beyond the linear\n",
" type seen here.\n",
"\n",
"It is unlikely, if not impossible, for the relationship we observe here\n",
"to completely explain everything about our data. We also need to include\n",
"a term which captures everything that is *not* described by the\n",
"relationship we described in the model. This is called the *residual*\n",
"term (meaning “leftover”).\n",
"\n",
"- The $\\epsilon_i$ is the **residual**: it is a component that\n",
" corresponds to the part of the data which is *not* described by the\n",
" model\n",
"- Residual terms will usually have certain assumed properties that\n",
" allow us to estimate the model\n",
"\n",
"Conceptually, you can think about a regression as two parts: the part of\n",
"the relationship explained by your model ($W_i = \\beta_0 + \\beta_1 I_i$)\n",
"and the part which is not explained ($\\epsilon_i$). The process of\n",
"“fitting” or estimating a regression model refers to finding values for\n",
"$\\beta_0$ and $\\beta_1$ such that as little as possible of the model is\n",
"explained by the residual term. We write the complete regression\n",
"equation by combining the two parts of the model:\n",
"\n",
"$$\n",
"W_i = \\beta_0 + \\beta_1 I_i + \\epsilon_i\n",
"$$\n",
"\n",
"The goal of regression analysis is to:\n",
"\n",
"1. Estimate this equation (and especially the model parameters) as\n",
" accurately as possible.\n",
"2. Learn about the relationship between $W_i$ and $I_i$ from the\n",
" results of that estimation\n",
"\n",
"There are many ways to define “as accurately as possible” and similarly\n",
"there are many ways to “estimate” the equation. In this course, we often\n",
"use *ordinary least squares* (OLS) as our estimation method which can be\n",
"understood as the following:\n",
"\n",
"$$\n",
"(\\hat{\\beta_0},\\hat{\\beta_1}) = \\arg \\min_{\\beta_0,\\beta_1} \\sum_{i=1}^{n} (W_i - \\beta_0 - \\beta_1 I_i)^2 =\\arg \\min_{b_0,b_1} \\sum_{i=1}^{n} (\\epsilon_i)^2\n",
"$$\n",
"\n",
"It is just the calculus way of writing “choose $\\beta_0$ and $\\beta_1$\n",
"(call them $\\hat{\\beta_0},\\hat{\\beta_1}$) such that they minimize the\n",
"sum of the squared residuals”. Ultimately, the goal of doing a\n",
"regression is to explain as much as possible using the parameters\n",
"($\\beta_0, \\beta_1$) and as little as possible using $\\epsilon_i$.\n",
"Through this equation, we have transformed our statistical problem into\n",
"a calculus problem, one that can can be solved, for example, by taking\n",
"derivatives.\n",
"\n",
"There are many, many ways to solve this estimation problem - most of\n",
"which are built into R. Before getting into how we can estimate using R\n",
"commands, we’ll discuss on how we can estimate manually.\n",
"\n",
"### Example: Manual Estimation\n",
"\n",
"If we think about the residuals as a gauge of error in our model\n",
"(remember we want to think about the error in absolute terms), we can\n",
"look at the scatterplot and guess how the model might perform based on\n",
"how small or large the residuals are from the regression line. As you\n",
"can probably imagine, this is not the most efficient nor the most\n",
"accurate way to solve our estimate problem!\n",
"\n",
"Try to get the best fit you can by playing around with the following\n",
"example."
],
"id": "9a9114a2-44ad-4804-80f3-917f217c3b89"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#set the value of B_0 and B_1 with these values\n",
"\n",
"B_0 <- 10000 #change me\n",
"B_1 <- 2 #change me\n",
"\n",
"# don't touch the rest of this code - but see if you can understand it!\n",
"SSE <- sum((SFS_data$wealth - B_0 - B_1*SFS_data$income_after_tax)^2) #sum of our squared errors\n",
"\n",
"SSE_rounded <- round(SSE/1000000,0) \n",
"print(paste(\"Your SSE is now,\", SSE_rounded,\", How low can you go?\")) #prints our SSE value\n",
"\n",
"options(repr.plot.width=10,repr.plot.height=8) #controls the image size\n",
"\n",
"fitted_line <- data.frame(income_before_tax = SFS_data$income_before_tax, wealth = B_0 + B_1*SFS_data$income_before_tax) #makes the regression line\n",
"\n",
"f <- ggplot(data = SFS_data,\n",
" aes(x = income_before_tax, y = wealth),\n",
" xlim=c(0,3000000),\n",
" ylim=c(0,30000000)) + \n",
" xlab(\"before tax income\") +\n",
" ylab(\"wealth\") +\n",
" scale_x_continuous(labels = comma) +\n",
" scale_y_continuous(labels = comma) \n",
"f <- f + geom_point() + geom_line(color = \"#330974\", data = fitted_line) #style preferences\n",
"\n",
"f #prints our graph with the line"
],
"id": "29bb6c87-2a75-48ac-8bf1-e22470705aed"
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As we change our $\\beta_0, \\beta_1$, notice how the best fit line\n",
"changes as well. The closer we fit our line to the data the lower SSE we\n",
"have.\n",
"\n",
"### Interactive Visualization of OLS\n",
"\n",
"Understanding OLS is fundamental to understanding regressions and other\n",
"topics in econometrics. Let’s try and understand the formula for OLS\n",
"above through a more visual approach.\n",
"\n",
"To demonstrate this, we will use a small scatter plot with just 4\n",
"points.\n",
"\n",
"- The straight line through the scatter plot is modelled by the simple\n",
" regression formula $B_0 + B_1X$.\n",
"\n",
"- Since it’s nearly impossible for a regression to perfectly predict\n",
" the relationship between two variables, we will almost always\n",
" include an **unobservable error** $e_i$ with our regression\n",
" estimation. This is the vertical distance between the regression\n",
" line and the actual data points.\n",
"\n",
"- Hence each of the points can be modelled by the equation\n",
" $Y_i = B_0 + B_1X + e_i$.\n",
"\n",
"- Instead of minimizing the error terms, we will try to minimize the\n",
" squared errors which are represented by the size of those red boxes.\n",
"\n",
"Try your own values for `beta_0` and `beta_1`. Make sure to try the\n",
"values only roughly within the specified range. The actual value of\n",
"`beta_0` and `beta_1` that minimize the residual sum of squares is 0.65\n",
"and 0.82 respectively. The code block below also displays the area of\n",
"the red boxes; deviation from these optimal values will increase the\n",
"area of the red boxes."
],
"id": "afeff29d-99e9-4193-8522-ddd85276e84e"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"beta_0 <- 0.65 #CHANGE THIS VALUE, TRY VALUES BETWEEN 0 - 1\n",
"beta_1 <- 0.82 #CHANGE THIS VALUE, TRY VALUES BETWEEN 0.6 - 1.4\n",
"\n",
"x <- c(1, 2, 3, 4)\n",
"y <- c(1.7, 1.5, 4, 3.6)\n",
"\n",
"# don't worry about this code, just run it!\n",
"dta <- data.frame(x, y)\n",
"example_df_graph <- dta %>%\n",
" ggplot(aes(x = x, y = y)) +\n",
" geom_point() +\n",
" geom_abline(intercept = beta_0, slope = beta_1) +\n",
" xlim(0, 5) +\n",
" ylim(0, 5) +\n",
" geom_rect(aes(xmin = (dta[1, \"x\"] + (beta_0 + (beta_1 * dta[1, \"x\"])) - dta[1, \"y\"]), xmax = dta[1, \"x\"], \n",
" ymin = (beta_0 + (beta_1 * dta[1, \"x\"])), ymax = dta[1, \"y\"]),\n",
" alpha = 0.1,\n",
" fill = \"red\") +\n",
" geom_rect(aes(xmin = dta[2, \"x\"], xmax = (dta[2, \"x\"] + (beta_0 + (beta_1 * dta[2, \"x\"])) - dta[2, \"y\"]), \n",
" ymin = dta[2, \"y\"], ymax = (beta_0 + (beta_1 * dta[2, \"x\"]))), \n",
" alpha = 0.1, \n",
" fill = \"red\") +\n",
" geom_rect(aes(xmin = (dta[3, \"x\"] + (beta_0 + (beta_1 * dta[3, \"x\"])) - dta[3, \"y\"]), xmax = dta[3, \"x\"], \n",
" ymin = (beta_0 + (beta_1 * dta[3, \"x\"])), ymax = dta[3, \"y\"]), \n",
" alpha = 0.1, \n",
" fill = \"red\") +\n",
" geom_rect(aes(xmin = dta[4, \"x\"], xmax = (dta[4, \"x\"] + (beta_0 + (beta_1 * dta[4, \"x\"])) - dta[4, \"y\"]), \n",
" ymin = dta[4, \"y\"], ymax = (beta_0 + (beta_1 * dta[4, \"x\"]))), \n",
" alpha = 0.1, \n",
" fill = \"red\")\n",
"example_df_graph\n",
"\n",
"area_1 <- ((dta[1, \"x\"] - (dta[1, \"x\"] + (beta_0 + (beta_1 * dta[1, \"x\"])) - dta[1, \"y\"])) * \n",
" ((beta_0 + (beta_1 * dta[2, \"x\"])) - dta[2, \"y\"]))\n",
"area_2 <- ((dta[2, \"x\"] + (beta_0 + (beta_1 * dta[2, \"x\"])) - dta[2, \"y\"]) - dta[2, \"x\"]) * \n",
" ((beta_0 + (beta_1 * dta[2, \"x\"])) - dta[2, \"y\"])\n",
"area_3 <- (dta[3, \"x\"] - (dta[3, \"x\"] + (beta_0 + (beta_1 * dta[3, \"x\"])) - dta[3, \"y\"])) * \n",
" (dta[3, \"y\"]) - (beta_0 + (beta_1 * dta[3, \"x\"]))\n",
"area_4 <- ((dta[4, \"x\"] + (beta_0 + (beta_1 * dta[4, \"x\"])) - dta[4, \"y\"]) - dta[4, \"x\"]) * \n",
" ((beta_0 + (beta_1 * dta[4, \"x\"])) - dta[4, \"y\"])\n",
"\n",
"area <- area_1 + area_2 + area_3 + area_4\n",
"print(\"Area of red boxes is: \")\n",
"area"
],
"id": "d71a1375-14e1-496e-8b28-05472b90e219"
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Simple Regressions in R\n",
"\n",
"Now, let’s see how we could use a regression in R to do this.\n",
"\n",
"- Regression models look like: `Y ~ X` where `Y` is regressed on `X`\n",
" and the `~` symbol is called “tilde” and is pronounced “TIL-duh” by\n",
" the way\n",
"\n",
"> For now you can ignore the residual terms and parameters when writing\n",
"> the model in R (the computer implicitly knows it’s there!) and just\n",
"> focus on the variables.\n",
"\n",
"So, for example, our regression model is\n",
"\n",
"$$\n",
"W_i = \\beta_0 + \\beta_1 I_i + \\epsilon_i\n",
"$$\n",
"\n",
"Which can be written in R as\n",
"\n",
"`wealth ~ income_before_tax`\n",
"\n",
"Regressions are estimated in R using the `lm` function which contains an\n",
"argument to specify the dataset.\n",
"\n",
"- This creates a *linear model* object, which can be used to calculate\n",
" things (using prediction) or perform tests.\n",
"- It also stores all of the information about the model, such as the\n",
" coefficient and fit.\n",
"- These models can also be printed and summarized to give important\n",
" basic information about a regression.\n",
"\n",
"Below are a few of the most important elements of a linear model. Let’s\n",
"say, for example, that we called the model `my_model.`\n",
"\n",
"- `my_model$coefficients`: gives us the parameter coefficients\n",
"- `my_model$residuals`: gives us the residuals\n",
"- `my_model$fitted.values`: gives us the predicted values\n",
"\n",
"Enough talk! Let’s see our model in action here."
],
"id": "cb8dbedb-4ee0-4f6b-bc08-2c5f95ea63f9"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"regression1 = lm(wealth ~ income_after_tax, data = SFS_data) # take note this is very important!\n",
"\n",
"summary(regression1)\n",
" \n",
"head(regression1$coefficients)"
],
"id": "c9527bc9-5bf2-428e-b78c-bcb14a9c237c"
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Take a close look at the results. Identify the following elements:\n",
"\n",
"- The values of the parameters\n",
"- The standard errors of the parameters\n",
"- The %-of the data explained by the model\n",
"\n",
"> **Test Your Knowledge**: What %-of the variance in wealth is explained\n",
"> by the model? \n",
"> Write the percentage in *decimal form* and include all decimals given\n",
"> by the model (example, x.xxx - where x are numbers)"
],
"id": "ce99cf2b-a23c-4150-8081-8e8f7acfa4e8"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"answer1 <- ??? #answer goes here\n",
"\n",
"test_1()"
],
"id": "8726aae6-75cc-43f6-a113-cbb63f049ca6"
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The underlying model and the parameters tells us about the relationship\n",
"between the different values:\n",
"\n",
"$$W_i = 169826.16 + 9.96 I_i + \\epsilon_i$$\n",
"\n",
"Notice, for example:\n",
"\n",
"$$\\frac{\\partial W_i}{\\partial I_i} = \\beta_1 = 9.96$$\n",
"\n",
"In other words, when incomes goes up by 1 dollar, we would expect that\n",
"the wealth accumulated for this given family will rise by 9.96 dollars.\n",
"This kind of analysis is key to interpreting what this model is telling\n",
"us.\n",
"\n",
"Finally, let’s visualize our fitted model on the scatterplot from\n",
"before. How does it compare to your original model?"
],
"id": "d18f9a11-96d5-445c-89d1-5fbcc1ed3373"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"options(repr.plot.width=10,repr.plot.height=8) #style preferences\n",
"\n",
"fitted_line2 = data.frame(income_before_tax = SFS_data$income_before_tax, wealth = regression1$coefficients[1] + regression1$coefficients[2]*SFS_data$income_before_tax)\n",
"#this is our estimated fitted line\n",
"\n",
"f <- ggplot(data = SFS_data, aes(x = income_before_tax, y = wealth)) + xlab(\"Wealth\") + ylab(\"Income before tax\")+scale_x_continuous() #defines our x and y\n",
"f <- f + geom_point() + geom_line(color = \"#070069\", data = fitted_line) + geom_line(color = \"#ff0000\", data = fitted_line2) #style preferences\n",
"\n",
"f #prints graph"
],
"id": "6da0ee33-b2d7-4f5d-9e6f-b87214bf414e"
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As you can see - there’s a very close relationship between\n",
"`after_tax_income` and `wealth`. The red line is a regression line of\n",
"wealth and after_tax_income.\n",
"\n",
"Notice as well we have negative values? Negative income and negative\n",
"wealth is weird. We will deal with this later.\n",
"\n",
"## Part 2: Simple Regressions and $t$-Tests\n",
"\n",
"Previously, we looked at the relationship between market income and\n",
"wages. However, these are both *quantitative* variables. What if we\n",
"wanted to work with a *qualitative* variable like `gender`?\n",
"\n",
"Regression models can still incorporate this kind of variable - which is\n",
"good, because (as the Census makes clear) this is the most common type\n",
"of variable in real-world data. How is this possible?\n",
"\n",
"Let’s start out with the simplest kind of qualitative variable: a\n",
"**dummy** (0 or 1) variable. Let’s use Male = $0$ and Female = $1$.\n",
"Consider the regression equation:\n",
"\n",
"$$W_i = \\beta_0 + \\beta_1 G_i + \\epsilon_i ~, \\text{where}\\ G_i \\ \\text{is Gender}$$\n",
"\n",
"Consider the conditional expectation:\n",
"\n",
"$$\n",
"E[W_i|G_i = \\text{Male}] = \\beta_0 + \\beta_1 \\cdot 1 + \\epsilon_i\n",
"$$\n",
"\n",
"$$\n",
"E[W_i|G_i = \\text{Female}] = \\beta_0 + \\beta_1 \\cdot 0 + \\epsilon_i\n",
"$$\n",
"\n",
"By the OLS reggression assumptions, we have that \\$E\\[\\_i\\|G_i\\] = 0 \\$,\n",
"so:\n",
"\n",
"$$\n",
"E[W_i|G_i = \\text{Female}] = \\beta_0 + \\beta_1\n",
"$$\n",
"\n",
"$$\n",
"E[W_i|G_i = \\text{Male}] = \\beta_0\n",
"$$\n",
"\n",
"Combining these two expressions:\n",
"\n",
"$$\n",
"\\beta_1 = E[W_i|G_i = \\text{Female}] - E[W_i|G_i = \\text{Male}] = \\beta_1-\\beta_0\n",
"$$\n",
"\n",
"What this tells us:\n",
"\n",
"1. We can include **dummy** variables in regressions just like\n",
" quantitative variables\n",
"2. The coefficients on the dummy variable have meaning in terms of the\n",
" regression model\n",
"3. The coefficients measure the (average) difference in the dependent\n",
" variable between the two levels of the dummy variable\n",
"\n",
"We can estimate this relationship of gender and wealth using R. As we\n",
"investigate the wealth gap between male and female lead households, we\n",
"might expect to see a negative sign on the coefficient - that is, if we\n",
"anticipate that female lead households will have less wealth than male\n",
"lead households."
],
"id": "1ae7f573-35b2-46ff-8d38-51a043f297de"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"regression2 <- lm(wealth ~ gender, data = SFS_data)\n",
"\n",
"summary(regression2)"
],
"id": "dc0cb389-f9dd-4ac0-9bba-0216010ecc88"
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What do you see here?\n",
"\n",
"> **Test Your Knowledge**: What is the difference in average wealth\n",
"> between male and female lead households?"
],
"id": "d15281e2-9ed1-48bf-ac7f-d5321e10720f"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# input the answer (to 1 decimal place, don't forget to add a negative sign, if relevant)\n",
"answer2 <- ??? # your answer here\n",
"\n",
"test_2()"
],
"id": "25e90ccf-3d2d-47c1-92cd-eda293549722"
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The number might seem familiar if you remember what we learned about a\n",
"$t$-test from earlier. Remember this result?"
],
"id": "67db4b1d-7ae5-486e-b499-9ca80e40aa59"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"t1 = t.test(\n",
" x = filter(SFS_data, gender == \"Male\")$wealth,\n",
" y = filter(SFS_data, gender == \"Female\")$wealth,\n",
" alternative = \"two.sided\",\n",
" mu = 0,\n",
" conf.level = 0.95)\n",
"\n",
"t1 \n",
"\n",
"t1$estimate[1] - t1$estimate[2]"
],
"id": "fe8ecf86-3a29-4337-9cba-0be50e434998"
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Look closely at this result, and the result above. What do you see? What\n",
"is the relationship here?\n",
"\n",
"This is a very important result because a dummy variable regression is\n",
"an example a two sample comparison. Why is this? Recall:\n",
"\n",
"$$\\beta_1 = E[W_i|G_i = \\text{Female}] - E[W_i|G_i = \\text{Male}]$$\n",
"\n",
"The regression coefficient of $\\beta_1$ can be interpreted as a\n",
"comparison of two means. This is exactly the same as what the $t$-test\n",
"is doing. Comparing two means by different groups - groups which are\n",
"specified by $G_i = \\text{Male}$ or $G_i = \\text{Female}$.\n",
"\n",
"In other words, another way of thinking about a regression is as a\n",
"`super` comparison of means test. However, regressions can handle\n",
"analysis using qualitative (dummy) variables as a well as quantitative\n",
"variables, which regular comparison of means tests cannot handle.\n",
"\n",
"### Multiple Levels\n",
"\n",
"Okay, but what if you have a qualitative variable that takes on *more*\n",
"than two levels? For example, the `education` variable includes four\n",
"different education classes."
],
"id": "3c543fa0-5ccc-4207-b49d-67fe2e2f5261"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"SFS_data %>%\n",
"group_by(education) %>%\n",
"summarize(number_of_observations = n())"
],
"id": "67470a61-4ac4-436a-baa8-d938b8476bed"
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this case, the idea is that you can replace a qualitative variable by\n",
"a *set of dummies*. Consider the following set of variables:\n",
"\n",
"- `d_1`: Is highest education less than high school? (Yes/No)\n",
"- `d_2`: Is highest education high school? (Yes/No)\n",
"- `d_3`: Is highest education non-university post-secondary? (Yes/No)\n",
"- `d_4`: Is highest education university? (Yes/No)\n",
"\n",
"These four dummy variables capture the same information as the\n",
"qualitative variable `education`. In other words, if we were told the\n",
"value of `education` we could discern which of these dummies were `Yes`\n",
"or `No`, and vice-versa. In fact, if wetake a closer look, we’ll notice\n",
"that we actually only need three of the four to figure our the value of\n",
"`education`. For example, if I told you that `d_4`, `d_3`, `d_2` were\n",
"all “No”, what would the value of `education` be?\n",
"\n",
"In other words, one of the dummies is redundant in helping us understand\n",
"the qualitative variable. This property is important; we usually will\n",
"omit one possible dummy to include only the minimum number of variables\n",
"needed to explain the qualitative variable in question. This omitted\n",
"dummy is called the **base level**. If we forget about this and still\n",
"add 4 dummy variables, we would be committing a dummy variable trap.\n",
"\n",
"- Which one should be the base level? It doesn’t matter, from a\n",
" technical perspective.\n",
"\n",
"> **Test Your Knowledge**: suppose you have a qualitative variable with\n",
"> $k$ distinct levels. What is the minimum number of *possible* ways to\n",
"> represent a set of dummies if you don’t want to include any redundant\n",
"> variables?\n",
"\n",
"- **A**: $k$\n",
"- **B**: $k-1$\n",
"- **C**: $k+1$\n",
"- **D**: $k^2$"
],
"id": "13024020-dbde-4564-83bc-d4991b7cfd37"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"answer2.5 <- ??? # type in your answer here \n",
"\n",
"test_2.5()"
],
"id": "402bf8c2-2db5-40b5-9610-a4caa9c6ed4a"
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In general, in R, most commands will automatically handle this process\n",
"of creating dummies from qualitative variables. As you saw with the\n",
"simple regression, R created them for you. You can also create dummies\n",
"using a variety of commands, if necessary - but in general, if you tell\n",
"R that your variables are factors, it will automatically handle the\n",
"creation of dummies properly.\n",
"\n",
"Technically, the example above which includes multiple variables is\n",
"called a **multiple regression** model, which we haven’t covered yet.\n",
"\n",
"Let’s explore regression some more, in the following series of\n",
"exercises.\n",
"\n",
"## Part 3: Weighting\n",
"\n",
"In the previous examples, we have assumed that all observations are\n",
"equally important. However, this is not always the case. Sometimes, we\n",
"may want to give more weight to some observations than others. This\n",
"could be for a wide variety of reasons, but a common one has to do with\n",
"the observations in your sample as compared to the true population\n",
"(other common reasons include data reliability and the need to correct\n",
"for selection bias).\n",
"\n",
"Consider the Survey of Financial Statistics versus the true population\n",
"of Canada. The SFS is a sample of the population, and it is possible\n",
"that the sample is not representative of the population. In this case,\n",
"we might want to weight the observations in the sample to make them more\n",
"representative of the population."
],
"id": "0e785ebe-0617-4e4e-b5c1-e4dc72c07da8"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"count(SFS_data)"
],
"id": "46069482-b2c2-4f4d-82fb-506e0a14d6a9"
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In our case, the SFS surveyed 10,141 individuals. However, the\n",
"[population of Canada in 2019 was 37.6 million\n",
"people](https://www150.statcan.gc.ca/n1/pub/11-627-m/11-627-m2019061-eng.htm).\n",
"While completing the census is a legal obligation which is typically\n",
"good at capturing practically all inhabitants, the SFS is voluntary and\n",
"may or may not represent the population well (e.g. perhaps the super\n",
"wealthy or those in destitute poverty are underrepresented, or\n",
"marginalized groups are less likely to respond to government surveys\n",
"about something as personal as your wealth).\n",
"\n",
"Whatever the case may be, we’ll want to adjust, “weight”, our sample to\n",
"make it more representative of the population (and we’ll use the census\n",
"as our benchmark). Luckily for us, Statcan is aware of this potential\n",
"issue, and has included pre-made/recommended weights for us to use!\n",
"\n",
"To see what happens when we include weights we simply add a\n",
"straightforward argument `weight = x` in the `lm()`function, where `x`\n",
"is the name of your weighting variable. (If for some reason your weights\n",
"are contained in a different dataset than the rest of your data you\n",
"could use `weight = other_dataset$x` to specify the correct weights).\n",
"\n",
"Let’s remind ourselves of what we found in `regression1`"
],
"id": "bfa2b436-6b0f-4f64-94ed-287c0336edbc"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"summary(regression1) #original non-weighted regression"
],
"id": "e14db9e4-455c-4162-a698-7c3fe48beec3"
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"How does this compare to a weighted regression?"
],
"id": "9f2aa8f6-4aff-4523-ac54-f6bdc71bcf74"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"regression3 <- lm(wealth ~ income_after_tax, data = SFS_data, weights = pweight) #the variable for weights in SFS_data is called pweight, short for population weight. \n",
"summary(regression3)"
],
"id": "1b274cc9-5b78-46a5-bc65-1643610723fb"
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"> *Think Deeper*: What is the difference in the coefficient on\n",
"> `income_after_tax` between the weighted and unweighted regressions?\n",
"> What does the change in the intercept imply?\n",
"\n",
"How about when we looked at the gender effect on wealth?"
],
"id": "93352e9f-631c-404e-b84d-6edad643a670"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"summary(regression2) #again, our original non-weighted regression"
],
"id": "fdff03f5-60c3-4e32-9bfa-593ec490d1e0"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"regression4 <- lm(wealth ~ gender, data = SFS_data, weights = pweight)\n",
"summary(regression4)"
],
"id": "85282902-a70f-47da-b896-7bfc97667670"
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"> *Think Deeper*: What is the difference in the coefficient on `gender`\n",
"> between the weighted and unweighted regressions? What does the change\n",
"> in the intercept imply? Did the unweighted regression overestimate or\n",
"> underestimate the gender-wealth gap?\n",
"\n",
"In general, weighting is a powerful tool that can help us correct for\n",
"biases in our data. However, it is important to remember that weighting\n",
"is not a panacea. It can only correct for biases that we know about and\n",
"can measure. If there are biases in our data that we are unaware of or\n",
"we don’t know what to benchmark against, weighting will not help us.\n",
"\n",
"## Part 4: Exercises\n",
"\n",
"### Activity 1\n",
"\n",
"Last week, we briefly explored the idea of the wealth gap and explored\n",
"the idea that it could be caused by some income related factors. We can\n",
"now examine this issue directly using regressions. Run a regression with\n",
"\\* before tax income \\* on male and female lead households.\n",
"\n",
"Tested objects: `regm` (the regression for males).Tested\n",
"objects: `regf` (the regression for females)."
],
"id": "a0cc7bf5-d53a-441c-b076-2c29e8539f6e"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Quiz 1\n",
"\n",
"# Regression for males\n",
"regm <- lm(??? ~ income_before_tax, filter(SFS_data, ??? == \"Male\")) \n",
"# Replace \"...\" with the appropriate variables \n",
"#remember answers are case sensitive!\n",
"\n",
"# Quiz 2\n",
"# Regression for females\n",
"regf <- lm(??? ~ income_before_tax, data = filter(SFS_data, ??? == \"Female\")) \n",
"#remember answers are case sensitive!\n",
"\n",
"summary(regm) # Allow us to view regm's coefficient estimates\n",
"summary(regf) # Same as above, but for regf\n",
"\n",
"test_3() # Quiz1\n",
"test_4() # Quiz2"
],
"id": "131abfea-f918-4bdc-bbff-3926f71b6ba5"
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Short Answer 1\n",
"\n",
"**Prompt:** How do we interpret the coefficient estimate on `income` in\n",
"each of these regressions?\n",
"\n",
"Answer in red here!"
],
"id": "414e96f5-8c51-49c3-930d-9fb7baf805ea"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"answer_1 <- #fill in your short answer"
],
"id": "5d0943e7-9c07-4b2e-b0c8-0af3571a154c"
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Activity 2\n",
"\n",
"We might think that income inequality between females and males might\n",
"depend on the educational gaps between these two groups. In this\n",
"activity, we will explore how the income gap varies by education. First,\n",
"let’s see the factor levels of the `education`:"
],
"id": "473536d8-28a4-4471-a172-67b2b9129651"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"levels(SFS_data$education) # Run this"
],
"id": "bb23321c-9831-44cd-9a93-0a40cea98b58"
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As we can see, there are a few education groups in this dataframe. Let’s\n",
"estimate the income gap (with no controls) for each of the four groups\n",
"separately:\n",
"\n",
"- Less than high school\n",
"- High school\n",
"- Non-university post-secondary\n",
"- University\n",
"\n",
"Tested objects: `rege2` (High School), `rege4` (University)\n",
"\n",
"Notice we don’t need to do 4 regressions we could just do three."
],
"id": "455fd550-fb61-4fd3-a5e8-a6e1ef5aa4f7"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#reg1 is a regression performed on people, with a less than high scchool education\n",
"reg1 <- lm(??? ~ ???, data = filter(SFS_data, education == \"Less than high school\")) #what should replace the ...\n",
"#reg2 is the same as rege1,but we are looking at people with a high school education\n",
"reg2 <- lm(??? ~ ???, data = filter(SFS_data, education == \"High school\")) #fill in the blanks\n",
"\n",
"reg3 <- lm(??? ~ ???, data = filter(SFS_data, education == \"Non-university post-secondary\")) #remember answers are case sensitive!\n",
"\n",
"reg4 <- lm(??? ~ ???, data = filter(SFS_data,education == \"University\"))\n",
"\n",
"# store the summaries (but don't show them! too many!)\n",
"sum20 <- summary(reg1)\n",
"sum30 <- summary(reg2)\n",
"sum40 <- summary(reg3)\n",
"sum50 <- summary(reg4)\n",
"\n",
"test_9() \n",
"test_10() \n",
"test_11() \n",
"test_12() "
],
"id": "55f77eb9-09f2-4e87-8c10-35f461d56ff4"
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The code below will tabulate a brief summary of each regression:"
],
"id": "4bf533bb-68b7-45ec-b131-36fa0b3c07b3"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# just run me. You don't need to edit this\n",
"\n",
"Educ_Group <- c(\"Less than high school\", \"High School\", \"Non-university post-secondary\", \"University\") #defines column 1\n",
"Income_Gap <- c(reg1$coefficients[2], reg2$coefficients[2], reg3$coefficients[2], reg4$coefficients[2]) #defines column 2\n",
"Std._Error <- c(sum20$coefficients[2,2], sum30$coefficients[2,2], sum40$coefficients[2,2], sum50$coefficients[2,2]) #defines column 3\n",
"t_Value <- c(sum20$coefficients[2,3], sum30$coefficients[2,3], sum40$coefficients[2,3], sum50$coefficients[2,3]) #defines column 4\n",
"p_Value <- c(sum20$coefficients[2,4], sum30$coefficients[2,4], sum40$coefficients[2,4], sum50$coefficients[2,4]) #defines column 5\n",
"\n",
"tibble(Educ_Group, Income_Gap, Std._Error, t_Value, p_Value) #it's like a table but a tibble"
],
"id": "dbe711fb-1b12-41a4-bf74-009f62083629"
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Short Answer 2\n",
"\n",
"**Prompt**: What happens to the income gap as we move across eduction\n",
"groups? What might explain these changes? (hint: think back to module\n",
"1!)\n",
"\n",
"Answer in red here!"
],
"id": "36315a8c-1809-487f-8346-d33f38df0adc"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"answer_3 <- #fill in your short answer"
],
"id": "f123143e-20f6-425b-91ff-05daf38b4a3c"
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Activity 3\n",
"\n",
"As we observed in last week’s worksheet, the income gap could differ by\n",
"education level. Since there are many education categories, however, we\n",
"may not want to examine this by running a regression for each education\n",
"level separately.\n",
"\n",
"Instead, we could run a single regression and add education level as a\n",
"second regressor, $E_i$:\n",
"\n",
"$$I_i = \\beta_0 + \\beta_1 G_i + \\beta_2 E_i + \\epsilon_i$$\n",
"\n",
"This is actually a **multiple regression**, which we will learn about\n",
"next week - but from the point of the this lesson, the idea is that it\n",
"is “run” in R essentially in the same way as a simple regression.\n",
"Estimate the regression model above without $E_i$, then re-estimate the\n",
"model with $E_i$ added. **USE INCOME BEFORE TAX**.\n",
"\n",
"Tested objects: `reg2A` (regression without controls), `reg2B`\n",
"(regression with controls)."
],
"id": "a6bb77f2-ee7e-4566-b872-4270413b308b"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Simple regression (just gender)\n",
"reg2A <- lm(income_before_tax ~ gender, data = SFS_data) # this one works already\n",
"\n",
"# Regression with controls\n",
"reg2B <- lm(income_before_tax ~ ??? + education, data = SFS_data) # replace the ...\n",
"\n",
"summary(reg2A)\n",
"summary(reg2B)\n",
"#this will look ugly; try to look carefully at the output\n",
"\n",
"test_7()\n",
"test_8() "
],
"id": "db270098-d040-492d-95de-0dc954fc7522"
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Short Answer 3\n",
"\n",
"Prompt: Compare the estimated income gap with and without $E_i$ in the\n",
"regression. What happens to the gap when we add $E_i$?\n",
"\n",
"Answer in red here!"
],
"id": "5aaf3550-5c3d-48a8-a092-90a5359b72ed"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"answer_4 <- #fill in your short answer"
],
"id": "1dc8dfc1-1830-4c26-b36c-356c0317dee8"
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Theoretical Activity 1\n",
"\n",
"When we deal with large quantitative variables, we often take the\n",
"natural log of it:"
],
"id": "713a3c72-a034-401b-a5db-247334b860a4"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"W = log(SFS_data$wealth[SFS_data$wealth>0]) "
],
"id": "a69818a5-0bf2-4d89-afd0-005e41d6e43e"
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You may recall that the derivative of the log of a variable is\n",
"approximately equal to percentage change in the variables:\n",
"\n",
"$$\\frac{dln(x)}{dx} \\approx \\frac{\\Delta x}{x}$$\n",
"\n",
"Thus, when we find the marginal effect of some continuous regressor\n",
"$X_i$ (say, `income`):\n",
"\n",
"$$ln(W_i) = \\beta_0 + \\beta_1 I_i + \\epsilon_i \\implies \\frac{\\Delta W_i}{W_i} \\approx \\beta_1 \\Delta I_{i}$$\n",
"\n",
"This allows us to interpret the changes in a continuous variable as\n",
"associated with a percentage change in wealth; for instance, if we\n",
"estimate a coefficient of $0.02$ on `income_before_tax`, we say that\n",
"when a family’s income before tax increases by 1 CAD, the corresponding\n",
"wealth increases by 2 percent on average.\n",
"\n",
"Notice as well we are now talking about percent changes, rather than\n",
"units.\n",
"\n",
"Let’s generate two variables that take the natural log of the wealth\n",
"and market income from the `SFS_data` dataframe (hint: use a\n",
"technique that we introduced last week). Then, estimate the effect of\n",
"logarithmic market income on logarithmic wealth.\n",
"\n",
"Tested Objects: `lnreg`"
],
"id": "95f00dd4-02b9-4367-89ab-8f9029b83fc5"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#Generate log wage variable\n",
"SFS_data <- SFS_data %>%\n",
" mutate(lnincome = log(SFS_data$income_before_tax)) %>% # what goes here?\n",
" mutate(lnwealth = log(SFS_data$wealth)) # what goes here?"
],
"id": "327d9ee5-f317-4232-a470-f322a7c163f4"
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Notice warning message “NaNs produced”. NaN means “Not a Number”. This\n",
"happens because we had negative income and negative wealth. No matter\n",
"how low our incomes are, the more we work, wealth and income should\n",
"increase."
],
"id": "a592095e-bee1-485d-b9a4-77662c353566"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# fix NANs\n",
"SFS_data_logged <- SFS_data %>%\n",
" filter(income_before_tax>0) %>% #removes negative values\n",
" filter(wealth>0) #removes negative values\n",
" \n",
"# Log Regression \n",
"lnreg <- lm(lnwealth ~ ???, data = SFS_data_logged) #the new and improved regression\n",
"\n",
"\n",
"summary(lnreg)\n",
"\n",
"test_5() #Quiz7"
],
"id": "d74a8b56-5d04-41ec-a18f-89482f3933e6"
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Short Answer 4\n",
"\n",
"Prompt: How do we interpret each of these estimates? (Hint: what does a\n",
"1-unit change in the explanatory variable mean here?)\n",
"\n",
" Answer here in red"
],
"id": "609a9114-3ef5-4b94-8623-7661315740ef"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"answer_5 <- #fill in your short answer"
],
"id": "e1817206-90fd-45ac-91f6-a8f5e858e3e6"
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Optional: Bonus Activity 4\n",
"\n",
"You have learned about a linear regression model of income; however,\n",
"income often follows a Pareto distribution. For now, using a linear\n",
"approximation to find the wage gap is fine. We may want to know stuff\n",
"about the underlying distribution of income in male and female lead\n",
"households, however. Here’s the PDF of pareto distribution:\n",
"\n",
"$$f(x) = {\\displaystyle {\\frac {\\alpha x_{\\mathrm {m} }^{\\alpha }}{x^{\\alpha +1}}}} $$\n",
"\n",
"Ok, now with regression remember we said that we estimate the parameter\n",
"given the data. To do this we said you could use Calcus or methods other\n",
"than OLS. Here the probability of the data can be approximated by\n",
"assuming independence between each $x_i$. If we do this, the probability\n",
"of the data is given by:\n",
"\n",
"$$\\Pi_{i=1}^n f(x)$$\n",
"\n",
"Now we can just make a function in r and optimize over it which performs\n",
"essentially the same operation as a linear regression."
],
"id": "003c327f-901b-4db9-b3ad-6f8e56bc75d4"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"x=filter(SFS_data,gender=='Female')\n",
"x <- filter(x, is.numeric(income_before_tax))\n",
"x <- x$income_before_tax"
],
"id": "1b51c18a-2c2b-48e2-ac76-3c87818f78f8"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"calc <- function (x){\n",
" q=0\n",
"for (i in x){\n",
" if (i >0){\n",
" a= log(i[1]) }\n",
" if (is.numeric(a)==TRUE){\n",
" q=q+a }\n",
" }\n",
"return (q)\n",
"}"
],
"id": "3e26b050-a98f-4f23-a98a-0ef6a1673b0f"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"calc(x)"
],
"id": "e1328300-a329-4c8f-8ffa-1eb51e6696b4"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ell <- function(a,q,xm,n) { # we use the log function of the pareto distrubtion instead\n",
" d=(n*log(a))\n",
" b=(-1)*(a+1)*q\n",
" c=a*log(xm)*n \n",
" return (d+b+c)\n",
"}"
],
"id": "ec5fd128-3355-4781-87de-50726514f8ab"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"a = optimize(ell,c(2,50),maximum=TRUE,q=43074.1853103325,xm=40000,n=length(x))\n",
"a\n",
"a_women=a$maximum "
],
"id": "2ebdca4a-f42e-4a61-8990-6c59b57717bc"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"y=filter(SFS_data,gender=='Male')\n",
"y <- filter(y, is.numeric(income_before_tax))\n",
"y <- y$income_before_tax\n",
"a_men = optimize(ell,c(2,1000),maximum=TRUE,q=calc(y),xm=65000,n=length(y))\n",
"a_men = a_men$maximum\n",
"a_men"
],
"id": "278adc54-cf93-49e0-8c57-cf3db587957e"
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The theoretical mean of the Pareto distribution is,\n",
"\n",
"$$ \\frac{\\alpha x_m}{\\alpha -1} $$ Can you calculate the expected income\n",
"gap with the Pareto distribution assumption?"
],
"id": "a4e4f1d7-a583-488d-998b-850b782e02c6"
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"xmw=40000\n",
"xmm=65000\n",
"income_gap =((a_women* xmw )/ (a_women-1)) - ((a_men* xmm )/ (a_men-1))\n",
"income_gap #note we set xm ourselves (I did this by playing around with xm, and doing a bit of research) see if you can get a better xm."
],
"id": "51df47b9-ae0f-4716-83ed-282f44254f88"
}
],
"nbformat": 4,
"nbformat_minor": 5,
"metadata": {
"kernelspec": {
"name": "ir",
"display_name": "R",
"language": "r"
}
}
}