One X Variable

OLS in R

In the last set of notes, we derived the OLS estimator, calculated statistics for hypothesis tests, and measures for goodness of fit. In this set of notes, we are going to talk more about application.

First, to estimate a linear regression model, we will use the lm() function.

• formula: Formulas in R take on the following form: y ~ x1 + x2. In the case for GPAs and SAT scores, we would write: sat_sum ~ hs_gpa.
• 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. For example, maybe we want to fit two lines: one for men and one for women. We would use subset = df$sex == 1 to estimate the model using only data from female students. Another reason to use this would be to eliminate outliers, etc.

OLS in R

We are going to practice on two sets of data. First, we are going to keep examining the GPA data (data; documentation). Second, we are going to explore some data on real estate transactions in Ames, Iowa (data; documentation). Let’s read both sets of data into R:

df <- read.csv("https://vincentarelbundock.github.io/Rdatasets/csv/openintro/satgpa.csv")
ames <- read.csv("https://vincentarelbundock.github.io/Rdatasets/csv/openintro/ames.csv")

There are a lot of columns in the ames data. We are only interested in a few (for now), so I am going to subset the data to keep only the columns we’ll use.

keep_columns <- c("price", "area", "Bedroom.AbvGr", "Year.Built", "Overall.Cond")
ames <- ames[,keep_columns]

OLS in R

Using the summary() function generates some important information about our regression and the individual coefficients.

summary(reg_sat_gpa)

