Introduction

Just like how the first half of the module was a quick (re-)primer on programming in R, this will be a short refresher on econometrics. To start, a simple model that uses some variable \(X\) to explain some outcome \(Y\) can be written as follows:

\[Y = \beta_0 + \beta_1 X + \epsilon\] If there are multiple variables that explain variation in \(Y\), we can rewrite this equation as follows. Econometrics is concerned with estimating, interpreting, and testing the \(\beta\)’s in this equation:

\[Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \epsilon\]

In this model, we would say that a one unit change in \(X_1\) is associated with a \(\beta_1\) unit change in \(Y\), holding \(X_2\) constant. For this interpretation to make sense, a one unit change in \(X_1\) must be tangible and concrete. For example, a one dollar increase in income would make sense, since increasing income from $50,000 to $50,001 is the same as increasing from $60,000 to $60,001. On the other hand, increasing by one unit on a likert scale of happiness, or some ranking system, does not make sense in the context of linear regression. Is increasing happiness from a 6 to a 7 the same as increasing from 7 to 8? Do 7s even mean the same thing for everyone? Clearly, interpreting the coefficient if this were the variable would become difficult.

In this course, we’ll often work with binary variables. These variables are binary in the sense that they can take on values of either one or zero. These are used to represent factors that are dichotomous such as sex, whether some is insured, or if they’ve attended some training. Back to our original model, if \(X_1\) was binary, we would interpret \(\beta_1\) as the difference in expected value of \(Y\) for when \(X_1 = 1\) compared to when \(X_1 = 0\), holding other factors constant. This, if you think about it, still represents a one unit increase in \(X_1\). You can also think of binary variables as on/off switches. For example, \(X_1\) could represent whether some policy has been enacted since this would otherwise be difficult to quantify.

For a real life example of this, consider the gender wage gap. Suppose \(X_1\) measured someone’s years of experience, \(X_2\) measured their sex (1 for men, 0 for women), and \(Y\) was wages. We would interpret \(\beta_1\) as the additional wages earned for one more year of experience for women. To calculate the additional wages earned for one more year of experience for men, we would need to add \(\beta_1\) and \(\beta_2\). \(\beta_2\) would then represent the additional wages earned for men relative to women for an additional year of experience. If this were zero, then the two groups would get similar returns for their experience. If this were positive, it would signal to us that men tend to get larger raises than women.

So, how do we estimate parameters (\(\beta\)s) in linear models? The most famous algorithm is called Ordinary Least Squares (OLS). How does it work? Intuitively, OLS finds the line of best fit for some set of data. However, this analogy only works with a single explanatory variable. If there are two variables, OLS finds the plane of best fit. If there are three variables, the analogy breaks down a bit, but hopefully you get the point. Intuitively, this means that OLS can determine the ceteris paribus relationship between each \(X\) variable and the outcome \(Y\). This is appealing to economists since we are always interested in how some thing \(X_1\) will effect \(Y\), holding everything else constant. This is precisely how OLS works.

Above, I claimed that OLS can pick out the relationship between \(X_1\) and \(Y\), while holding everything else constant. This statement needs a small disclaimer: OLS can only hold constant the factors that are included in the estimation. In other words, if OLS doesn’t know about some factor \(X_2\), it’s unable to partial out its effects. This is fine when \(X_1\) and \(X_2\) are uncorrelated or independent, but it becomes an issue as the relationship between \(X_1\) and \(X_2\) becomes stronger. The example that I think does a good job of demonstrating this uses housing prices, bedrooms, and square footage. We will talk through the example here and then code it up below.

According to rockethomes.com, there were about 1,200 homes for sale in November 2023. The median sale price for homes was about $320,000. The median sale price, by number of bedrooms, was as follows:

What’s missing from this table is the size of each house. Houses with \(n+1\) bedrooms are (generally) bigger than houses with \(n\) bedrooms. So, is the relationship between bedrooms and price actually driven by bedrooms, or is it driven by the property’s square footage? In other words, what happens when we “control for”, “condition on”, or “partial out” square footage? Let’s write down our empirical model of sale price for property \(i\) below:

