Introduction

Here you can find lab notes and resources for Econ 415. These will be updated after our in-class lab sessions. These notes are not a substitute for attending lab but serve as an additional resource.

Much of the lab content will be drawn from R for Data Science

Useful R Resources

Getting Started with R – A collection of resources for those getting started with R

TidyR Cheatsheet – A useful cheatsheet for data cleaning and tidy data using the tidyverse functions

ggplot2 Cheatsheet – A useful cheatsheet for various ggplot geoms

Tidy Tuesday Repo – A weekly data science project in R to test your tidyverse skills!

R Setup and Workflow Basics

Before R: File Management

Effective file folder management is crucial for maintaining an organized and efficient digital workspace. Setting up organized folders will make your life significantly easier in the future!

  • Use your computer’s file management tools or RStudio’s “Files” tab.
  • Create a main folder for your project or course.
  • Go to the desired folder and create a new folder, e.g., “Lab1”
  • Save all downloaded data files and R Scripts in this folder

A Look Around R Studio

In the bottom left of R Studio, you will see the console. The console executes code. You can type code and execute it using the console but the code is not saved when you close R Studio. It is recommended that you do not use the console in your regular workflow.

To save your work, you should code in an R Script. Open a script using the button that looks like a piece of paper with a green plus sign in the top left corner of R Studio.

R scripts will open here. You can code, comment, and run the code from your script. To run the code, either click the “Run” button or by pressing CMD+Enter (Mac) or Ctrl+Enter (Windows). R scripts will be saved to the folder you are currently working in.

In the top left corner, we have the workspace/environment panes.

The workspace/enviroment tab tells you what objects are stored in R (i.e. what is loaded or stored in memory). The History tab which shows previous commands you have run.

Last, on the bottom right, we have several tabs including:

  • Files - shows the files on your computer in the directory you are working in
  • Viewer - can vew data or R objects
  • Help - shows help documentations for R functions and datasets
  • Plots - can see current and previous plots generated in your R session, save, and export them to png/pdf formats.
  • Packages - list of R packages you have installed. You can also install packages directly from this tab.

Coding Basics

R uses object-oriented programming. If you have never used this type of programming before, it can be a bit confusing at first. Essentially, R uses functions, which we apply to objects. More on this shortly, but if you aren’t sure what a function does, or how it works, you can use ? before the function to get the documentation. Ex: ?mean will bring up the help page for the mean() function. Try typing ?mean in the console and looking at the help page.

Objects

An object is an assignment between a name and a value. You assign values to names using <- or =. The first assignment symbol consists of a < next to a dash - to make it look like an arrow.

x <- 5 #assign the value of 5 to a variable called x
# notice that this x is now in your global environment
x # print x
## [1] 5
y = 10
y
## [1] 10

You can combine objects together as well which lets us do some basic math operations.

# create a new object called z that is equal to x*y
z <- x * y
#print z
z
## [1] 50

If you do not create an object, R will not save it to the global environment. If an object is not in the global environment and you try to reference it later, R will not know what you are referring to.

Math Operations

a <- 2+3
a
## [1] 5
b<-4-5
b
## [1] -1
c<-4*2
c
## [1] 8
d<-6/3
d
## [1] 2
e<-7^2
e
## [1] 49

Vectors

You can create a vector (a list) of items in R.

# create a vector of 1 through 10
vector1 <- 1:10
vector1
##  [1]  1  2  3  4  5  6  7  8  9 10

If we want specific items, we use the c() function and separate the items with a comma.

vector2 <- c(1,3,5,7,9)
vector2
## [1] 1 3 5 7 9

Mathematical operations work on vectors too!

vector2^2
## [1]  1  9 25 49 81

Classes

Objects in R have different classes. Check the class of a few objects we have already created:

class(x)
## [1] "numeric"
class(vector1)
## [1] "integer"

There are other classes too!

# create a string
my_string <- "Econ is cool!"
class(my_string)
## [1] "character"
# logical class
class(2>3)
## [1] "logical"

What happens if we have a vector of characters and numbers?

char_vector <- c(1:5, "banana", "apple")
char_vector
## [1] "1"      "2"      "3"      "4"      "5"      "banana" "apple"
#cant use mathematical operations on characters
# why?? because the entire vector is a character class!
class(char_vector)
## [1] "character"

Functions

Functions are operations that can transform your created object in a ton of different ways. We have actually already used two functions, c() and class(). Here are a few other useful ones:

#print the first few objects in vector1
head(vector1)
## [1] 1 2 3 4 5 6
#print the first 2 objects in vector1
head(vector1, 2)
## [1] 1 2
#print the last few objects in vector1
tail(vector1)
## [1]  5  6  7  8  9 10
#print last two objects in vector1
tail(vector1, 2)
## [1]  9 10
#find the mean of vector1
mean(vector1)
## [1] 5.5
#median 
median(vector1)
## [1] 5.5
#standard deviation
sd(vector1)
## [1] 3.02765
#Summary() prints summary stats
summary(vector1)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00    3.25    5.50    5.50    7.75   10.00
min(vector1)
## [1] 1
max(vector1)
## [1] 10

Code Style

Coding style is the punctuation and grammer of the coding world. Making sure that your code is formatting in a readable, standard format is helpful for yourself and others to understand the code. We will follow guidelines from the tidyverse style guide

Spaces: Put spaces on either side of mathematical operators (i.e. +, -, ==, <, …), and around the assignment operator (<-). The exception to this is the ^ symbol.

# Strive for
z <- (a + b)^2 / d

# Avoid
z<-( a + b ) ^ 2/d

Don’t put spaces inside or outside parentheses for regular function calls. Always put a space after a comma.

# Strive for
mean(x, na.rm = TRUE)
## [1] 5
# Avoid
mean (x ,na.rm=TRUE)
## [1] 5

Adding extra spaces is fine if it helps with alignment. For example:

example_data_frame <- 
  data.frame(
    variable1      = c(1:10),
    variable_name2 = c(2:11),
    var_name       = c(3:12)
  )

Naming Conventions: Object names must start with a letter and can only contain letters, numbers, _, and .. The names should be descriptive–snake case is the recommended naming convention (separating lowercase wrods with _).

really_long_variable_name <- 1

Commenting: You can comment your code with #. It is strongly recommended to leave comments in your code so that others, and future you, can keep track of your thought process.

# Good code is well-commented code!!

You can also create section comments that will be collapsable. This is incredibly helpful when you have a really long R script! Any comment line which includes at least four trailing dashes (-) will create a section.

# This is section 1 ----

# ---- This is section 2 ----

Packages

R is really useful because of its ability to use packages. Pacman is a package for “package management” - it helps us load multiple packages at onc.

# if you have not previously installed the package, include the line:
#install.packages("pacman")

# you only have to do this once. you can also install packages from the "Packages" side panel tab

We need to load the a package after installing it to use it by using library().

library(pacman)

Now we use the p_load function to load other packages we want to use. We will use the tidyverse() package throughout the course, so let’s load that one.

p_load(tidyverse)

Tidyverse Functions

Tidyverse is used for data wrangling. It allows you to manipulate data frames in a rather intuitive way. Tidyverse is a huge package so today we will be focusing on functions from the dplyr package (comes with tidyverse).

  • select(): subset columns
  • filter(): subset rows on conditions
  • arrange(): sort results
  • mutate(): create new columns by using information from other columns
  • group_by() and summarize(): create summary statisitcs on grouped data
  • count(): count discrete values

Let’s use the starwars dataset that is built into the tidyverse package.

names(starwars)
##  [1] "name"       "height"     "mass"       "hair_color" "skin_color"
##  [6] "eye_color"  "birth_year" "sex"        "gender"     "homeworld" 
## [11] "species"    "films"      "vehicles"   "starships"
head(starwars)
## # A tibble: 6 × 14
##   name      height  mass hair_color skin_color eye_color birth_year sex   gender
##   <chr>      <int> <dbl> <chr>      <chr>      <chr>          <dbl> <chr> <chr> 
## 1 Luke Sky…    172    77 blond      fair       blue            19   male  mascu…
## 2 C-3PO        167    75 <NA>       gold       yellow         112   none  mascu…
## 3 R2-D2         96    32 <NA>       white, bl… red             33   none  mascu…
## 4 Darth Va…    202   136 none       white      yellow          41.9 male  mascu…
## 5 Leia Org…    150    49 brown      light      brown           19   fema… femin…
## 6 Owen Lars    178   120 brown, gr… light      blue            52   male  mascu…
## # ℹ 5 more variables: homeworld <chr>, species <chr>, films <list>,
## #   vehicles <list>, starships <list>

Select and Filter

Let’s select only the name, gender, and homeworld variables

select(starwars, c(name, gender, homeworld))
## # A tibble: 87 × 3
##    name               gender    homeworld
##    <chr>              <chr>     <chr>    
##  1 Luke Skywalker     masculine Tatooine 
##  2 C-3PO              masculine Tatooine 
##  3 R2-D2              masculine Naboo    
##  4 Darth Vader        masculine Tatooine 
##  5 Leia Organa        feminine  Alderaan 
##  6 Owen Lars          masculine Tatooine 
##  7 Beru Whitesun Lars feminine  Tatooine 
##  8 R5-D4              masculine Tatooine 
##  9 Biggs Darklighter  masculine Tatooine 
## 10 Obi-Wan Kenobi     masculine Stewjon  
## # ℹ 77 more rows

Notice that this didn’t save anything in our global environment! If you want to save this new dataframe, you have to give it a name!

To select all columns except a certain one, use a minus sign

select(starwars, c(-homeworld, -gender))
## # A tibble: 87 × 12
##    name    height  mass hair_color skin_color eye_color birth_year sex   species
##    <chr>    <int> <dbl> <chr>      <chr>      <chr>          <dbl> <chr> <chr>  
##  1 Luke S…    172    77 blond      fair       blue            19   male  Human  
##  2 C-3PO      167    75 <NA>       gold       yellow         112   none  Droid  
##  3 R2-D2       96    32 <NA>       white, bl… red             33   none  Droid  
##  4 Darth …    202   136 none       white      yellow          41.9 male  Human  
##  5 Leia O…    150    49 brown      light      brown           19   fema… Human  
##  6 Owen L…    178   120 brown, gr… light      blue            52   male  Human  
##  7 Beru W…    165    75 brown      light      blue            47   fema… Human  
##  8 R5-D4       97    32 <NA>       white, red red             NA   none  Droid  
##  9 Biggs …    183    84 black      light      brown           24   male  Human  
## 10 Obi-Wa…    182    77 auburn, w… fair       blue-gray       57   male  Human  
## # ℹ 77 more rows
## # ℹ 3 more variables: films <list>, vehicles <list>, starships <list>

Filter the data frame to include only droids

filter(starwars, species == "Droid")
## # A tibble: 6 × 14
##   name   height  mass hair_color skin_color  eye_color birth_year sex   gender  
##   <chr>   <int> <dbl> <chr>      <chr>       <chr>          <dbl> <chr> <chr>   
## 1 C-3PO     167    75 <NA>       gold        yellow           112 none  masculi…
## 2 R2-D2      96    32 <NA>       white, blue red               33 none  masculi…
## 3 R5-D4      97    32 <NA>       white, red  red               NA none  masculi…
## 4 IG-88     200   140 none       metal       red               15 none  masculi…
## 5 R4-P17     96    NA none       silver, red red, blue         NA none  feminine
## 6 BB8        NA    NA none       none        black             NA none  masculi…
## # ℹ 5 more variables: homeworld <chr>, species <chr>, films <list>,
## #   vehicles <list>, starships <list>

