More on formats: Numeric vs. double

Someone asked about double versus numeric. It turns out that numeric is a more general class/mode. Numeric objects come in different modes. Specifically, numeric objects can be either double-precision or integer (single-precision is not really an option in R, unless you are calling C or Fortran at a lower level).

In practice:

# Does as.numeric() create integers or doubles?
is.double(as.numeric(1))
## [1] TRUE
is.integer(as.numeric(1))
## [1] FALSE
# Are integers and doubles numeric?
is.numeric(as.double(1))
## [1] TRUE
is.numeric(as.integer(1))
## [1] TRUE

Follow up: Vectorized operations

I want to point out that I probably did not give vectors a fair shake. While seemingly simple, R allows you to do a lot of things with vectors that might be much more difficult in other languages. Specifically, R allows you to apply (most) functions to vectors, rather than individual elements of a vector.

For an underwhelming example, if you want to square each number in a vector vec in R, you can simply write vec^2. The alternative that many languages use requires iterating through the vector and squaring the individual elements while simultaneously storing the results.

# Define the vector
vec <- 1:5
# Square the elements of the vector
vec2 <- vec^2
# Look at the result
vec2
## [1]  1  4  9 16 25

Speeding up knitr

You may have noticed that knitr can take a little while to compile your document (.Rnw or .Rmd) once you’ve written a bunch of R code. This slowness is often due to the fact that R is re-running each of your code chunks every time you compile your file. To speed things up, you can tell R to cache the code chunk either—store the code chunk in memory until something inside the chunk changes. You have two options here:

Option 1: At the beginning of your document, add a code chunk that has only the following two¹ lines of code:

library(knitr)
opts_chunk$set(cache = T)

These two lines load the knitr package and then tell R to default to caching the chunks.

Option 2: Use the cache = true option within any/all code chunks in your document.²

For more on knitr chunks and options, check out the knitr website or my handy summary of LaTeX and knitr.

Missing data (NA) in R

As you dive into your problem set, you may notice that there are some missing values. In R, missing values typically take the form of NA (R remains agnostic as to why the datum is missing). NAs are very special and can give you some challenges when coding.

First, what class is NA?

# Class of NA
class(NA)

## [1] "logical"

# Class of NA from a vector of numbers
class(c(1, NA)[2])

## [1] "numeric"

# Class of NA from a vector of characters
class(c("hi", NA)[2])

## [1] "character"

# Class of NA from a vector of logicals
class(c(T, NA)[2])

## [1] "logical"

Sort of makes sense, right? Any class could be missing a value. Also notice the funny behavior of NA with some other functions:

# Addition
2 + NA

## [1] NA

# Pasting
paste("Test", NA)

## [1] "Test NA"

Luckily, R has a special function, just for NA, that tests whether an object is NA. What’s it’s name? is.na()

# Demonstrate is.na()
is.na(NA)

## [1] TRUE

is.na(1)

## [1] FALSE

is.na(T)

## [1] FALSE

In addition to simply missing data, R treats not-a-number (NaN) as NA:

# What is NaN?
0 / 0

## [1] NaN

# Is NaN NA?
is.na(NaN)

## [1] TRUE

# Is NaN NA?
is.na(0 / 0)

## [1] TRUE

However, NA and NaN are not truly identical

# Are NaN and NA identical?
identical(NA, NaN)

## [1] FALSE

# Are they equal?
NA == NaN

## [1] NA

So, what does all of this NA stuff mean for your data work? One of the most common things you will do when working with data is subsetting/filtering. For instance, let’s define a (very simple) data frame called test_df.

# Create the data frame
test_df <- data.frame(
  x = c(NA, "A", "B", NA, "A"),
  y = c(1:4, NA))
# Print test_df to the screen
test_df

##      x  y
## 1 <NA>  1
## 2    A  2
## 3    B  3
## 4 <NA>  4
## 5    A NA

Now, let’s use the dplyr function filter() to grab all observations whose value for the variable x is equal to "A".³

# Subset to x == A
dplyr::filter(test_df, x == "A")

##   x  y
## 1 A  2
## 2 A NA

Notice that we get only the values of x that are equal to "A"—meaning we do not get values of x that are equal to "B" or NA.

What if we take the opposite—those values not equal to "A"?

# Subset to x != A
dplyr::filter(test_df, x != "A")

##   x y
## 1 B 3

Notice here that we get only observations with x equal to "B"—we still do not get values of x equal to NA. Why? It’s because

NA == "A"

## [1] NA