\[P_i = \beta_0 + \beta_1 SF_i + \beta_2 BR_i + \epsilon\] where \(P\) represents price, \(SF\) represents square footage, and \(BR\) represents bedrooms. How can we interpret each parameter?

• \(\beta_0\) is difficult to interpret since not only do properties with zero square feet and zero bedrooms not typically transact on the market, they do not really exist. For this reason, we would avoid thinking too hard about how to interpret the intercept.
• \(\beta_1\) is the effect of increasing the square footage of a home, holding the number of bedrooms constant. In words, if we increase the size of a property by one square foot, we expect the property to sell for an additional $\(\beta_1\).
• \(\beta_2\) is the effect of increasing the number of bedrooms in the home, holding the size constant.

Omitted Variable Bias is defined as when some variable not in the regression is driving, or contributing to, an estimated relationship. We need to be very careful to consider omitted variable bias when we interpret our coefficients. Unless we can be sure that we’ve controlled for all potential confounding variables (i.e. variables that would bias our results if they were omitted from the model), we must hesitate before we can claim causality. Much of this course can be boiled down to different ways of addressing concerns about omitted variable bias. You should keep this in the back of your mind as you move through the course.

The first function we’ll use to estimate regressions in R is lm(), which stands for “linear model”. The syntax for this function’s arguments are as follows:

• formula: Formulas in R take on the following form: y ~ x1 + x2.
• data: This is the data you plan to use. This tells R where to look for the variables in your formula.
• subset: An optional vector that specifies a subset of observations used to fit the line. This can be either a logical vector or a vector of indices.

Coefficient estimates and standard errors from these three models can be seen in the table above. In the first column, we find a coefficient of 111.694 on square footage. This suggests that for every square foot increase in the size of a property, the price increases by $111.69. In the second column, our estimate for bedrooms suggests that every additional bedroom is “worth” $13889.5. Let’s think hard about this for a second – that’s not really what this coefficient means. This regression isn’t controlling for anything, so this coefficient is picking up everything associated with extra bedrooms. Remember, this usually means larger homes, and larger homes mean garages and backyards and pools. The third column is trying to address this by holding property size constant. Now, we find a negative effect of additional bedrooms. Adding an extra bedroom, while holding square footage constant, either makes every other bedroom smaller or takes away from other rooms in the home.

Hopefully, this demonstration has illustrated the potential perils of omitted variable bias. It should also be noted that you can “sign” the bias associated with omitted variables if you can make certain assumptions about the relationship between the omitted variable, the explanatory variable and the outcome variable. If you have a model \(Y = \delta_0 + \delta_1\times A + \epsilon\) and \(B\) is your omitted variable, the bias in \(\delta_1\) will be:

Positive bias would indicate that the magnitude of the coefficient is pushed upwards, in the positive direction, artificially. To think through this, we can look at column (2) in the table above. We find a coefficient of 13889.5, though we believe this to be effected by omitted variable bias. Since we think square footage and number of bedrooms are positively correlated, and that the effect of square footage on price is positive, we would say that this coefficient is biased upwards. Once we “correct” for this bias by including square footage into the model in column (3), we see the coefficient for bedrooms plummet. Therefore, as we predicted, the initial, unconditional coefficient estimate was too positive, or biased upwards, due to square footage being omitted from the model.

So, what are the differences? Actually, there are none! The output is very similar and the coefficients are identical. So what’s the catch? Let me introduce some additional capabilities, arguments, and functions.