Call:
lm(formula = sat_sum ~ hs_gpa, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-42.055  -8.669  -0.351   8.803  34.240 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  66.9885     2.4441   27.41   <2e-16 ***
hs_gpa       11.3632     0.7535   15.08   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 12.9 on 998 degrees of freedom
Multiple R-squared:  0.1856,    Adjusted R-squared:  0.1848 
F-statistic: 227.4 on 1 and 998 DF,  p-value: < 2.2e-16

In the table portion, each coefficient has an Estimate, Std. Errors, t value, and Pr(>|t|) In addition, you can find the Multiple R-Squared value, which is the measure for goodness of fit.

OLS in R

reg_sat_gpa_m <- lm(sat_sum ~ hs_gpa,
                    data = df, subset = df$sex == 2)
summary(reg_sat_gpa_m)

Call:
lm(formula = sat_sum ~ hs_gpa, data = df, subset = df$sex == 
    2)

Residuals:
    Min      1Q  Median      3Q     Max 
-31.997  -8.305  -0.455   8.574  33.287 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   55.851      3.579   15.61   <2e-16 ***
hs_gpa        13.715      1.079   12.71   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 12.33 on 482 degrees of freedom
Multiple R-squared:  0.2512,    Adjusted R-squared:  0.2496 
F-statistic: 161.7 on 1 and 482 DF,  p-value: < 2.2e-16

reg_sat_gpa_f <- lm(sat_sum ~ hs_gpa,
                    data = df, subset = df$sex == 1)
summary(reg_sat_gpa_f)

Call:
lm(formula = sat_sum ~ hs_gpa, data = df, subset = df$sex == 
    1)

Residuals:
    Min      1Q  Median      3Q     Max 
-45.497  -8.323   0.153   8.877  31.153 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   71.280      3.183   22.39   <2e-16 ***
hs_gpa        11.019      1.003   10.98   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 12.55 on 514 degrees of freedom
Multiple R-squared:  0.1901,    Adjusted R-squared:  0.1885 
F-statistic: 120.6 on 1 and 514 DF,  p-value: < 2.2e-16

OLS in R

Take a look at the Pr(>|t|) column. This contains a p-value corresponding to the hypothesis test that \(\beta_i \neq 0\). The value <2e-16 means that the p-value is smaller than 0.0000000000000002. R doesn’t report p-values that are smaller than this for brevity. In addition, R includes asterisks next to p-values that are statistically significant.

• Nothing next to the p-value indicates \(0.1 < p \leq 1\).
• . next to the p-value indicates \(0.05 < p \leq 0.1\)
• * next to the p-value indicates \(0.01 < p \leq 0.05\)
• ** next to the p-value indicates \(0.001 < p \leq 0.01\)
• *** next to the p-value indicates \(0 < p \leq 0.001\)

This is helpful when looking to see if your coefficients are significant at different levels.

OLS in R

How can we interpret these results? Let’s start with \(\beta_0\). For a male student with a GPA of 0, the expected SAT score is 55.85. For a female student with a GPA of 0, the expected SAT score is 71.28. In addition, for any GPA, we can compute the expected SAT score. We can write a function to do this for us.

OLS in R

Thankfully, R has a pre-defined function called predict(). This function accepts two arguments: a lm object and newdata. newdata must have the same column names as variables in the lm object.

# Create a new data.frame of GPAs from 2 to 4 by 0.25 increments.
pred_df <- data.frame(hs_gpa = seq(2, 4, by = 0.25))
pred_df$sat_m <- predict(reg_sat_gpa_m, newdata = pred_df)
pred_df$sat_f <- predict(reg_sat_gpa_f, newdata = pred_df)
pred_df

  hs_gpa     sat_m     sat_f
1   2.00  83.28204  93.31803
2   2.25  86.71088  96.07282
3   2.50  90.13971  98.82762
4   2.75  93.56855 101.58242
5   3.00  96.99738 104.33721
6   3.25 100.42622 107.09201
7   3.50 103.85505 109.84680
8   3.75 107.28389 112.60160
9   4.00 110.71272 115.35639

It seems that female students consistently outperform male students on the SAT, according to these data.

OLS in R

Another helpful way to look at this would be to plot both lines together with the underlying data.

In this plot, it’s very clear that the biggest difference between male and female students is when GPAs are relatively low. The gap then decreases as GPAs get larger.

OLS in R

library("modelsummary")
options("modelsummary_factory_default" = "kableExtra")

regz <- list(`All Students` = reg_sat_gpa,
             `Male` = reg_sat_gpa_m,
             `Female` = reg_sat_gpa_f)
coefz <- c("hs_gpa" = "High School GPA",
           "(Intercept)" = "Constant")
gofz <- c("nobs", "r.squared")
modelsummary(regz,
             title = "Effect of GPA on SAT Scores",
             estimate = "{estimate}{stars}",
             coef_map = coefz,
             gof_map = gofz)

Ames, Iowa

Location, Location, Location

Real estate agents will tell you that the most important factor in determining the price of a property is location. Whether or not this is true, there are obviously other things that contriibute to a property’s value. For example, square footage, number of bedrooms, age, condition, etc. all play a role in determining a home’s value.

We are going to estimate a bunch of regressions using the ames data to demonstrate economic interpretations of OLS coefficients.

Summary Statistics

Let’s rename some columns and take a peak at our data.

colnames(ames) <- c("price", "sqft", "bedrooms", "yr_built", "condition")
View(head(ames, 100)) # View() in RStudio will allow you to see the data in a tab.

Summary Statistics

Interpreting the summary statistics table:

• price: This appears to be measured in dollars with a lot of variation (e.g. about 1000 unique values out of about 3000 observations). The average sale price is $180,000 with a standard deviation of $80,000. The standard deviation is quite high relative to the mean, meaning the distribution is very wide. It’s likely there are a lot of outliers in the right tail, which is typical of housing price data.
• sqft: This exhibits very similar characteristics as the price variable. The mean is 1,500 with a long right tail, which means there are a few very large homes.
• bedrooms: This variable only takes on 8 unique values, meaning there is not much variation. The average home has about three bedrooms, which matches my prior expectations of average houses. Once again, there is likely a long right tail, evidenced by the maximum of eight bedrooms.

Summary Statistics

Interpreting the summary statistics table (continued):

• yr_built: This variable has a maximum of 2010, which would be brand new construction, but a mean of about 1970. The minimum value is 1870, which would be quite an old home.
• condition: This variable has a mean and median of about 5. If 5 means “average condition”, then this would make sense. However, any quantitative interpretation of this variable is effectively meaningless. What does it mean for a home to improve from a 6 to a 7? Is this the same as moving from a 3 to a 4? This is an ordinal variable, and should not be considered in a linear regression. We will explore this variable nonetheless.

Summary Statistics

It’s always helpful to create distribution plots for your main variables. Sometimes, this step can help inform you when making modeling decisions. For example, when variables have long right tails, you’ll often see people use a log transformation.

par(mfrow = c(1, 2))
plot(table(round(ames$price/10000)) / nrow(ames),
     xlab = "Price", ylab = "Rel. Freq.")
# plot(table(round(ames$sqft/100)))
plot(table(round(log(ames$price), 1)) / nrow(ames),
     xlab = "log(Price)", ylab = "Rel. Freq.")
par(mfrow = c(1, 1))

Log Models

What would the log transformation do to our models? Let’s write down four models, visualize the data, and examine estimated parameters:

Of course, log() is a non-linear transformation. However, that does not go against our definition of a linear model. For a model to be linear, all variables, no matter their transformation, must enter into the model linearly. Therefore, we can have anything of the form \(Y = \alpha + \beta X + \epsilon\), even if \(X\) is non-linear. In fact, we could estimate a model like: \(Y = \alpha + \beta^X + \epsilon\) following a (relatively) simple log transformation.