Filter the data frame to include droids OR humans

filter(starwars, species == "Droid" | species == "Human")
## # A tibble: 41 × 14
##    name     height  mass hair_color skin_color eye_color birth_year sex   gender
##    <chr>     <int> <dbl> <chr>      <chr>      <chr>          <dbl> <chr> <chr> 
##  1 Luke Sk…    172    77 blond      fair       blue            19   male  mascu…
##  2 C-3PO       167    75 <NA>       gold       yellow         112   none  mascu…
##  3 R2-D2        96    32 <NA>       white, bl… red             33   none  mascu…
##  4 Darth V…    202   136 none       white      yellow          41.9 male  mascu…
##  5 Leia Or…    150    49 brown      light      brown           19   fema… femin…
##  6 Owen La…    178   120 brown, gr… light      blue            52   male  mascu…
##  7 Beru Wh…    165    75 brown      light      blue            47   fema… femin…
##  8 R5-D4        97    32 <NA>       white, red red             NA   none  mascu…
##  9 Biggs D…    183    84 black      light      brown           24   male  mascu…
## 10 Obi-Wan…    182    77 auburn, w… fair       blue-gray       57   male  mascu…
## # ℹ 31 more rows
## # ℹ 5 more variables: homeworld <chr>, species <chr>, films <list>,
## #   vehicles <list>, starships <list>

Filter the data frame to include characters taler than 100 cm and a mass over 100

filter(starwars, height > 100 & mass > 100)
## # A tibble: 10 × 14
##    name     height  mass hair_color skin_color eye_color birth_year sex   gender
##    <chr>     <int> <dbl> <chr>      <chr>      <chr>          <dbl> <chr> <chr> 
##  1 Darth V…    202   136 none       white      yellow          41.9 male  mascu…
##  2 Owen La…    178   120 brown, gr… light      blue            52   male  mascu…
##  3 Chewbac…    228   112 brown      unknown    blue           200   male  mascu…
##  4 Jabba D…    175  1358 <NA>       green-tan… orange         600   herm… mascu…
##  5 Jek Ton…    180   110 brown      fair       blue            NA   <NA>  <NA>  
##  6 IG-88       200   140 none       metal      red             15   none  mascu…
##  7 Bossk       190   113 none       green      red             53   male  mascu…
##  8 Dexter …    198   102 none       brown      yellow          NA   male  mascu…
##  9 Grievous    216   159 none       brown, wh… green, y…       NA   male  mascu…
## 10 Tarfful     234   136 brown      brown      blue            NA   male  mascu…
## # ℹ 5 more variables: homeworld <chr>, species <chr>, films <list>,
## #   vehicles <list>, starships <list>

Piping

What if we want to do those things all in one step??? The tidyverse allows us to chain functions together with %>%. This is called a pipe and the pipe connects the LHS to the RHS (like reading a book). Let’s make a new dataframe where we select the name, height, and mass. Filter out those who are shorter than 100 cm.

new_df <- starwars %>% select(name, height, mass) %>% filter(height >= 100)
new_df
## # A tibble: 74 × 3
##    name               height  mass
##    <chr>               <int> <dbl>
##  1 Luke Skywalker        172    77
##  2 C-3PO                 167    75
##  3 Darth Vader           202   136
##  4 Leia Organa           150    49
##  5 Owen Lars             178   120
##  6 Beru Whitesun Lars    165    75
##  7 Biggs Darklighter     183    84
##  8 Obi-Wan Kenobi        182    77
##  9 Anakin Skywalker      188    84
## 10 Wilhuff Tarkin        180    NA
## # ℹ 64 more rows

Self check: make a new data frame where you select all columns except gender and has characters that have green skin color

example_df <- starwars %>% select(-gender) %>% filter(skin_color == "green")
example_df
## # A tibble: 6 × 13
##   name   height  mass hair_color skin_color eye_color birth_year sex   homeworld
##   <chr>   <int> <dbl> <chr>      <chr>      <chr>          <dbl> <chr> <chr>    
## 1 Greedo    173    74 <NA>       green      black             44 male  Rodia    
## 2 Yoda       66    17 white      green      brown            896 male  <NA>     
## 3 Bossk     190   113 none       green      red               53 male  Trandosha
## 4 Rugor…    206    NA none       green      orange            NA male  Naboo    
## 5 Kit F…    196    87 none       green      black             NA male  Glee Ans…
## 6 Poggl…    183    80 none       green      yellow            NA male  Geonosis 
## # ℹ 4 more variables: species <chr>, films <list>, vehicles <list>,
## #   starships <list>

Arrange

Let’s arrange all of the characters by their height

starwars %>% arrange(height)
## # A tibble: 87 × 14
##    name     height  mass hair_color skin_color eye_color birth_year sex   gender
##    <chr>     <int> <dbl> <chr>      <chr>      <chr>          <dbl> <chr> <chr> 
##  1 Yoda         66    17 white      green      brown            896 male  mascu…
##  2 Ratts T…     79    15 none       grey, blue unknown           NA male  mascu…
##  3 Wicket …     88    20 brown      brown      brown              8 male  mascu…
##  4 Dud Bolt     94    45 none       blue, grey yellow            NA male  mascu…
##  5 R2-D2        96    32 <NA>       white, bl… red               33 none  mascu…
##  6 R4-P17       96    NA none       silver, r… red, blue         NA none  femin…
##  7 R5-D4        97    32 <NA>       white, red red               NA none  mascu…
##  8 Sebulba     112    40 none       grey, red  orange            NA male  mascu…
##  9 Gasgano     122    NA none       white, bl… black             NA male  mascu…
## 10 Watto       137    NA black      blue, grey yellow            NA male  mascu…
## # ℹ 77 more rows
## # ℹ 5 more variables: homeworld <chr>, species <chr>, films <list>,
## #   vehicles <list>, starships <list>