• Multiple Outcomes: There are many cases where you might want to estimate the same model for multiple outcomes. For example, maybe one model where you use price as an outcome, and another where you use log(price). To do this, put c(Y1, Y2) as the left side of your formula. Specifically, c(price, log(price)) ~ sqft + bedrooms.
• Stepwise Estimation (sw(), sw0()): Similarly, you might want to estimate different specifications for each outcome. Putting control variables into the sw() function will estimate the model multiple times with different sets of controls. For example, we can estimate the three models before in a single feols() call by using sw(sqft, bedrooms, sqft + bedrooms) on the right side of the formula. sw0() will do the same thing, except it is a shortcut for sw(0, sqft, bedrooms, sqft + bedrooms). Keep whatever variables you’d like in every version of the model outside of the sw() function.
• Subsampling: (split, fsplit): This argument allows one to estimate the model(s) on different subsets of data. For example, suppose these housing sales data were for multiple cities. If we wanted a unique regression for each city, we could add split = ~city to the function call. This will estimate a model using only data from each unique city. The difference between split and fsplit is that fsplit will estimate the model on the full sample in addition to each subsample.

Another major advantage of using fixest is the convenience it provides when estimating fixed effects. Fixed effects control for unobservable factors that cause differences across groups. For example, let’s suppose we had housing data for Manhattan, New York and Ames, Iowa. Homes in Manhattan are very expensive and very small. Homes in Iowa are relatively bigger and cheaper. Therefore, if you try to estimate a relationship between price and square footage, you might end up finding that smaller homes sell for higher prices. This is, in some ways, another manifestation of omitted variable bias! We would then be able to control for the differences in price levels and property sizes by estimating New York and Iowa fixed effects. See below for a demonstration.

The titular function, modelsummary(), can accept a lot of different arguments, but for simplicity, we will just discuss a few:

• models: A named list of models we want to include in the table.
• title: The title of our table.
• estimate: This is a bit complicated, but we are going to set this equal to "{estimate}{stars}". This way, the table will display coefficients, standard errors, and stars to denote statistical significance. This is the norm in economics. If you do not like the stars, you can remove them by omitting the estimate argument in the modelsummary function call.
• coef_map: A named vector of variable names. This will rename the variable names from what they are labeled as in the data to something more human readable.
• gof_map: A vector of desired goodness of fit measures for each model. We are going to stick with "nobs" and "r.squared" for the time being.

The last topic we will discuss in these notes is statistical inference. Of course, we are interested in estimating models to generate and interpret \(\beta\)s, but we are also interested in testing hypotheses. You might have noticed that both summary() and modelsummary() generate standard errors. Standard errors are a measure of uncertainty for coefficients. Smaller standard errors indicate higher levels of precision, which is obviously desired. However, sometimes standard errors can be misleadingly small (or large), so we need to be sure we take care in addressing this. We will first talk about what we do with standard errors, and then we will discuss how we might want to adjust our standard errors.

Instead of testing whether zero is inside or outside of a set confidence interval, we can also calculate \(p\)-values. These measure the probability of observing as large of a \(\widehat{\beta}\) given the true \(\beta = 0\), or some other hypothesized value. Once your sample size is over a couple hundred, and you’re testing \(\widehat{\beta}\) against 0, you can simply divide \(\widehat{\beta}\) by its standard error and check out the result. If this result is equal to 1.96 in absolute value, then 0 is lying right on the border of the 95% confidence interval. If the result is larger, then you can reject the hypothesis that \(\beta = 0\). \(t\)-statistics of 1.96 correspond to 95% confidence intervals which then correspond to \(p\)-values of 0.05. As the \(t\)-statistic increases, the \(p\)-value will decrease. R will calculate \(p\)-values for you, and a typical threshold is 0.05.

For a more thorough discussion of hypothesis testing, visit this link.

When we use feols() without fixed effects, the default is "iid" standard errors, which are also called analytical standard errors. Again, heteroskedasticity robust standard errors are preferred to analytical standard errors. However, when you use feols() with fixed effects, the default is clustered standard errors at the level of the first fixed effect. Clustering is outside of the scope of the course, but the idea is that standard errors should be clustered at the level of treatment. For example, if we had some policy changing at the state level, and we could observe multiple individuals within each state across time, we would cluster our standard errors at the state level. In other words, we would not use vcov = "hetero", we would just make sure that the state fixed effect comes first in our set of fixed effects, or we could use cluster = ~state in feols().