NA == "B"

## [1] NA

So what do we do if we want values of x equal to both "A" and NA? is.na()! (Also using the logical operator for or, i.e., |).

# x equal to A or NA
dplyr::filter(test_df, x == "A" | is.na(x))

##      x  y
## 1 <NA>  1
## 2    A  2
## 3 <NA>  4
## 4    A NA

There we go.

Finally, note that when you read a .csv file, you might need to tell R which characters should be considered as NA. By default, read_csv() (in the readr package) reads "" and "NA" as NA. However, you might know that the World Bank also uses ".." for missing data. Thus, you would want the following code to read in a World Bank file:

wb_df <- read_csv(
  file = "world_bank.csv",
  na = c("", "NA", ".."))

Aside: Load your data

Our OLS function will need some data. Load the auto.dta data from Section 1 (also in this section’s zip file). (Remember: you will need the haven package to load the .dta file.) We’re not loading the data inside our function because we’ll probably want to use the function on different datasets.

# Setup ----
# Options
options(stringsAsFactors = F)
# Packages
library(pacman)
p_load(haven, dplyr)
# Define directories
dir_class <- "/Users/edwardarubin/Dropbox/Teaching/ARE212/"
dir_sec3 <- paste0(dir_class, "Section03/")

# Load the data ----
cars <- read_dta(
  file = paste0(dir_sec3, "auto.dta"))

required packages

Spoiler: Our function is going to make use of the dplyr package. So let’s tell our function to make sure the dplyr package is loaded. The function require() is the standard way to have a function make sure a package is loaded. You use it just like the library()⁹ function. Since we know that we plan to use the dplyr package, let’s require it within our function:

b_ols <- function(data, y, X) {
  # Require the 'dplyr' package
  require(dplyr)
  # Put some code here...
  return(beta_hat)
}

select_ing variables

Let’s take an inventory of which objects we have, once we are inside the function. We have data, which is a tibble with columns that represent various variables. We have y, the name of our outcome variable (e.g., weight). And we have X, a vector of the names of our independent variables (e.g. c("mpg", "weight")).¹⁰

The first step for our function is to grab the data for y and X from data. For this task, we will use a variation of the select() function introduced in Section 1: select_(). The difference between select() and select_() (besides the underscore) is that select() wants the variable names without quotes (non-standard evaluation), e.g. select(cars, mpg, weight). This notation is pretty convenient except when you are writing your own function. Generally, you will have variable/column names in a character vector, and select(cars, "mpg", "weight") does not work. Here is where select_() comes in: it wants you to use characters (standard evaluation). There is one more complexity: while select_(cars, "mpg", "weight") works, select_(cars, c("mpg", "weight")) does not. So if you have a vector of variable names, like our X, you need a slightly different way to use select_(). The solution is the .dots argument in select_(): select_(cars, .dots = c("mpg", "weight")) works!

So… we now want to select the y and X variables from data. Let’s do it.

# Select y variable data from 'data'
y_data <- select_(data, .dots = y)
# Select X variable data from 'data'
X_data <- select_(data, .dots = X)

This code should do the trick. To test it, you’ll need to define y and X (e.g., y = "price" and X = c("mpg", "weight")).

Exercise: Finish the function

The function now looks like

b_ols <- function(data, y, X) {
  # Require the 'dplyr' package
  require(dplyr)
  # Select y variable data from 'data'
  y_data <- select_(data, .dots = y)
  # Select X variable data from 'data'
  X_data <- select_(data, .dots = X)
  # Put some code here...
  return(beta_hat)
}

Fill in the # Put some code here... section of our new function with the code needed to produce OLS estimates via matrix operations. More kudos for fewer lines.

Matrices

We have a few tasks left:

1. Add an intercept (column of ones) to X_data
2. Convert the data objects to matrices
3. Calculate \(\widehat{\boldsymbol{\beta}}_{ols}\) via matrix operations

First, let’s add a column of ones to X_data. We will use mutate_().¹¹ The mutate() and mutate_() functions allow us to add new columns/variables to an existing data object. Often the new variables will be a combination of existing variables, but in our case, we just want a column of ones, so all we need to do is write mutate_(X_data, "ones" = 1).

It is customary to have the intercept column be the first column in the matrix. We can use select_() again to change the order of the columns: select_(X_data, "ones", .dots = X).

We will use the as.matrix() function to convert our tibbles to matrices.