Notice this does lowest to highest, we can do the other way too

starwars %>% arrange(desc(height))
## # A tibble: 87 × 14
##    name     height  mass hair_color skin_color eye_color birth_year sex   gender
##    <chr>     <int> <dbl> <chr>      <chr>      <chr>          <dbl> <chr> <chr> 
##  1 Yarael …    264    NA none       white      yellow          NA   male  mascu…
##  2 Tarfful     234   136 brown      brown      blue            NA   male  mascu…
##  3 Lama Su     229    88 none       grey       black           NA   male  mascu…
##  4 Chewbac…    228   112 brown      unknown    blue           200   male  mascu…
##  5 Roos Ta…    224    82 none       grey       orange          NA   male  mascu…
##  6 Grievous    216   159 none       brown, wh… green, y…       NA   male  mascu…
##  7 Taun We     213    NA none       grey       black           NA   fema… femin…
##  8 Rugor N…    206    NA none       green      orange          NA   male  mascu…
##  9 Tion Me…    206    80 none       grey       black           NA   male  mascu…
## 10 Darth V…    202   136 none       white      yellow          41.9 male  mascu…
## # ℹ 77 more rows
## # ℹ 5 more variables: homeworld <chr>, species <chr>, films <list>,
## #   vehicles <list>, starships <list>

Self check: Arrange the characters names in alphabetical order

starwars %>% arrange(name)
## # A tibble: 87 × 14
##    name     height  mass hair_color skin_color eye_color birth_year sex   gender
##    <chr>     <int> <dbl> <chr>      <chr>      <chr>          <dbl> <chr> <chr> 
##  1 Ackbar      180    83 none       brown mot… orange          41   male  mascu…
##  2 Adi Gal…    184    50 none       dark       blue            NA   fema… femin…
##  3 Anakin …    188    84 blond      fair       blue            41.9 male  mascu…
##  4 Arvel C…     NA    NA brown      fair       brown           NA   male  mascu…
##  5 Ayla Se…    178    55 none       blue       hazel           48   fema… femin…
##  6 BB8          NA    NA none       none       black           NA   none  mascu…
##  7 Bail Pr…    191    NA black      tan        brown           67   male  mascu…
##  8 Barriss…    166    50 black      yellow     blue            40   fema… femin…
##  9 Ben Qua…    163    65 none       grey, gre… orange          NA   male  mascu…
## 10 Beru Wh…    165    75 brown      light      blue            47   fema… femin…
## # ℹ 77 more rows
## # ℹ 5 more variables: homeworld <chr>, species <chr>, films <list>,
## #   vehicles <list>, starships <list>

Mutate

We can create new variables with mutate(). Let’s create a new variable that measures height in inches instead of centimeters (2.54cm per inch)

starwars %>% mutate(height_inches = height/2.54)
## # A tibble: 87 × 15
##    name     height  mass hair_color skin_color eye_color birth_year sex   gender
##    <chr>     <int> <dbl> <chr>      <chr>      <chr>          <dbl> <chr> <chr> 
##  1 Luke Sk…    172    77 blond      fair       blue            19   male  mascu…
##  2 C-3PO       167    75 <NA>       gold       yellow         112   none  mascu…
##  3 R2-D2        96    32 <NA>       white, bl… red             33   none  mascu…
##  4 Darth V…    202   136 none       white      yellow          41.9 male  mascu…
##  5 Leia Or…    150    49 brown      light      brown           19   fema… femin…
##  6 Owen La…    178   120 brown, gr… light      blue            52   male  mascu…
##  7 Beru Wh…    165    75 brown      light      blue            47   fema… femin…
##  8 R5-D4        97    32 <NA>       white, red red             NA   none  mascu…
##  9 Biggs D…    183    84 black      light      brown           24   male  mascu…
## 10 Obi-Wan…    182    77 auburn, w… fair       blue-gray       57   male  mascu…
## # ℹ 77 more rows
## # ℹ 6 more variables: homeworld <chr>, species <chr>, films <list>,
## #   vehicles <list>, starships <list>, height_inches <dbl>

Notice that this new variable is not in the original data frame. Why? Because we didn’t assign it to our object.

starwars <- starwars %>% mutate(height_inches = height/2.54)

Self check: Create a new variable that is the sum of person’s mass and height

starwars %>% mutate(total = height + mass)
## # A tibble: 87 × 16
##    name     height  mass hair_color skin_color eye_color birth_year sex   gender
##    <chr>     <int> <dbl> <chr>      <chr>      <chr>          <dbl> <chr> <chr> 
##  1 Luke Sk…    172    77 blond      fair       blue            19   male  mascu…
##  2 C-3PO       167    75 <NA>       gold       yellow         112   none  mascu…
##  3 R2-D2        96    32 <NA>       white, bl… red             33   none  mascu…
##  4 Darth V…    202   136 none       white      yellow          41.9 male  mascu…
##  5 Leia Or…    150    49 brown      light      brown           19   fema… femin…
##  6 Owen La…    178   120 brown, gr… light      blue            52   male  mascu…
##  7 Beru Wh…    165    75 brown      light      blue            47   fema… femin…
##  8 R5-D4        97    32 <NA>       white, red red             NA   none  mascu…
##  9 Biggs D…    183    84 black      light      brown           24   male  mascu…
## 10 Obi-Wan…    182    77 auburn, w… fair       blue-gray       57   male  mascu…
## # ℹ 77 more rows
## # ℹ 7 more variables: homeworld <chr>, species <chr>, films <list>,
## #   vehicles <list>, starships <list>, height_inches <dbl>, total <dbl>

Group_by and Summarize

This will group data together and can make summary statistics. Let’s find the average height for each species.