Finally, once we have our matrices, we can use the basic matrix functions discussed in Section 2—namely %*%, t(), and solve()—to calculate \(\widehat{\boldsymbol{\beta}}_{ols} = \left(\mathbf{X}'\mathbf{X}\right)^{-1}\mathbf{X}'\mathbf{y}\).

Putting these steps together, we can finish our function:

b_ols <- function(data, y, X) {
  # Require the 'dplyr' package
  require(dplyr)

  # Select y variable data from 'data'
  y_data <- select_(data, .dots = y)
  # Convert y_data to matrices
  y_data <- as.matrix(y_data)

  # Select X variable data from 'data'
  X_data <- select_(data, .dots = X)
  # Add a column of ones to X_data
  X_data <- mutate_(X_data, "ones" = 1)
  # Move the intercept column to the front
  X_data <- select_(X_data, "ones", .dots = X)
  # Convert X_data to matrices
  X_data <- as.matrix(X_data)

  # Calculate beta hat
  beta_hat <- solve(t(X_data) %*% X_data) %*% t(X_data) %*% y_data
  # Change the name of 'ones' to 'intercept'
  rownames(beta_hat) <- c("intercept", X)
  # Return beta_hat
  return(beta_hat)
}

lapply() meets b_ols()

What if we want to regress each of the numerical variables in the cars data on mpg and weight (with the exception of rep78, because I don’t really understand what “Repair Record 1978” means)? Surprise, surprise: we can use lapply().

What should our X value be? The numeric variables excluding rep78, mpg, and weight. Let’s create a vector for it.

target_vars <- c("price", "headroom", "trunk", "length", "turn",
  "displacement", "gear_ratio", "foreign")

With respect to the FUN argument, keep in mind that lapply() plugs the X values into the first argument of the function. For b_ols(), the first argument is data, which is not what we currently want to vary. We want to vary y, which is the second argument. Rather than redefining the b_ols() function, we can augment it by wrapping another function around it.¹⁶ For example,

function(i) b_ols(data = cars, y = i, X = c("mpg", "weight"))

This line of code creates a new, unnamed function with one argument i. The argument i is then fed to our b_ols() function as its y argument. Let’s put it all together…

# The 'lapply' call
results_list <- lapply(
  X = target_vars,
  FUN = function(i) b_ols(data = cars, y = i, X = c("mpg", "weight"))
  )
# The results
results_list

## [[1]]
##                 price
## intercept  -49.512221
## mpg          1.746559
## weight    1946.068668
## 
## [[2]]
##                headroom
## intercept -0.0098904309
## mpg        0.0004668253
## weight     1.7943225731
## 
## [[3]]
##                  trunk
## intercept -0.082739270
## mpg        0.003202433
## weight     5.849262628
## 
## [[4]]
##                 length
## intercept  -0.35546594
## mpg         0.02496695
## weight    120.11619444
## 
## [[5]]
##                   turn
## intercept -0.059092537
## mpg        0.004498541
## weight    27.323996368
## 
## [[6]]
##           displacement
## intercept    0.7604918
## mpg          0.1103151
## weight    -151.9910285
## 
## [[7]]
##              gear_ratio
## intercept  0.0007521123
## mpg       -0.0004412382
## weight     4.3311476331
## 
## [[8]]
##                 foreign
## intercept -0.0194295266
## mpg       -0.0004677698
## weight     2.1235056112

These results are a bit of a mess. Let’s change the list into a more legible data structure. We will use lapply() to apply the function data.frame() to each of the results (each of the elements of results_list). Finally, we will use the bind_cols() function from dplyr to bind all of the results together (so we don’t end up with another list).¹⁷

# Cleaning up the results list
results_df <- lapply(X = results_list, FUN = data.frame) %>% bind_cols()
# We lose the row names in the process; add them back
rownames(results_df) <- c("intercept", "mpg", "weight")
# Check out results_df
results_df

##                 price      headroom        trunk       length         turn
## intercept  -49.512221 -0.0098904309 -0.082739270  -0.35546594 -0.059092537
## mpg          1.746559  0.0004668253  0.003202433   0.02496695  0.004498541
## weight    1946.068668  1.7943225731  5.849262628 120.11619444 27.323996368
##           displacement    gear_ratio       foreign
## intercept    0.7604918  0.0007521123 -0.0194295266
## mpg          0.1103151 -0.0004412382 -0.0004677698
## weight    -151.9910285  4.3311476331  2.1235056112

Exercise: Check your work

Check the results in results_df using lapply() and felm(). Hint: remember to check the class of the object returned felm(). You might want to try the coef() function on the object returned by felm().