starwars %>% group_by(species) %>% summarize(avg_height = mean(height))
## # A tibble: 38 × 2
##    species   avg_height
##    <chr>          <dbl>
##  1 Aleena           79 
##  2 Besalisk        198 
##  3 Cerean          198 
##  4 Chagrian        196 
##  5 Clawdite        168 
##  6 Droid            NA 
##  7 Dug             112 
##  8 Ewok             88 
##  9 Geonosian       183 
## 10 Gungan          209.
## # ℹ 28 more rows

Notice we have NA’s! We can get rid of those.

# Using na.omit()
starwars %>% na.omit() %>% group_by(species) %>% summarize(avg_height = mean(height))
## # A tibble: 11 × 2
##    species      avg_height
##    <chr>             <dbl>
##  1 Cerean             198 
##  2 Ewok                88 
##  3 Gungan             196 
##  4 Human              179.
##  5 Kel Dor            188 
##  6 Mirialan           168 
##  7 Mon Calamari       180 
##  8 Trandoshan         190 
##  9 Twi'lek            178 
## 10 Wookiee            228 
## 11 Zabrak             175
# Using na.rm = T inside the mean function
starwars %>%  group_by(species) %>% summarize(avg_height = mean(height, na.rm = T))
## # A tibble: 38 × 2
##    species   avg_height
##    <chr>          <dbl>
##  1 Aleena           79 
##  2 Besalisk        198 
##  3 Cerean          198 
##  4 Chagrian        196 
##  5 Clawdite        168 
##  6 Droid           131.
##  7 Dug             112 
##  8 Ewok             88 
##  9 Geonosian       183 
## 10 Gungan          209.
## # ℹ 28 more rows

Count

Count the number of each species.

starwars %>% count(species)
## # A tibble: 38 × 2
##    species       n
##    <chr>     <int>
##  1 Aleena        1
##  2 Besalisk      1
##  3 Cerean        1
##  4 Chagrian      1
##  5 Clawdite      1
##  6 Droid         6
##  7 Dug           1
##  8 Ewok          1
##  9 Geonosian     1
## 10 Gungan        3
## # ℹ 28 more rows

Graphing

Let’s use the penguins dataset that is in the palmerpenguins package. Also load the tidyverse package.

library(pacman)
p_load(tidyverse, palmerpenguins)

The dataset contains 344 observations of penguins at the Palmer Station in Antartica.

head(penguins)
## # A tibble: 6 × 8
##   species island    bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
##   <fct>   <fct>              <dbl>         <dbl>             <int>       <int>
## 1 Adelie  Torgersen           39.1          18.7               181        3750
## 2 Adelie  Torgersen           39.5          17.4               186        3800
## 3 Adelie  Torgersen           40.3          18                 195        3250
## 4 Adelie  Torgersen           NA            NA                  NA          NA
## 5 Adelie  Torgersen           36.7          19.3               193        3450
## 6 Adelie  Torgersen           39.3          20.6               190        3650
## # ℹ 2 more variables: sex <fct>, year <int>

Let’s use this dataset to create plots, following the code in “R for Data Science”.

Creating a ggplot

The basic setup of making a ggplot requires three things: the data, the aesthetic mapping, and a geom. The aesthetic mappings describe how variables in the data are mapped to visual properties (aesthetics) of geoms, like which variables are on the axes, the variable to color or fill by, etc. The geoms tell R how to draw the data like points, lines, columns, etc.

In general, we can make a ggplot by typing the following: \[\text{ggplot(data = <DATA>, mapping = aes(<MAPPING>)) + <geom_function>()}\]

The way ggplot works is by adding layers. We can add a new layer with the + sign. Let’s build a ggplot step by step. First, start with ggplot() and tell R what data we are using.

ggplot(data = penguins)

Why did this make a blank graph? Well, we haven’t given R the aesthetic mapping yet so it doesn’t know what to put on top of the base layer. Let’s add the x and y variables, flipper length and body mass.

ggplot(data = penguins, mapping = aes(x = flipper_length_mm, y = body_mass_g)) # note you can put the mapping here or in the geom

Now we have a graph with axes and gridlines but no information on the graph. To get data on the graph, we need to tell R how we want to draw the data with a geom. To make a scatterplot, we use geom_point().

ggplot(data = penguins, mapping = aes(x = flipper_length_mm, y = body_mass_g)) + geom_point()

Adding Aesthetics and Layers

Suppose we want to see how the relationship between flipper length and body mass differs by species of penguin. We can represent each species on our scatterplot with different colored points.

ggplot(data = penguins,
       mapping = aes(x = flipper_length_mm, y = body_mass_g, color = species)) +
  geom_point()

Let’s add another layer: a curve displaying the relationship using geom_smooth(). Using method = "lm" gives the line of best fit.

ggplot(data = penguins,
       mapping = aes(x = flipper_length_mm, y = body_mass_g, color = species)) +
  geom_point() +
  geom_smooth(method = "lm")

Note that this creates three separate lines, one for each species. When the aesthetic mapping is defined in the `ggplot function, the mappings are passed down to all of the subsequent layers. You can also define aesthetic mappings at the local level. For example, what if we wanted the points to be colored by species but only one line of best fit?

ggplot(data = penguins,
       mapping = aes(x = flipper_length_mm, y = body_mass_g)) +
  geom_point(mapping = aes(color = species)) +
  geom_smooth(method = "lm")

We can also change the shapes of the points.

ggplot(data = penguins,
       mapping = aes(x = flipper_length_mm, y = body_mass_g)) +
  geom_point(mapping = aes(color = species, shape = species)) +
  geom_smooth(method = "lm")

Now let’s make it pretty! We can add titles, axis labels, themes and more! You can also get fancy and change the color of your points and lines.

ggplot(data = penguins,
       mapping = aes(x = flipper_length_mm, y = body_mass_g)) +
  geom_point(mapping = aes(color = species, shape = species)) +
  geom_smooth(method = "lm") +
  labs(title = "Body Mass and Flipper Length of Penguins",
       x = "Flipper Length (mm)",
       y = "Body Mass (g)") +
  theme_minimal()

Visualizing Distributions

For categorical variables, we can use a bar chart. The height of the bars is how many observations occurred for each x value.

ggplot(penguins, aes(x = species)) + geom_bar()

For numerical variables, we can create histograms and density plots. For histograms, changing the binwidth sets the intervals for the x variable.

ggplot(penguins, aes(x = body_mass_g)) + geom_histogram(binwidth = 20)

ggplot(penguins, aes(x = body_mass_g)) + geom_histogram(binwidth = 200)

ggplot(penguins, aes(x = body_mass_g)) + geom_histogram(binwidth = 2000)

Density plots are smooth versions of histograms and are useful for continuous data.

ggplot(penguins, aes(x = body_mass_g)) + geom_density()

Test Your Knowledge

Try to reproduce this graph:

Importing Data

We need to get some data before we can start working with it! Download the “lab1.csv” file from Canvas. It should then be in your Downloads folder. Now, we need to make sure that we are working in the same place as our data is, so that R knows where to get the csv file from. Move the csv file to your current directory (the folder your are currently working out of.)

If you are unsure which directory you are in, the function getwd() will give you the file path of your current location.

getwd()
## [1] "/Users/jenniputz/Dropbox/Econ415/Winter 2025/R"

Now, if this is not the same location as your desired folder, you need to tell R to set a new working directory. You can use setwd() to tell R where to go.

For example, if I want to go to my Downloads folder on my Mac, I can write setwd("/Users/jenniputz/Downloads"). If you are on windows, your file path would start with a C:\ and use \.

Once you are in your correct directory, we can read in our data using read_csv(). Remember that everything in R is an object, so you must give your data frame a name.

cps <- read_csv("lab1.csv")
## Rows: 8891 Columns: 5
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): educ
## dbl (4): employed, black, female, exper
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Non-Experimental Data

Today you will investigate racial disparities in the labor market using data from the Current Population Survey (CPS), a large survey administered by the US Census Bureau and the Bureau of Labor Statistics. The federal government uses the CPS to estimate the unemployment rate. Economists use the CPS to study a variety of topics in labor economics, including the effect of binding minimum wages, the gender pay gap, and returns to schooling. You will use a CPS sample of workers from Boston and Chicago to study employment patterns by race.

Import data

If you have not already loaded the tidyverse package, please do so.

The data file is lab1.csv and is available on Canvas. Import the data using read_csv if you have not already done so.

By looking at the dataset, you can see that most of the variables are binary: they take values of either 1 or 0.

head(cps)
## # A tibble: 6 × 5
##   employed black female educ              exper
##      <dbl> <dbl>  <dbl> <chr>             <dbl>
## 1        1     0      1 HS Graduate          20
## 2        0     0      0 HS Graduate          20
## 3        1     1      1 HS Dropout           12
## 4        1     1      0 HS Dropout           17
## 5        0     1      1 HS Graduate          21
## 6        1     0      1 College or Higher    13

For example, individuals with employed == 1 have a job while those with employed == 0 do not.

Employment Rates

What percentage of individuals in the sample are employed?

The mean of a binary variable gives the fraction of observations with values equal to 1.

mean(cps$employed)
## [1] NA

Something went wrong. If you use the summary function, you’ll see that there are missing values (NAs) of employed.

summary(cps$employed)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##  0.0000  1.0000  1.0000  0.7811  1.0000  1.0000      18

When there are missing values, some functions, like mean, will return a missing value as output. To circumvent this, you can specify na.rm = TRUE in the mean function:

mean(cps$employed, na.rm = TRUE)
## [1] 0.7811338

The employment rate is 78 percent.

What are the employment rates by race and gender?

cps %>%
  group_by(black, female) %>%
  summarize(employed = mean(employed, na.rm = TRUE))
## `summarise()` has grouped output by 'black'. You can override using the
## `.groups` argument.
## # A tibble: 4 × 3
## # Groups:   black [2]
##   black female employed
##   <dbl>  <dbl>    <dbl>
## 1     0      0    0.868
## 2     0      1    0.723
## 3     1      0    0.718
## 4     1      1    0.693

You can see that the employment rate is

  • 87 percent for white males (black == 0 and female == 0)
  • 72 percent for white females (black == 0 and female == 1)
  • 72 percent for black males (black == 1 and female == 0)
  • 69 percent for black females (black == 1 and female == 1).

Racial Disparities

What is the average difference in employment status between black individuals and white individuals?

To find the difference,

  1. Use the filter function to restrict the sample to one group (black or white)
  2. Use mean to calculate the group mean
  3. Repeat for the other group
  4. Take the difference in means.
black_emp_df <- cps %>% filter(black == 1)
black_emp <- mean(black_emp_df$employed, na.rm = T)
black_emp
## [1] 0.7035928
white_emp_df <- cps %>% filter(black == 0)
white_emp <- mean(white_emp_df$employed, na.rm = T)
white_emp
## [1] 0.7948786
black_emp - white_emp
## [1] -0.09128578

The employment rate is 9 percentage points lower for black individuals than for white individuals.

Does this mean that there is racial disparity?

Not yet. We still don’t know if the difference is statistically significant. You can find out by conducting a \(t\)-test of the null hypothesis that the true difference-in-means is zero against the alternative hypothesis that the difference is nonzero.

To conduct the test, you need to calculate the \(t\)-statistic for the difference-in-means, which is given by

\[t = \dfrac{\overline{\text{Employed}}_\text{Black} - \overline{\text{Employed}}_\text{White}}{\sqrt{\frac{S^2_\text{Black}}{N_\text{Black}} + \frac{S^2_\text{White}}{N_\text{White}}}}\] We can conduct the \(t\)-test using the t.test() function. Here, you input the variable names for x and y, and tell R if you want a two-sided or one-sided test.

t.test(black_emp_df$employed, white_emp_df$employed, alternative = "two.sided", var.equal = FALSE)
## 
##  Welch Two Sample t-test
## 
## data:  black_emp_df$employed and white_emp_df$employed
## t = -6.845, df = 1724.5, p-value = 1.059e-11
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.11744252 -0.06512905
## sample estimates:
## mean of x mean of y 
## 0.7035928 0.7948786

To conclude your test, compare your \(t\)-stat to 2 (an approximation for the critical value of \(t\) in 5 percent test–we will do more rigorous hypothesis testing later in the course). If \(|t| > 2\), then you can reject the null hypothesis. Your \(t\)-stat of -6.86 is certainly more extreme than 2, so you can reject the null hypothesis. This means that the difference in employment rates is statistically significant. There is a racial disparity in employment.

Does the disparity in employment rates provide causal evidence of racial discrimination in hiring? Does the comparison of employment rates by race hold all else constant? What else could explain the gap?

Simple OLS Regression

library(pacman)
p_load(tidyverse, broom, stargazer, AER)
# Note: broom package contains the tidy() function

For this lab, we will use data from the AER package on California schools. To get this data and use it we can:

# Note: if you install and load the AER package, the dataset is included in the package. Then run the line below:
data("CASchools")


# look at a  snapshot
head(CASchools, 10)
##    district                          school      county grades students
## 1     75119              Sunol Glen Unified     Alameda  KK-08      195
## 2     61499            Manzanita Elementary       Butte  KK-08      240
## 3     61549     Thermalito Union Elementary       Butte  KK-08     1550
## 4     61457 Golden Feather Union Elementary       Butte  KK-08      243
## 5     61523        Palermo Union Elementary       Butte  KK-08     1335
## 6     62042         Burrel Union Elementary      Fresno  KK-08      137
## 7     68536           Holt Union Elementary San Joaquin  KK-08      195
## 8     63834             Vineland Elementary        Kern  KK-08      888
## 9     62331        Orange Center Elementary      Fresno  KK-08      379
## 10    67306     Del Paso Heights Elementary  Sacramento  KK-06     2247
##    teachers calworks    lunch computer expenditure    income   english  read
## 1     10.90   0.5102   2.0408       67    6384.911 22.690001  0.000000 691.6
## 2     11.15  15.4167  47.9167      101    5099.381  9.824000  4.583333 660.5
## 3     82.90  55.0323  76.3226      169    5501.955  8.978000 30.000002 636.3
## 4     14.00  36.4754  77.0492       85    7101.831  8.978000  0.000000 651.9
## 5     71.50  33.1086  78.4270      171    5235.988  9.080333 13.857677 641.8
## 6      6.40  12.3188  86.9565       25    5580.147 10.415000 12.408759 605.7
## 7     10.00  12.9032  94.6237       28    5253.331  6.577000 68.717949 604.5
## 8     42.50  18.8063 100.0000       66    4565.746  8.174000 46.959461 605.5
## 9     19.00  32.1900  93.1398       35    5355.548  7.385000 30.079157 608.9
## 10   108.00  78.9942  87.3164        0    5036.211 11.613333 40.275921 611.9
##     math
## 1  690.0
## 2  661.9
## 3  650.9
## 4  643.5
## 5  639.9
## 6  605.4
## 7  609.0
## 8  612.5
## 9  616.1
## 10 613.4
names(CASchools)
##  [1] "district"    "school"      "county"      "grades"      "students"   
##  [6] "teachers"    "calworks"    "lunch"       "computer"    "expenditure"
## [11] "income"      "english"     "read"        "math"

Running a Regression

To do a regression in R, we use lm(). The basic steup: name <- lm(y ~ x, data = name_of_df).

Let’s regress reading scores on student expenditure.

lm(read ~ expenditure, data = CASchools)
## 
## Call:
## lm(formula = read ~ expenditure, data = CASchools)
## 
## Coefficients:
## (Intercept)  expenditure  
##   6.182e+02    6.912e-03

The output from lm() gives us an intercept coefficient, \(\hat{\beta}_0\), and a slope coefficient, \(\hat{\beta}_1\).

Let’s run another regression. Regress math scores on student expenditure.

lm(math ~ expenditure, data = CASchools)
## 
## Call:
## lm(formula = math ~ expenditure, data = CASchools)
## 
## Coefficients:
## (Intercept)  expenditure  
##   6.290e+02    4.585e-03

On your own, try to code a regression for \(Math_i = \beta_0 + \beta_1 Income_i + u_i\).

lm(math ~ income,  data = CASchools)
## 
## Call:
## lm(formula = math ~ income, data = CASchools)
## 
## Coefficients:
## (Intercept)       income  
##     625.539        1.815

Making tables

Now that we know how to run a regression, let’s talk about how to look at the output. The output above wasn’t super informative…

Using summary()

The first option is to use the summary function in R. There are two ways to do this:

Save your regression as an object in R so that we can then use that object later.

reg1 <- lm(math ~ income, data = CASchools)
summary(reg1)
## 
## Call:
## lm(formula = math ~ income, data = CASchools)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -39.045  -8.997   0.308   8.416  34.246 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 625.53948    1.53627  407.18   <2e-16 ***
## income        1.81523    0.09073   20.01   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.42 on 418 degrees of freedom
## Multiple R-squared:  0.4892, Adjusted R-squared:  0.4879 
## F-statistic: 400.3 on 1 and 418 DF,  p-value: < 2.2e-16

We have the intercept coefficient, \(\hat{\beta}_0 = 625\), and the slope coefficient \(\hat{\beta}_1 = 1.8\). summary() also gives us other information that we didn’t have before like the standard error, the t-score, the p-value, and the \(R^2\). The stars on the p-value tell us if the coefficient is statistically significant (more on this after the midterm).

Using tidy()

Another way to make nice regression output is to use the tidy() function from the broom package. To use this, you must have loaded the broom package in p_load(). The process is similar:

# since we have already created reg1 as an object, we can just call it without having to redo the regression
tidy(reg1)
## # A tibble: 2 × 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)   626.      1.54       407.  0       
## 2 income          1.82    0.0907      20.0 5.99e-63

tidy() puts the information from summary() into a much nicer looking table.

Stargazer

By far, the most powerful tool for making amazing tables in R is the stargazer package. I would encourage you to try this out if you are feeling up for it! There is a very helpful cheatsheet linked here.

First, let’s try to make a simple table with stargazer() using our reg1 object.

stargazer(reg1)
## 
## % Table created by stargazer v.5.2.3 by Marek Hlavac, Social Policy Institute. E-mail: marek.hlavac at gmail.com
## % Date and time: Thu, Feb 12, 2026 - 13:21:07
## \begin{table}[!htbp] \centering 
##   \caption{} 
##   \label{} 
## \begin{tabular}{@{\extracolsep{5pt}}lc} 
## \\[-1.8ex]\hline 
## \hline \\[-1.8ex] 
##  & \multicolumn{1}{c}{\textit{Dependent variable:}} \\ 
## \cline{2-2} 
## \\[-1.8ex] & math \\ 
## \hline \\[-1.8ex] 
##  income & 1.815$^{***}$ \\ 
##   & (0.091) \\ 
##   & \\ 
##  Constant & 625.539$^{***}$ \\ 
##   & (1.536) \\ 
##   & \\ 
## \hline \\[-1.8ex] 
## Observations & 420 \\ 
## R$^{2}$ & 0.489 \\ 
## Adjusted R$^{2}$ & 0.488 \\ 
## Residual Std. Error & 13.420 (df = 418) \\ 
## F Statistic & 400.257$^{***}$ (df = 1; 418) \\ 
## \hline 
## \hline \\[-1.8ex] 
## \textit{Note:}  & \multicolumn{1}{r}{$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01} \\ 
## \end{tabular} 
## \end{table}

Huh, weird output right? Stargazer is defaulting to TeX output. We can change the type of output we want.

1: As text (use this for your problem set)

stargazer(reg1, type = "text")
## 
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                                math            
## -----------------------------------------------
## income                       1.815***          
##                               (0.091)          
##                                                
## Constant                    625.539***         
##                               (1.536)          
##                                                
## -----------------------------------------------
## Observations                    420            
## R2                             0.489           
## Adjusted R2                    0.488           
## Residual Std. Error      13.420 (df = 418)     
## F Statistic          400.257*** (df = 1; 418)  
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

Stargazer can also do html output and MS word output (although word output is not recommended).

Stargazer gave us some stats that we don’t really care about… let’s get rid of them…

stargazer(reg1, keep.stat =c("rsq", "n"), type = "text")
## 
## ========================================
##                  Dependent variable:    
##              ---------------------------
##                         math            
## ----------------------------------------
## income                1.815***          
##                        (0.091)          
##                                         
## Constant             625.539***         
##                        (1.536)          
##                                         
## ----------------------------------------
## Observations             420            
## R2                      0.489           
## ========================================
## Note:        *p<0.1; **p<0.05; ***p<0.01

Looks great!

Unlike the tables we made using summary() and broom(), stargazer can combine multiple regressions into one table!

# do another regression and save it as an object
reg2 <- lm(read ~ income, data = CASchools)
stargazer(reg1, reg2, keep.stat =c("rsq", "n"), type = "text")
## 
## =========================================
##                  Dependent variable:     
##              ----------------------------
##                   math          read     
##                   (1)            (2)     
## -----------------------------------------
## income          1.815***      1.942***   
##                 (0.091)        (0.097)   
##                                          
## Constant       625.539***    625.228***  
##                 (1.536)        (1.651)   
##                                          
## -----------------------------------------
## Observations      420            420     
## R2               0.489          0.487    
## =========================================
## Note:         *p<0.1; **p<0.05; ***p<0.01

The variables in the regressions don’t even have to be the same!

# do another regression and save it as an object
reg3 <- lm(read ~ expenditure, data = CASchools)
stargazer(reg1, reg2, reg3, keep.stat =c("rsq", "n"), type = "text")
## 
## =============================================
##                    Dependent variable:       
##              --------------------------------
##                 math            read         
##                 (1)        (2)        (3)    
## ---------------------------------------------
## income        1.815***   1.942***            
##               (0.091)    (0.097)             
##                                              
## expenditure                         0.007*** 
##                                     (0.002)  
##                                              
## Constant     625.539*** 625.228*** 618.249***
##               (1.536)    (1.651)    (8.101)  
##                                              
## ---------------------------------------------
## Observations    420        420        420    
## R2             0.489      0.487      0.047   
## =============================================
## Note:             *p<0.1; **p<0.05; ***p<0.01

We can get really fancy with stargazer(). For more info, see the cheatsheet on Canvas.

Visualizing your regression

Last, we can use ggplot and fit a regression line through our data. Let’s start by making a scatterplot with math scores on the y axis and income on the x axis.

ggplot(data = CASchools, aes(x = income, y = math)) + geom_point()

To add a regression line, we use the stat_smooth() function:

# method = "lm" tells R that we are using OLS. se = FALSE removes the standard error bars.
ggplot(data = CASchools, aes(x = income, y = math)) + geom_point() + stat_smooth(method = "lm", se = FALSE)
## `geom_smooth()` using formula = 'y ~ x'

Test Your Knowledge

Instructions:

  1. Use tidyverse functions to create a new variable for the Student-Teacher Ratio (students/teachers). Call this variable STR.

  2. Estimate the regression models: \[ math_i = \beta_0 + \beta_1 STR_i + u_i\] and \[ read_i = \beta_0 + \beta_1 STR_i + u_i\]

  3. Summarize your regression output from the above regressions in a table.

  4. Please upload your .R file to today’s Canvas assignment.