Introduction to Econometrics with R

Florian Oswald, Jean-Marc Robin and Vincent Viers

2019-10-10
2
Contents

Syllabus 5

1 Introduction to R 9
1.1 Getting Started . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.2 Starting R and RStudio . . . . . . . . . . . . . . . . . . . . . . . 10
1.3 Basic Calculations . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.4 Getting Help . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.5 Installing Packages . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.6 Code vs Output in this Book . . . . . . . . . . . . . . . . . . . . 14
1.7 ScPoEconometrics Package . . . . . . . . . . . . . . . . . . . . . 15
1.8 Data Types . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.9 Data Structures . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.10 Data Frames . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
1.11 Programming Basics . . . . . . . . . . . . . . . . . . . . . . . . . 35

2 Working With Data 39

2.1 Summary Statistics . . . . . . . . . . . . . . . . . . . . . . . . . . 39
2.2 Plotting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
2.3 Summarizing Two Variables . . . . . . . . . . . . . . . . . . . . . 49
2.4 The tidyverse . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

3 Linear Regression 67
3.1 How are x and y related? . . . . . . . . . . . . . . . . . . . . . . 67
3.2 Ordinary Least Squares (OLS) Estimator . . . . . . . . . . . . . 73
3.3 Predictions and Residuals . . . . . . . . . . . . . . . . . . . . . . 77
3.4 Correlation, Covariance and Linearity . . . . . . . . . . . . . . . 78
3.5 Analysing V ar(y) . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
3.6 Assessing the Goodness of Fit . . . . . . . . . . . . . . . . . . . . 82
3.7 An Example: A Log Wage Equation . . . . . . . . . . . . . . . . 82
3.8 Scaling Regressions . . . . . . . . . . . . . . . . . . . . . . . . . . 86
3.9 A Particular Rescaling: The log Transform . . . . . . . . . . . . 88

4 Multiple Regression 93
4.1 All Else Equal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94

3
4 CONTENTS

4.2 Multicolinearity . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
4.3 Log Wage Equation . . . . . . . . . . . . . . . . . . . . . . . . . 97
4.4 How To Make Predictions . . . . . . . . . . . . . . . . . . . . . . 100

5 Categorial Variables 103

5.1 The Binary Regressor Case . . . . . . . . . . . . . . . . . . . . . 103
5.2 Dummy and Continuous Variables . . . . . . . . . . . . . . . . . 106
5.3 Categorical Variables in R: factor . . . . . . . . . . . . . . . . . 107
5.4 Interactions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
5.5 (Unobserved) Individual Heterogeneity . . . . . . . . . . . . . . . 112

6 Standard Errors 115

6.1 Sampling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
6.2 Taking Eleven Samples From The Population . . . . . . . . . . . 119
6.3 Handover to Moderndive . . . . . . . . . . . . . . . . . . . . . . . 121
6.4 Inference in Theory . . . . . . . . . . . . . . . . . . . . . . . . . . 121
6.5 Uncertainty in Regression Estimates . . . . . . . . . . . . . . . . 128
6.6 What is true? What are Statistical Models? . . . . . . . . . . . . 128
6.7 The Classical Regression Model (CRM) . . . . . . . . . . . . . . 129
6.8 Standard Errors in Theory . . . . . . . . . . . . . . . . . . . . . . 131
6.9 What’s in my model? (And what is not?) . . . . . . . . . . . . . 135

7 Instrumental Variables 141

7.1 Simultaneity Bias . . . . . . . . . . . . . . . . . . . . . . . . . . . 141

8 Projects 143
8.1 Trade Exercise . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143
Syllabus

Welcome to Introductory Econometrics for 2nd year undergraduates at ScPo!

On this page we outline the course and present the Syllabus. 2018/2019 is the
first time this course will be taught, so we are still in a beta release stage - you
should expect a couple of loose ends here and there, but we think the overall
experience is going to be pleasant!

Objective

This course aims to teach you the basics of data analysis needed in a Social
Sciences oriented University like SciencesPo. We purposefully start at a level
that assumes no prior knowledge about statistics whatsoever. Our objective is
to have you understand and be able to interpret linear regression analysis. We
will not rely on maths and statistics, but practical learning in order to teach
the main concepts.

5
6 CONTENTS

Syllabus and Requirements

You can find the topics we want to go over in the left panel of this page. The
later chapters are optional and depend on the speed with which we will proceed
eventually. Chapters 1-4 are the core material of the course.

The only requirement is that you bring your own personal computer to
each session. We will be using the free statistical computing language R very
intensively. Before coming to the first session, please install R and RStudio as
explained at the beginning of chapter 1.

Course Structure

This course is taught in several different groups across various campuses of

SciencesPo. All groups will go over the same material, do the same exercises,
and will have the same assessments.

Groups meet once per week for 2 hours. The main purpose of the weekly meet-
ings is to clarify any questions, and to work together through tutorials. The
little theory we need will be covered in this book, and you are expected to
read through this in your own time before coming to class.

This Book and Other Material

What you are looking at is an online textbook. You can therefore look at it in
your browser (as you are doing just now), on your mobile phone or tablet, but
you can also download it as a pdf file or as an epub file for your ebook-reader.
We don’t have any ambition to actually produce and publish a book for now, so
you should just see this as a way to disseminate our lecture notes to you. The
second part of course material next to the book is an extensive suite of tutorials
and interactive demonstrations, which are all contained in the R package that
builds this book (and which you installed by issuing the above commands).

Open Source

The book and all other content for this course are hosted under an open source
license on github. You can contribute to the book by just clicking on the appro-
priate edit symbol in the top bar of this page. Other teachers who want to use
our material can freely do so, observing the terms of the license on the github
repository.
CONTENTS 7

Assessments

We will assess participation in class and conduct a final exam.

Communication

We will communicate exclusively on our slack group. You will get an invitation
email to join from your instructor in due course.
8 CONTENTS
Chapter 1

Introduction to R

1.1 Getting Started

R is both a programming language and software environment for statistical com-

puting, which is free and open-source. To get started, you will need to install
two pieces of software:

1. R, the actual programming language.

• Chose your operating system, and select the most recent version.
2. RStudio, an excellent IDE for working with R.
• Note, you must have R installed to use RStudio. RStudio is simply
an interface used to interact with R.

The popularity of R is on the rise, and everyday it becomes a better tool for
statistical analysis. It even generated this book!

The following few chapters will serve as a whirlwind introduction to R. They are
by no means meant to be a complete reference for the R language, but simply an
introduction to the basics that we will need along the way. Several of the more
important topics will be re-stressed as they are actually needed for analyses.

This introductory R chapter may feel like an overwhelming amount of informa-

tion. You are not expected to pick up everything the first time through. You
should try all of the code from this chapter, then return to it a number of times
as you return to the concepts when performing analyses. We only present the
most basic aspects of R. If you want to know more, there are countless online
tutorials, and you could start with the official CRAN sample session or have a
look at the resources at Rstudio or on this github repo.

9
10 CHAPTER 1. INTRODUCTION TO R

Figure 1.1: R GUI symbol and R in a MacOS Terminal

Figure 1.2: R GUI symbol and R in a MacOS Terminal

1.2 Starting R and RStudio

A key difference for you to understand is the one between R, the actual program-
ming language, and RStudio, a popular interface to R which allows you to work
efficiently and with greater ease with R.
The best way to appreciate the value of RStudio is to start using R without
RStudio. To do this, double-click on the R GUI that you should have down-
loaded on your computer following the steps above (on windows or Mac), or
start R in your terminal (on Linux or Mac) by just typing R in a terminal, see
figure ??. You’ve just opened the R console which allows you to start typing
code right after the > sign, called prompt. Try typing 2 + 2 or print("Your
Name") and hit the return key. And voilà, your first R commands!
Typing one command after the other into the console is not very convenient
as our analysis becomes more involved. Ideally, we would like to collect all
1.3. BASIC CALCULATIONS 11

command statements in a file and run them one after the other, automatically.
We can do this by writing so-called script files or just scripts, i.e. simple text
files with extension .R or .r which can be inserted (or sourced) into an R session.
RStudio makes this process very easy.
Open RStudio by clicking on the RStudio application on your computer, and
notice how different the whole environment is from the basic R console – in fact,
that very same R console is running in your bottom left panel. The upper-left
panel is a space for you to write scripts – that is to say many lines of codes
which you can run when you choose to. To run a single line of code, simply
highlight it and hit Command + Return.

We highly recommend that you use RStudio for everything related to

this course (in particular, to launch our apps and tutorials).

RStudio has a large number of useful keyboard shortcuts. A list of these can be
found using a keyboard shortcut – the keyboard shortcut to rule them all:
• On Windows: Alt + Shift + K
• On Mac: Option + Shift + K
The RStudio team has developed a number of “cheatsheets” for working with
both R and RStudio. This particular cheatseet for Base R will summarize many
of the concepts in this document. 1

1.2.1 First Glossary

• R: a statistical programming language

• RStudio: an integrated development environment (IDE) to work with R
• command: user input (text or numbers) that R understands.
• script: a list of commands collected in a text file, each separated by a new
line, to be run one after the other.

1.3 Basic Calculations

To get started, we’ll use R like a simple calculator. Run the following code either
directly from your RStudio console, or in RStudio by writting them in a script
and running them using Command + Return.
1 When programming, it is often a good practice to follow a style guide. (Where do spaces

go? Tabs or spaces? Underscores or CamelCase when naming variables?) No style guide is
“correct” but it helps to be aware of what others do. The more import thing is to be consistent
within your own code. Here are two guides: Hadley Wickham Style Guide, and the Google
Style Guide. For this course, our main deviation from these two guides is the use of = in place
of <-. For all practical purposes, you should think = whenever you see <-.
12 CHAPTER 1. INTRODUCTION TO R

Addition, Subtraction, Multiplication and Division

Math R code Result

3+2 3 + 2 5
3−2 3 - 2 1
3·2 3 * 2 6
3/2 3 / 2 1.5

Exponents

Math R code Result

32 3 ^ 2 9
(−3)
2 2 ^ (-3) 0.125
1/2
100
√ 100 ^ (1 / 2) 10
100 sqrt(100) 10

Mathematical Constants

Math R code Result

π pi 3.1415927
e exp(1) 2.7182818

Logarithms

Note that we will use ln and log interchangeably to mean the natural logarithm.
There is no ln() in R, instead it uses log() to mean the natural logarithm.

Math R code Result

log(e) log(exp(1)) 1
log10 (1000) log10(1000) 3
log2 (8) log2(8) 3
log4 (16) log(16, base = 4) 2

Trigonometry

Math R code Result

sin(π/2) sin(pi / 2) 1
1.4. GETTING HELP 13

Math R code Result

cos(0) cos(0) 1

1.4 Getting Help

In using R as a calculator, we have seen a number of functions: sqrt(),
exp(), log() and sin(). To get documentation about a function in R,
simply put a question mark in front of the function name, or call the function
help(function) and RStudio will display the documentation, for example:
?log
?sin
?paste
?lm
help(lm) # help() is equivalent
help(ggplot,package="ggplot2") # show help from a certain package

Frequently one of the most difficult things to do when learning R is asking for
help. First, you need to decide to ask for help, then you need to know how
to ask for help. Your very first line of defense should be to Google your error
message or a short description of your issue. (The ability to solve problems
using this method is quickly becoming an extremely valuable skill.) If that fails,
and it eventually will, you should ask for help. There are a number of things
you should include when contacting an instructor, or posting to a help website
such as Stack Overflow.
• Describe what you expect the code to do.
• State the end goal you are trying to achieve. (Sometimes what you expect
the code to do, is not what you want to actually do.)
• Provide the full text of any errors you have received.
• Provide enough code to recreate the error. Often for the purpose of this
course, you could simply post your entire .R script or .Rmd to slack.
• Sometimes it is also helpful to include a screenshot of your entire RStudio
window when the error occurs.
If you follow these steps, you will get your issue resolved much quicker, and
possibly learn more in the process. Do not be discouraged by running into
errors and difficulties when learning R. (Or any other technical skill.) It is
simply part of the learning process.

1.5 Installing Packages

R comes with a number of built-in functions and datasets, but one of the main
strengths of R as an open-source project is its package system. Packages add
14 CHAPTER 1. INTRODUCTION TO R

additional functions and data. Frequently if you want to do something in R, and

it is not available by default, there is a good chance that there is a package that
will fulfill your needs.
To install a package, use the [Link]() function. Think of this as
buying a recipe book from the store, bringing it home, and putting it on your
shelf (i.e. into your library):
[Link]("ggplot2")

Once a package is installed, it must be loaded into your current R session before
being used. Think of this as taking the book off of the shelf and opening it up
to read.
library(ggplot2)

Once you close R, all the packages are closed and put back on the imaginary
shelf. The next time you open R, you do not have to install the package again,
but you do have to load any packages you intend to use by invoking library().

1.6 Code vs Output in this Book

A quick note on styling choices in this book. We had to make a decision how to
visually separate R code and resulting output in this book. We decided to prefix
all output lines with #OUT> to make the distinction. A typical code snippet with
output is thus going to look like this:
1 + 3

#OUT> [1] 4
# everything after a # is a comment, i.e. R disregards it.

where you see on the first line the R code, and on the second line the output.
As mentioned, that line starts with #OUT> to say this is an output, followed by
[1] (indicating this is a vector of length one - more on this below!), followed by
the actual result - 1 + 3 = 4!
Notice that you can simply copy and paste all the code you see into your R
console. In fact, you are strongly encouraged to actually do this and try out all
the code you see in this book.
Finally, please note that this way of showing output is fully our choice in this
textbook, and that you should expect other output formats elsewhere. For
example, in my RStudio console, the above code and output looks like this:
> 1 + 3
[1] 4
1.7. SCPOECONOMETRICS PACKAGE 15

1.7 ScPoEconometrics Package

To fully take advantage of our course, please install the associated R package
directly from its online code repository. You can do this by copy and pasting
the following three lines into your R console:
if (!require("devtools")) [Link]("devtools")
library(devtools)
install_github(repo = "ScPoEcon/ScPoEconometrics")

In order to check whether everything works fine, you could load the library, and
check it’s current version:
library(ScPoEconometrics)
packageVersion("ScPoEconometrics")

#OUT> [1] '0.2.4'

1.8 Data Types

R has a number of basic data types. While R is not a strongly typed language
(i.e. you can be agnostic about types most of the times), it is useful to know
what data types are available to you:

• Numeric
– Also known as Double. The default type when dealing with numbers.
– Examples: 1, 1.0, 42.5
• Integer
– Examples: 1L, 2L, 42L
• Complex
– Example: 4 + 2i
• Logical
– Two possible values: TRUE and FALSE
– You can also use T and F, but this is not recommended.
– NA is also considered logical.
• Character
– Examples: "a", "Statistics", "1 plus 2."
• Categorical or factor
– A mixture of integer and character. A factor variable assigns a label
to a numeric value.
– For example factor(x=c(0,1),labels=c("male","female"))
assigns the string male to the numeric values 0, and the string
female to the value 1.
16 CHAPTER 1. INTRODUCTION TO R

1.9 Data Structures

R also has a number of basic data structures. A data structure is either homoge-
neous (all elements are of the same data type) or heterogeneous (elements can
be of more than one data type).

Dimension Homogeneous Heterogeneous

1 Vector List
2 Matrix Data Frame
3+ Array nested Lists

1.9.1 Vectors

Many operations in R make heavy use of vectors. A vector is a container for

objects of identical type (see 1.8 above). Vectors in R are indexed starting at
1. That is what the [1] in the output is indicating, that the first element of
the row being displayed is the first element of the vector. Larger vectors will
start additional rows with something like [7] where 7 is the index of the first
element of that row.
Possibly the most common way to create a vector in R is using the c() function,
which is short for “combine”. As the name suggests, it combines a list of elements
separated by commas. (Are you busy typing all of those examples into your R
console? :-) )
c(1, 3, 5, 7, 8, 9)

#OUT> [1] 1 3 5 7 8 9
Here R simply outputs this vector. If we would like to store this vector in
a variable we can do so with the assignment operator =. In this case the
variable x now holds the vector we just created, and we can access the vector
by typing x.
x = c(1, 3, 5, 7, 8, 9)
x

#OUT> [1] 1 3 5 7 8 9
As an aside, there is a long history of the assignment operator in R, partially
due to the keys available on the keyboards of the creators of the S language.
(Which preceded R.) For simplicity we will use =, but know that often you will
see <- as the assignment operator.
Because vectors must contain elements that are all the same type, R will au-
tomatically coerce (i.e. convert) to a single type when attempting to create a
vector that combines multiple types.
1.9. DATA STRUCTURES 17

c(42, "Statistics", TRUE)

#OUT> [1] "42" "Statistics" "TRUE"

c(42, TRUE)

#OUT> [1] 42 1
Frequently you may wish to create a vector based on a sequence of numbers.
The quickest and easiest way to do this is with the : operator, which creates a
sequence of integers between two specified integers.
(y = 1:100)

#OUT> [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
#OUT> [18] 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
#OUT> [35] 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51
#OUT> [52] 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68
#OUT> [69] 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85
#OUT> [86] 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
Here we see R labeling the rows after the first since this is a large vector. Also,
we see that by putting parentheses around the assignment, R both stores the
vector in a variable called y and automatically outputs y to the console.
Note that scalars do not exists in R. They are simply vectors of length 1.
2

#OUT> [1] 2
If we want to create a sequence that isn’t limited to integers and increasing by
1 at a time, we can use the seq() function.
seq(from = 1.5, to = 4.2, by = 0.1)

#OUT> [1] 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 3.1
#OUT> [18] 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 4.0 4.1 4.2
We will discuss functions in detail later, but note here that the input labels
from, to, and by are optional.
seq(1.5, 4.2, 0.1)

#OUT> [1] 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 3.1
#OUT> [18] 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 4.0 4.1 4.2
Another common operation to create a vector is rep(), which can repeat a
single value a number of times.
rep("A", times = 10)

#OUT> [1] "A" "A" "A" "A" "A" "A" "A" "A" "A" "A"
18 CHAPTER 1. INTRODUCTION TO R

The rep() function can be used to repeat a vector some number of times.
rep(x, times = 3)

#OUT> [1] 1 3 5 7 8 9 1 3 5 7 8 9 1 3 5 7 8 9
We have now seen four different ways to create vectors:
• c()
• :
• seq()
• rep()
So far we have mostly used them in isolation, but they are often used together.
c(x, rep(seq(1, 9, 2), 3), c(1, 2, 3), 42, 2:4)

#OUT> [1] 1 3 5 7 8 9 1 3 5 7 9 1 3 5 7 9 1 3 5 7 9 1 2
#OUT> [24] 3 42 2 3 4
The length of a vector can be obtained with the length() function.
length(x)

#OUT> [1] 6
length(y)

#OUT> [1] 100

Let’s try this out! Your turn:

[Link] Task 1

1. Create a vector of five ones, i.e. [1,1,1,1,1]

2. Notice that the colon operator a:b is just short for construct a sequence
from a to b. Create a vector the counts down from 10 to 0, i.e. it looks
like [10,9,8,7,6,5,4,3,2,1,0]!
3. the rep function takes additional arguments times (as above), and each,
which tells you how often each element should be repeated (as opposed to
the entire input vector). Use rep to create a vector that looks like this:
[1 1 1 2 2 2 3 3 3 1 1 1 2 2 2 3 3 3]

[Link] Subsetting

To subset a vector, i.e. to choose only some elements of it, we use square brackets,
[]. Here we see that x[1] returns the first element, and x[3] returns the third
element:
1.9. DATA STRUCTURES 19

#OUT> [1] 1 3 5 7 8 9
x[1]

#OUT> [1] 1
x[3]

#OUT> [1] 5
We can also exclude certain indexes, in this case the second element.
x[-2]

#OUT> [1] 1 5 7 8 9
Lastly we see that we can subset based on a vector of indices.
x[1:3]

#OUT> [1] 1 3 5
x[c(1,3,4)]

#OUT> [1] 1 5 7
All of the above are subsetting a vector using a vector of indexes. (Remember a
single number is still a vector.) We could instead use a vector of logical values.
z = c(TRUE, TRUE, FALSE, TRUE, TRUE, FALSE)
z

#OUT> [1] TRUE TRUE FALSE TRUE TRUE FALSE

x[z]

#OUT> [1] 1 3 7 8
R is able to perform many operations on vectors and scalars alike:
x = 1:10 # a vector
x + 1 # add a scalar

#OUT> [1] 2 3 4 5 6 7 8 9 10 11
2 * x # multiply all elements by 2

#OUT> [1] 2 4 6 8 10 12 14 16 18 20
2 ^ x # take 2 to the x as exponents

#OUT> [1] 2 4 8 16 32 64 128 256 512 1024

sqrt(x) # compute the square root of all elements in x
20 CHAPTER 1. INTRODUCTION TO R

#OUT> [1] 1.000000 1.414214 1.732051 2.000000 2.236068 2.449490 2.645751

#OUT> [8] 2.828427 3.000000 3.162278
log(x) # take the natural log of all elements in x

#OUT> [1] 0.0000000 0.6931472 1.0986123 1.3862944 1.6094379 1.7917595 1.9459101

#OUT> [8] 2.0794415 2.1972246 2.3025851
x + 2*x # add vector x to vector 2x

#OUT> [1] 3 6 9 12 15 18 21 24 27 30
We see that when a function like log() is called on a vector x, a vector is
returned which has applied the function to each element of the vector x.

1.9.2 Logical Operators

Operator Summary Example Result

x < y x less than y 3 < 42 TRUE
x > y x greater than y 3 > 42 FALSE
x <= y x less than or equal to y 3 <= 42 TRUE
x >= y x greater than or equal to y 3 >= 42 FALSE
x == y xequal to y 3 == 42 FALSE
x != y x not equal to y 3 != 42 TRUE
!x not x !(3 > 42) TRUE
x | y x or y (3 > 42) | TRUE TRUE
x & y x and y (3 < 4) & ( 42 > 13) TRUE

In R, logical operators also work on vectors:

x = c(1, 3, 5, 7, 8, 9)
x > 3

#OUT> [1] FALSE FALSE TRUE TRUE TRUE TRUE

x < 3

#OUT> [1] TRUE FALSE FALSE FALSE FALSE FALSE

x == 3

#OUT> [1] FALSE TRUE FALSE FALSE FALSE FALSE

x != 3

#OUT> [1] TRUE FALSE TRUE TRUE TRUE TRUE

1.9. DATA STRUCTURES 21

x == 3 & x != 3

#OUT> [1] FALSE FALSE FALSE FALSE FALSE FALSE

x == 3 | x != 3

#OUT> [1] TRUE TRUE TRUE TRUE TRUE TRUE

This is quite useful for subsetting.
x[x > 3]

#OUT> [1] 5 7 8 9
x[x != 3]

#OUT> [1] 1 5 7 8 9
sum(x > 3)

#OUT> [1] 4
[Link](x > 3)

#OUT> [1] 0 0 1 1 1 1
Here we saw that using the sum() function on a vector of logical TRUE and
FALSE values that is the result of x > 3 results in a numeric result: you just
counted for how many elements of x, the condition > 3 is TRUE. During the call
to sum(), R is first automatically coercing the logical to numeric where TRUE
is 1 and FALSE is 0. This coercion from logical to numeric happens for most
mathematical operations.
# which(condition of x) returns true/false
# each index of x where condition is true
which(x > 3)

#OUT> [1] 3 4 5 6
x[which(x > 3)]

#OUT> [1] 5 7 8 9
max(x)

#OUT> [1] 9
which(x == max(x))

#OUT> [1] 6
[Link](x)

#OUT> [1] 6
22 CHAPTER 1. INTRODUCTION TO R

[Link] Task 2

1. Create a vector filled with 10 numbers drawn from the uniform distribution
(hint: use function runif) and store them in x.
2. Using logical subsetting as above, get all the elements of x which are larger
than 0.5, and store them in y.
3. using the function which, store the indices of all the elements of x which
are larger than 0.5 in iy.
4. Check that y and x[iy] are identical.

1.9.3 Matrices

R can also be used for matrix calculations. Matrices have rows and columns
containing a single data type. In a matrix, the order of rows and columns is
important. (This is not true of data frames, which we will see later.)
Matrices can be created using the matrix function.
x = 1:9
x

#OUT> [1] 1 2 3 4 5 6 7 8 9
X = matrix(x, nrow = 3, ncol = 3)
X

#OUT> [,1] [,2] [,3]

#OUT> [1,] 1 4 7
#OUT> [2,] 2 5 8
#OUT> [3,] 3 6 9
Notice here that R is case sensitive (x vs X).
By default the matrix function fills your data into the matrix column by column.
But we can also tell R to fill rows instead:
Y = matrix(x, nrow = 3, ncol = 3, byrow = TRUE)
Y

#OUT> [,1] [,2] [,3]

#OUT> [1,] 1 2 3
#OUT> [2,] 4 5 6
#OUT> [3,] 7 8 9
We can also create a matrix of a specified dimension where every element is the
same, in this case 0.
Z = matrix(0, 2, 4)
Z
1.9. DATA STRUCTURES 23

#OUT> [,1] [,2] [,3] [,4]

#OUT> [1,] 0 0 0 0
#OUT> [2,] 0 0 0 0
Like vectors, matrices can be subsetted using square brackets, []. However,
since matrices are two-dimensional, we need to specify both a row and a column
when subsetting.
X

#OUT> [,1] [,2] [,3]

#OUT> [1,] 1 4 7
#OUT> [2,] 2 5 8
#OUT> [3,] 3 6 9
X[1, 2]

#OUT> [1] 4
Here we accessed the element in the first row and the second column. We could
also subset an entire row or column.
X[1, ]

#OUT> [1] 1 4 7
X[, 2]

#OUT> [1] 4 5 6
We can also use vectors to subset more than one row or column at a time. Here
we subset to the first and third column of the second row:
X[2, c(1, 3)]

#OUT> [1] 2 8
Matrices can also be created by combining vectors as columns, using cbind, or
combining vectors as rows, using rbind.
x = 1:9
rev(x)

#OUT> [1] 9 8 7 6 5 4 3 2 1
rep(1, 9)

#OUT> [1] 1 1 1 1 1 1 1 1 1
rbind(x, rev(x), rep(1, 9))

#OUT> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
#OUT> x 1 2 3 4 5 6 7 8 9
#OUT> 9 8 7 6 5 4 3 2 1
24 CHAPTER 1. INTRODUCTION TO R

#OUT> 1 1 1 1 1 1 1 1 1
cbind(col_1 = x, col_2 = rev(x), col_3 = rep(1, 9))

#OUT> col_1 col_2 col_3

#OUT> [1,] 1 9 1
#OUT> [2,] 2 8 1
#OUT> [3,] 3 7 1
#OUT> [4,] 4 6 1
#OUT> [5,] 5 5 1
#OUT> [6,] 6 4 1
#OUT> [7,] 7 3 1
#OUT> [8,] 8 2 1
#OUT> [9,] 9 1 1
When using rbind and cbind you can specify “argument” names that will be
used as column names.
R can then be used to perform matrix calculations.
x = 1:9
y = 9:1
X = matrix(x, 3, 3)
Y = matrix(y, 3, 3)
X

#OUT> [,1] [,2] [,3]

#OUT> [1,] 1 4 7
#OUT> [2,] 2 5 8
#OUT> [3,] 3 6 9
Y

#OUT> [,1] [,2] [,3]

#OUT> [1,] 9 6 3
#OUT> [2,] 8 5 2
#OUT> [3,] 7 4 1
X + Y

#OUT> [,1] [,2] [,3]

#OUT> [1,] 10 10 10
#OUT> [2,] 10 10 10
#OUT> [3,] 10 10 10
X - Y

#OUT> [,1] [,2] [,3]

#OUT> [1,] -8 -2 4
#OUT> [2,] -6 0 6
#OUT> [3,] -4 2 8
1.9. DATA STRUCTURES 25

X * Y

#OUT> [,1] [,2] [,3]

#OUT> [1,] 9 24 21
#OUT> [2,] 16 25 16
#OUT> [3,] 21 24 9
X / Y

#OUT> [,1] [,2] [,3]

#OUT> [1,] 0.1111111 0.6666667 2.333333
#OUT> [2,] 0.2500000 1.0000000 4.000000
#OUT> [3,] 0.4285714 1.5000000 9.000000
Note that X * Y is not matrix multiplication. It is element by element mul-
tiplication. (Same for X / Y). Matrix multiplication uses %*%. Other matrix
functions include t() which gives the transpose of a matrix and solve() which
returns the inverse of a square matrix if it is invertible.
X %*% Y

#OUT> [,1] [,2] [,3]

#OUT> [1,] 90 54 18
#OUT> [2,] 114 69 24
#OUT> [3,] 138 84 30
t(X)

#OUT> [,1] [,2] [,3]

#OUT> [1,] 1 2 3
#OUT> [2,] 4 5 6
#OUT> [3,] 7 8 9

1.9.4 Arrays

A vector is a one-dimensional array. A matrix is a two-dimensional array. In R

you can create arrays of arbitrary dimensionality N. Here is how:
d = 1:16
d3 = array(data = d,dim = c(4,2,2))
d4 = array(data = d,dim = c(4,2,2,3)) # will recycle 1:16
d3

#OUT> , , 1
#OUT>
#OUT> [,1] [,2]
#OUT> [1,] 1 5
#OUT> [2,] 2 6
#OUT> [3,] 3 7
26 CHAPTER 1. INTRODUCTION TO R

#OUT> [4,] 4 8
#OUT>
#OUT> , , 2
#OUT>
#OUT> [,1] [,2]
#OUT> [1,] 9 13
#OUT> [2,] 10 14
#OUT> [3,] 11 15
#OUT> [4,] 12 16

You can see that d3 are simply two (4,2) matrices laid on top of each other,
as if there were two pages. Similary, d4 would have two pages, and another 3
registers in a fourth dimension. And so on. You can subset an array like you
would a vector or a matrix, taking care to index each dimension:
d3[ ,1,1] # all elements from col 1, page 1

#OUT> [1] 1 2 3 4
d3[2:3, , ] # rows 2:3 from all pages

#OUT> , , 1
#OUT>
#OUT> [,1] [,2]
#OUT> [1,] 2 6
#OUT> [2,] 3 7
#OUT>
#OUT> , , 2
#OUT>
#OUT> [,1] [,2]
#OUT> [1,] 10 14
#OUT> [2,] 11 15
d3[2,2, ] # row 2, col 2 from both pages.

#OUT> [1] 6 14

[Link] Task 3

1. Create a vector containing 1,2,3,4,5 called v.

2. Create a (2,5) matrix m containing the data 1,2,3,4,5,6,7,8,9,10. The
first row should be 1,2,3,4,5.
3. Perform matrix multiplication of m with v. Use the command %*%. What
dimension does the output have?
4. Why does v %*% m not work?
1.9. DATA STRUCTURES 27

1.9.5 Lists

A list is a one-dimensional heterogeneous data structure. So it is indexed like

a vector with a single integer value (or with a name), but each element can
contain an element of any type. Lists are similar to a python or julia Dict
object. Many R structures and outputs are lists themselves. Lists are extremely
useful and versatile objects, so make sure you understand their useage:
# creation without fieldnames
list(42, "Hello", TRUE)

#OUT> [[1]]
#OUT> [1] 42
#OUT>
#OUT> [[2]]
#OUT> [1] "Hello"
#OUT>
#OUT> [[3]]
#OUT> [1] TRUE
# creation with fieldnames
ex_list = list(
a = c(1, 2, 3, 4),
b = TRUE,
c = "Hello!",
d = function(arg = 42) {print("Hello World!")},
e = diag(5)
)

Lists can be subset using two syntaxes, the $ operator, and square brackets [].
The $ operator returns a named element of a list. The [] syntax returns a
list, while the [[]] returns an element of a list.
• ex_list[1] returns a list contain the first element.
• ex_list[[1]] returns the first element of the list, in this case, a vector.
# subsetting
ex_list$e

#OUT> [,1] [,2] [,3] [,4] [,5]

#OUT> [1,] 1 0 0 0 0
#OUT> [2,] 0 1 0 0 0
#OUT> [3,] 0 0 1 0 0
#OUT> [4,] 0 0 0 1 0
#OUT> [5,] 0 0 0 0 1
ex_list[1:2]

#OUT> $a
28 CHAPTER 1. INTRODUCTION TO R

#OUT> [1] 1 2 3 4
#OUT>
#OUT> $b
#OUT> [1] TRUE
ex_list[1]

#OUT> $a
#OUT> [1] 1 2 3 4
ex_list[[1]]

#OUT> [1] 1 2 3 4
ex_list[c("e", "a")]

#OUT> $e
#OUT> [,1] [,2] [,3] [,4] [,5]
#OUT> [1,] 1 0 0 0 0
#OUT> [2,] 0 1 0 0 0
#OUT> [3,] 0 0 1 0 0
#OUT> [4,] 0 0 0 1 0
#OUT> [5,] 0 0 0 0 1
#OUT>
#OUT> $a
#OUT> [1] 1 2 3 4
ex_list["e"]

#OUT> $e
#OUT> [,1] [,2] [,3] [,4] [,5]
#OUT> [1,] 1 0 0 0 0
#OUT> [2,] 0 1 0 0 0
#OUT> [3,] 0 0 1 0 0
#OUT> [4,] 0 0 0 1 0
#OUT> [5,] 0 0 0 0 1
ex_list[["e"]]

#OUT> [,1] [,2] [,3] [,4] [,5]

#OUT> [1,] 1 0 0 0 0
#OUT> [2,] 0 1 0 0 0
#OUT> [3,] 0 0 1 0 0
#OUT> [4,] 0 0 0 1 0
#OUT> [5,] 0 0 0 0 1
ex_list$d

#OUT> function(arg = 42) {print("Hello World!")}

#OUT> <environment: 0x7fc5d141ee78>
1.10. DATA FRAMES 29

ex_list$d(arg = 1)

#OUT> [1] "Hello World!"

[Link] Task 4

1. Copy and paste the above code for ex_list into your R session. Re-
member that list can hold any kind of R object. Like…another list! So,
create a new list new_list that has two fields: a first field called “this”
with string content "is awesome", and a second field called “ex_list” that
contains ex_list.
2. Accessing members is like in a plain list, just with several layers now. Get
the element c from ex_list in new_list!
3. Compose a new string out of the first element in new_list, the element
under label this. Use the function paste to print R is awesome to your
screen.

1.10 Data Frames

We have previously seen vectors and matrices for storing data as we introduced
R. We will now introduce a data frame which will be the most common way
that we store and interact with data in this course. A [Link] is similar to
a python [Link] or a julia DataFrame. (But the R version was the
first! :-) )
example_data = [Link](x = c(1, 3, 5, 7, 9, 1, 3, 5, 7, 9),
y = c(rep("Hello", 9), "Goodbye"),
z = rep(c(TRUE, FALSE), 5))

Unlike a matrix, which can be thought of as a vector rearranged into rows and
columns, a data frame is not required to have the same data type for each
element. A data frame is a list of vectors, and each vector has a name. So,
each vector must contain the same data type, but the different vectors can
store different data types. Note, however, that all vectors must have the same
length (differently from a list)!

A [Link] is similar to a typical Spreadsheet. There are rows, and

there are columns. A row is typically thought of as an observation, and
each column is a certain variable, characteristic or feature of that obser-
vation.

Let’s look at the data frame we just created above:

30 CHAPTER 1. INTRODUCTION TO R

example_data

#OUT> x y z
#OUT> 1 1 Hello TRUE
#OUT> 2 3 Hello FALSE
#OUT> 3 5 Hello TRUE
#OUT> 4 7 Hello FALSE
#OUT> 5 9 Hello TRUE
#OUT> 6 1 Hello FALSE
#OUT> 7 3 Hello TRUE
#OUT> 8 5 Hello FALSE
#OUT> 9 7 Hello TRUE
#OUT> 10 9 Goodbye FALSE

Unlike a list, which has more flexibility, the elements of a data frame must all
be vectors. Again, we access any given column with the $ operator:
example_data$x

#OUT> [1] 1 3 5 7 9 1 3 5 7 9
[Link](length(example_data$x),
length(example_data$y),
length(example_data$z))

#OUT> [1] TRUE

str(example_data)

#OUT> '[Link]': 10 obs. of 3 variables:

#OUT> $ x: num 1 3 5 7 9 1 3 5 7 9
#OUT> $ y: Factor w/ 2 levels "Goodbye","Hello": 2 2 2 2 2 2 2 2 2 1
#OUT> $ z: logi TRUE FALSE TRUE FALSE TRUE FALSE ...
nrow(example_data)

#OUT> [1] 10
ncol(example_data)

#OUT> [1] 3
dim(example_data)

#OUT> [1] 10 3
names(example_data)

#OUT> [1] "x" "y" "z"

1.10. DATA FRAMES 31

1.10.1 Working with [Link]

The [Link]() function above is one way to create a data frame. We can
also import data from various file types in into R, as well as use data stored in
packages.
To read this data back into R, we will use the built-in function [Link]:
path = [Link](package="ScPoEconometrics","datasets","[Link]")
example_data_from_disk = [Link](path)

This particular line of code assumes that you installed the associated R
package to this book, hence you have this dataset stored on your computer at
[Link](package = "ScPoEconometrics","datasets","[Link]").
example_data_from_disk

#OUT> x y z
#OUT> 1 1 Hello TRUE
#OUT> 2 3 Hello FALSE
#OUT> 3 5 Hello TRUE
#OUT> 4 7 Hello FALSE
#OUT> 5 9 Hello TRUE
#OUT> 6 1 Hello FALSE
#OUT> 7 3 Hello TRUE
#OUT> 8 5 Hello FALSE
#OUT> 9 7 Hello TRUE
#OUT> 10 9 Goodbye FALSE
When using data, there are three things we would generally like to do:
• Look at the raw data.
• Understand the data. (Where did it come from? What are the variables?
Etc.)
• Visualize the data.
To look at data in a [Link], we have two useful commands: head() and
str().
# we are working with the built-in mtcars dataset:
mtcars

#OUT> mpg cyl disp hp drat wt qsec vs am gear carb

#OUT> Mazda RX4 21.0 6 160.0 110 3.90 2.620 16.46 0 1 4 4
#OUT> Mazda RX4 Wag 21.0 6 160.0 110 3.90 2.875 17.02 0 1 4 4
#OUT> Datsun 710 22.8 4 108.0 93 3.85 2.320 18.61 1 1 4 1
#OUT> Hornet 4 Drive 21.4 6 258.0 110 3.08 3.215 19.44 1 0 3 1
#OUT> Hornet Sportabout 18.7 8 360.0 175 3.15 3.440 17.02 0 0 3 2
#OUT> Valiant 18.1 6 225.0 105 2.76 3.460 20.22 1 0 3 1
#OUT> Duster 360 14.3 8 360.0 245 3.21 3.570 15.84 0 0 3 4
32 CHAPTER 1. INTRODUCTION TO R

#OUT> Merc 240D 24.4 4 146.7 62 3.69 3.190 20.00 1 0 4 2

#OUT> Merc 230 22.8 4 140.8 95 3.92 3.150 22.90 1 0 4 2
#OUT> Merc 280 19.2 6 167.6 123 3.92 3.440 18.30 1 0 4 4
#OUT> Merc 280C 17.8 6 167.6 123 3.92 3.440 18.90 1 0 4 4
#OUT> Merc 450SE 16.4 8 275.8 180 3.07 4.070 17.40 0 0 3 3
#OUT> Merc 450SL 17.3 8 275.8 180 3.07 3.730 17.60 0 0 3 3
#OUT> Merc 450SLC 15.2 8 275.8 180 3.07 3.780 18.00 0 0 3 3
#OUT> Cadillac Fleetwood 10.4 8 472.0 205 2.93 5.250 17.98 0 0 3 4
#OUT> Lincoln Continental 10.4 8 460.0 215 3.00 5.424 17.82 0 0 3 4
#OUT> Chrysler Imperial 14.7 8 440.0 230 3.23 5.345 17.42 0 0 3 4
#OUT> Fiat 128 32.4 4 78.7 66 4.08 2.200 19.47 1 1 4 1
#OUT> Honda Civic 30.4 4 75.7 52 4.93 1.615 18.52 1 1 4 2
#OUT> Toyota Corolla 33.9 4 71.1 65 4.22 1.835 19.90 1 1 4 1
#OUT> Toyota Corona 21.5 4 120.1 97 3.70 2.465 20.01 1 0 3 1
#OUT> Dodge Challenger 15.5 8 318.0 150 2.76 3.520 16.87 0 0 3 2
#OUT> AMC Javelin 15.2 8 304.0 150 3.15 3.435 17.30 0 0 3 2
#OUT> Camaro Z28 13.3 8 350.0 245 3.73 3.840 15.41 0 0 3 4
#OUT> Pontiac Firebird 19.2 8 400.0 175 3.08 3.845 17.05 0 0 3 2
#OUT> Fiat X1-9 27.3 4 79.0 66 4.08 1.935 18.90 1 1 4 1
#OUT> Porsche 914-2 26.0 4 120.3 91 4.43 2.140 16.70 0 1 5 2
#OUT> Lotus Europa 30.4 4 95.1 113 3.77 1.513 16.90 1 1 5 2
#OUT> Ford Pantera L 15.8 8 351.0 264 4.22 3.170 14.50 0 1 5 4
#OUT> Ferrari Dino 19.7 6 145.0 175 3.62 2.770 15.50 0 1 5 6
#OUT> Maserati Bora 15.0 8 301.0 335 3.54 3.570 14.60 0 1 5 8
#OUT> Volvo 142E 21.4 4 121.0 109 4.11 2.780 18.60 1 1 4 2
You can see that this prints the entire [Link] to screen. The function head()
will display the first n observations of the data frame.
head(mtcars,n=2)

#OUT> mpg cyl disp hp drat wt qsec vs am gear carb

#OUT> Mazda RX4 21 6 160 110 3.9 2.620 16.46 0 1 4 4
#OUT> Mazda RX4 Wag 21 6 160 110 3.9 2.875 17.02 0 1 4 4
head(mtcars) # default

#OUT> mpg cyl disp hp drat wt qsec vs am gear carb

#OUT> Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4
#OUT> Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4
#OUT> Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1
#OUT> Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1
#OUT> Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2
#OUT> Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1
The function str() will display the “structure” of the data frame. It will display
the number of observations and variables, list the variables, give the type of
each variable, and show some elements of each variable. This information can
1.10. DATA FRAMES 33

also be found in the “Environment” window in RStudio.

str(mtcars)

#OUT> '[Link]': 32 obs. of 11 variables:

#OUT> $ mpg : num 21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
#OUT> $ cyl : num 6 6 4 6 8 6 8 4 4 6 ...
#OUT> $ disp: num 160 160 108 258 360 ...
#OUT> $ hp : num 110 110 93 110 175 105 245 62 95 123 ...
#OUT> $ drat: num 3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
#OUT> $ wt : num 2.62 2.88 2.32 3.21 3.44 ...
#OUT> $ qsec: num 16.5 17 18.6 19.4 17 ...
#OUT> $ vs : num 0 0 1 1 0 1 0 1 1 1 ...
#OUT> $ am : num 1 1 1 0 0 0 0 0 0 0 ...
#OUT> $ gear: num 4 4 4 3 3 3 3 4 4 4 ...
#OUT> $ carb: num 4 4 1 1 2 1 4 2 2 4 ...

In this dataset an observation is for a particular model of a car, and the variables
describe attributes of the car, for example its fuel efficiency, or its weight.

To understand more about the data set, we use the ? operator to pull up the
documentation for the data.
?mtcars

R has a number of functions for quickly working with and extracting basic in-
formation from data frames. To quickly obtain a vector of the variable names,
we use the names() function.
names(mtcars)

#OUT> [1] "mpg" "cyl" "disp" "hp" "drat" "wt" "qsec" "vs" "am" "gear"
#OUT> [11] "carb"

To access one of the variables as a vector, we use the $ operator.

mtcars$mpg

#OUT> [1] 21.0 21.0 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 17.8 16.4 17.3 15.2
#OUT> [15] 10.4 10.4 14.7 32.4 30.4 33.9 21.5 15.5 15.2 13.3 19.2 27.3 26.0 30.4
#OUT> [29] 15.8 19.7 15.0 21.4
mtcars$wt

#OUT> [1] 2.620 2.875 2.320 3.215 3.440 3.460 3.570 3.190 3.150 3.440 3.440
#OUT> [12] 4.070 3.730 3.780 5.250 5.424 5.345 2.200 1.615 1.835 2.465 3.520
#OUT> [23] 3.435 3.840 3.845 1.935 2.140 1.513 3.170 2.770 3.570 2.780

We can use the dim(), nrow() and ncol() functions to obtain information
about the dimension of the data frame.
34 CHAPTER 1. INTRODUCTION TO R

dim(mtcars)

#OUT> [1] 32 11
nrow(mtcars)

#OUT> [1] 32
ncol(mtcars)

#OUT> [1] 11
Here nrow() is also the number of observations, which in most cases is the
sample size.
Subsetting data frames can work much like subsetting matrices using square
brackets, [ , ]. Here, we find vehicles with mpg over 25 miles per gallon and
only display columns cyl, disp and wt.
# mpg[row condition, col condition]
mtcars[mtcars$mpg > 20, c("cyl", "disp", "wt")]

#OUT> cyl disp wt

#OUT> Mazda RX4 6 160.0 2.620
#OUT> Mazda RX4 Wag 6 160.0 2.875
#OUT> Datsun 710 4 108.0 2.320
#OUT> Hornet 4 Drive 6 258.0 3.215
#OUT> Merc 240D 4 146.7 3.190
#OUT> Merc 230 4 140.8 3.150
#OUT> Fiat 128 4 78.7 2.200
#OUT> Honda Civic 4 75.7 1.615
#OUT> Toyota Corolla 4 71.1 1.835
#OUT> Toyota Corona 4 120.1 2.465
#OUT> Fiat X1-9 4 79.0 1.935
#OUT> Porsche 914-2 4 120.3 2.140
#OUT> Lotus Europa 4 95.1 1.513
#OUT> Volvo 142E 4 121.0 2.780
An alternative would be to use the subset() function, which has a much more
readable syntax.
subset(mtcars, subset = mpg > 25, select = c("cyl", "disp", "wt"))

[Link] Task 5

1. How many observations are there in mtcars?

2. How many variables?
3. What is the average value of mpg?
1.11. PROGRAMMING BASICS 35

4. What is the average value of mpg for cars with more than 4 cylinders,
i.e. with cyl>4?

1.11 Programming Basics

In this section we illustrate some general concepts related to programming.

1.11.1 Variables

We encountered the term variable already several times, but mainly in the con-
text of a column of a [Link]. In programming, a variable is denotes an object.
Another way to say it is that a variable is a name or a label for something:
x = 1
y = "roses"
z = function(x){sqrt(x)}

Here x refers to the value √1, y holds the string “roses”, and z is the name of
a function that computes x. Notice that the argument x of the function is
different from the x we just defined. It is local to the function:
x

#OUT> [1] 1
z(9)

#OUT> [1] 3

1.11.2 Control Flow

Control Flow relates to ways in which you can adapt your code to different
circumstances. Based on a condition being TRUE, your program will do one
thing, as opposed to another thing. This is most widely known as an if/else
statement. In R, the if/else syntax is:
if (condition = TRUE) {
some R code
} else {
some other R code
}

For example,
x = 1
y = 3
36 CHAPTER 1. INTRODUCTION TO R

if (x > y) { # test if x > y

# if TRUE
z = x * y
print("x is larger than y")
} else {
# if FALSE
z = x + 5 * y
print("x is less than or equal to y")
}

#OUT> [1] "x is less than or equal to y"

z

#OUT> [1] 16

1.11.3 Loops

Loops are a very important programming construct. As the name suggests, in

a loop, the programming repeatedly loops over a set of instructions, until some
condition tells it to stop. A very powerful, yet simple, construction is that
the program can count how many steps it has done already - which may be
important to know for many algorithms. The syntax of a for loop (there are
others), is
for (ix in 1:10){ # does not have to be 1:10!
# loop body: gets executed each time
# the value of ix changes with each iteration
}

For example, consider this simple for loop, which will simply print the value of
the iterator (called i in this case) to screen:
for (i in 1:5){
print(i)
}

#OUT> [1] 1
#OUT> [1] 2
#OUT> [1] 3
#OUT> [1] 4
#OUT> [1] 5
Notice that instead of 1:5, we could have any kind of iterable collection:
for (i in c("mangos","bananas","apples")){
print(paste("I love",i)) # the paste function pastes together strings
}
1.11. PROGRAMMING BASICS 37

#OUT> [1] "I love mangos"

#OUT> [1] "I love bananas"
#OUT> [1] "I love apples"
We often also see nested loops, which are just what its name suggests:
for (i in 2:3){
# first nest: for each i
for (j in c("mangos","bananas","apples")){
# second nest: for each j
print(paste("Can I get",i,j,"please?"))
}
}

#OUT> [1] "Can I get 2 mangos please?"

#OUT> [1] "Can I get 2 bananas please?"
#OUT> [1] "Can I get 2 apples please?"
#OUT> [1] "Can I get 3 mangos please?"
#OUT> [1] "Can I get 3 bananas please?"
#OUT> [1] "Can I get 3 apples please?"
The important thing to note here is that you can do calculations with the
iterators while inside a loop.

1.11.4 Functions

So far we have been using functions, but haven’t actually discussed some of
their details. A function is a set of instructions that R executes for us, much like
those collected in a script file. The good thing is that functions are much more
flexible than scripts, since they can depend on input arguments, which change
the way the function behaves. Here is how to define a function:
function_name <- function(arg1,arg2=default_value){
# function body
# you do stuff with arg1 and arg2
# you can have any number of arguments, with or without defaults
# any valid `R` commands can be included here
# the last line is returned
}

And here is a trivial example of a function definition:

hello <- function(your_name = "Lord Vader"){
paste("You R most welcome,",your_name)
# we could also write:
# return(paste("You R most welcome,",your_name))
}
38 CHAPTER 1. INTRODUCTION TO R

# we call the function by typing it's name with round brackets

hello()

#OUT> [1] "You R most welcome, Lord Vader"

You see that by not specifying the argument your_name, R reverts to the default
value given. Try with your own name now!
Just typing the function name returns the actual definition to us, which is handy
sometimes:
hello

#OUT> function(your_name = "Lord Vader"){

#OUT> paste("You R most welcome,",your_name)
#OUT> # we could also write:
#OUT> # return(paste("You R most welcome,",your_name))
#OUT> }
#OUT> <environment: 0x7fc5d141ee78>
It’s instructive to consider that before we defined the function hello above, R
did not know what to do, had you called hello(). The function did not exist!
In this sense, we taught R a new trick. This feature to create new capabilities on
top of a core language is one of the most powerful characteristics of programming
languages. In general, it is good practice to split your code into several smaller
functions, rather than one long script file. It makes your code more readable,
and it is easier to track down mistakes.

[Link] Task 6

1. Write a for loop that counts down from 10 to 1, printing the value of the
iterator to the screen.
2. Modify that loop to write “i iterations to go” where i is the iterator
3. Modify that loop so that each iteration takes roughly one second. You can
achieve that by adding the command [Link](1) below the line that
prints “i iterations to go”.
Chapter 2

Working With Data

In this chapter we will first learn some basic concepts that help summarizing
data. Then, we will tackle a real-world task and read, clean, and summarize
data from the web.

2.1 Summary Statistics

R has built in functions for a large number of summary statistics. For numeric
variables, we can summarize data by looking at their center and spread, for
example.
# for the mpg dataset, we load:
library(ggplot2)

Central Tendency

Suppose we want to know the mean and median of all the values stored in the
[Link] column mpg$cty:

Measure R Result
Mean mean(mpg$cty) 16.8589744
Median median(mpg$cty) 17

Spread

How do the values in that column vary? How far spread out are they?

39
40 CHAPTER 2. WORKING WITH DATA

Measure R Result
Variance var(mpg$cty) 18.1130736
Standard Deviation sd(mpg$cty) 4.2559457
IQR IQR(mpg$cty) 5
Minimum min(mpg$cty) 9
Maximum max(mpg$cty) 35
Range range(mpg$cty) 9, 35

Categorical

For categorical variables, counts and percentages can be used for summary.
table(mpg$drv)

#OUT>
#OUT> 4 f r
#OUT> 103 106 25
table(mpg$drv) / nrow(mpg)

#OUT>
#OUT> 4 f r
#OUT> 0.4401709 0.4529915 0.1068376

2.2 Plotting
Now that we have some data to work with, and we have learned about the
data at the most basic level, our next tasks will be to visualize it. Often, a
proper visualization can illuminate features of the data that can inform further
analysis.
We will look at four methods of visualizing data by using the basic plot facilities
built-in with R:
• Histograms
• Barplots
• Boxplots
• Scatterplots

2.2.1 Histograms

When visualizing a single numerical variable, a histogram is useful. It summa-

rizes the distribution of values in a vector. In R you create one using the hist()
function:
2.2. PLOTTING 41

hist(mpg$cty)

Histogram of mpg$cty
80
60
Frequency

40
20
0

5 10 15 20 25 30 35

mpg$cty

The histogram function has a number of parameters which can be changed to

make our plot look much nicer. Use the ? operator to read the documentation
for the hist() to see a full list of these parameters.
hist(mpg$cty,
xlab = "Miles Per Gallon (City)",
main = "Histogram of MPG (City)", # main title
breaks = 12, # how many breaks?
col = "red",
border = "blue")
42 CHAPTER 2. WORKING WITH DATA

Histogram of MPG (City)

40
30
Frequency

20
10
0

10 15 20 25 30 35

Miles Per Gallon (City)

Importantly, you should always be sure to label your axes and give the plot a
title. The argument breaks is specific to hist(). Entering an integer will give
a suggestion to R for how many bars to use for the histogram. By default R will
attempt to intelligently guess a good number of breaks, but as we can see here,
it is sometimes useful to modify this yourself.

2.2.2 Barplots

Somewhat similar to a histogram, a barplot can provide a visual summary of a

categorical variable, or a numeric variable with a finite number of values, like a
ranking from 1 to 10.
barplot(table(mpg$drv))
2.2. PLOTTING 43
100
80
60
40
20
0

4 f r
barplot(table(mpg$drv),
xlab = "Drivetrain (f = FWD, r = RWD, 4 = 4WD)",
ylab = "Frequency",
main = "Drivetrains",
col = "dodgerblue",
border = "darkorange")

Drivetrains
100
80
Frequency

60
40
20
0

4 f r

Drivetrain (f = FWD, r = RWD, 4 = 4WD)

44 CHAPTER 2. WORKING WITH DATA

2.2.3 Boxplots

To visualize the relationship between a numerical and categorical variable, once

could use a boxplot. In the mpg dataset, the drv variable takes a small, finite
number of values. A car can only be front wheel drive, 4 wheel drive, or rear
wheel drive.
unique(mpg$drv)

#OUT> [1] "f" "4" "r"

First note that we can use a single boxplot as an alternative to a histogram for
visualizing a single numerical variable. To do so in R, we use the boxplot()
function. The box shows the interquartile range, the solid line in the middle is
the value of the median, the wiskers show 1.5 times the interquartile range, and
the dots are outliers.
boxplot(mpg$hwy)
45
40
35
30
25
20
15

However, more often we will use boxplots to compare a numerical variable for
different values of a categorical variable.
boxplot(hwy ~ drv, data = mpg)
2.2. PLOTTING 45

45
40
35
30
hwy

25
20
15

4 f r

drv

Here used the boxplot() command to create side-by-side boxplots. However,

since we are now dealing with two variables, the syntax has changed. The R
syntax hwy ~ drv, data = mpg reads “Plot the hwy variable against the drv
variable using the dataset mpg.” We see the use of a ~ (which specifies a formula)
and also a data = argument. This will be a syntax that is common to many
functions we will use in this course.
boxplot(hwy ~ drv, data = mpg,
xlab = "Drivetrain (f = FWD, r = RWD, 4 = 4WD)",
ylab = "Miles Per Gallon (Highway)",
main = "MPG (Highway) vs Drivetrain",
pch = 20,
cex = 2,
col = "darkorange",
border = "dodgerblue")
46 CHAPTER 2. WORKING WITH DATA

MPG (Highway) vs Drivetrain

45
40
Miles Per Gallon (Highway)

35
30
25
20
15

4 f r

Drivetrain (f = FWD, r = RWD, 4 = 4WD)

Again, boxplot() has a number of additional arguments which have the ability
to make our plot more visually appealing.

2.2.4 Scatterplots

Lastly, to visualize the relationship between two numeric variables we will use a
scatterplot. This can be done with the plot() function and the ~ syntax we
just used with a boxplot. (The function plot() can also be used more generally;
see the documentation for details.)
plot(hwy ~ displ, data = mpg)
2.2. PLOTTING 47

45
40
35
30
hwy

25
20
15

2 3 4 5 6 7

displ
plot(hwy ~ displ, data = mpg,
xlab = "Engine Displacement (in Liters)",
ylab = "Miles Per Gallon (Highway)",
main = "MPG (Highway) vs Engine Displacement",
pch = 20,
cex = 2,
col = "dodgerblue")

MPG (Highway) vs Engine Displacement

45
40
Miles Per Gallon (Highway)

35
30
25
20
15

2 3 4 5 6 7

Engine Displacement (in Liters)

48 CHAPTER 2. WORKING WITH DATA

2.2.5 ggplot

All of the above plots could also have been generated using the ggplot function
from the already loaded ggplot2 package. Which function you use is up to you,
but sometimes a plot is easier to build in base R (like in the boxplot example
maybe), sometimes the other way around.
ggplot(data = mpg,mapping = aes(x=displ,y=hwy)) + geom_point()

40

30
hwy

20

2 3 4 5 6 7
displ

ggplot is impossible to describe in brief terms, so please look at the package’s

website which provides excellent guidance. We will from time to time use ggplot
in this book, so you could familiarize yourself with it. Let’s quickly demonstrate
how one could further customize that first plot:
ggplot(data = mpg, mapping = aes(x=displ,y=hwy)) + # ggplot() makes base plot
geom_point(color="blue",size=2) + # how to show x and y?
scale_y_continuous(name="Miles Per Gallon (Highway)") + # name of y axis
scale_x_continuous(name="Engine Displacement (in Liters)") + # x axis
theme_bw() + # change the background
ggtitle("MPG (Highway) vs Engine Displacement") # add a title
2.3. SUMMARIZING TWO VARIABLES 49

MPG (Highway) vs Engine Displacement

40
Miles Per Gallon (Highway)

30

20

2 3 4 5 6 7
Engine Displacement (in Liters)

If you want to see ggplot in action, you could start with this and then look at
that very nice tutorial? It’s fun!

2.3 Summarizing Two Variables

We often are interested in how two variables are related to each other. The core
concepts here are covariance and correlation. Let’s generate some data on x
and y and plot them against each other:
Taking as example the data in this plot, the concepts covariance and correlation
relate to the following type of question:

Given we observe value of something like x = 2, say, can we expect a

high or a low value of y, on average? Something like y = 2 or rather
something like y = −2?

The answer to this type of question can be addressed by computing the covari-
ance of both variables:
cov(x,y)

#OUT> [1] 1.041195

Here, this gives a positive number, 1.04, indicating that as one variable lies above
it’s average, the other one does as well. In other words, it indicates a positive
50 CHAPTER 2. WORKING WITH DATA

3
2
1
y

0
−1
−2

−2 −1 0 1 2

x
Figure 2.1: How are x and y related?

relationship. What is less clear, however, how to interpret the magnitude of

1.04. Is that a strong or a weak positive association?
In fact, we cannot tell. This is because the covariance is measured in the same
units as the data, and those units often differ between both variables. There is
a better measure available to us though, the correlation, which is obtained by
standardizing each variable. By standardizing a variable x one means to divide
x by its standard deviation σx :

x
z=
σx

The correlation coefficient between x and y, commonly denoted rx,y , is then

defined as

cov(x, y)
rx,y = ,
σx σy

and we get rid of the units problem. In R, you can call directly
cor(x,y)

#OUT> [1] 0.9142495

2.3. SUMMARIZING TWO VARIABLES 51

−2.20 3.00 10.45

Figure 2.2: visual estimation on σ. The x-axis labels min and max as well as
mean of x.

Now this is better. Given that the correlation has to lie in [−1, 1], a value of
0.91 is indicative of a rather strong positive relationship for the data in figure
2.1

Note that x, y being drawn from a continuous distribution (they are joint nor-
mally distributed) had no implication for covariance and correlation: We can
compute those measures also for discrete random variables (like the throws of
two dice, as you will see in one of our tutorials).

2.3.1 Visually estimating σ

Sometimes it is useful to estimate the standard deviation of some data without

the help of a computer (for example during an exam ;-) ). If x is approximately
normally distributed, 95% of its observations will lie within a range of x̄± two
standard deviations of x. That is to say, four standard deviations of x cover
95% of its observations. Hence, a simple way to estimate the standard deviation
for a variable is to look at the range of x, and simply divide that number by
four.

This is illustrated in figure 2.2. Here we see that range(x)/4 gives 3.16 which
compares favourably to the actual standard deviation 3.
52 CHAPTER 2. WORKING WITH DATA

2.4 The tidyverse

Hadley Wickham is the author of R packages ggplot2 and also of dplyr (and
also a myriad of others). With ggplot2 he introduced what is called the gram-
mar of graphics (hence, gg) to R. Grammar in the sense that there are nouns
and verbs and a syntax, i.e. rules of how nouns and verbs are to be put to-
gether to construct an understandable sentence. He has extended the grammar
idea into various other packages. The tidyverse package is a collection of those
packages.
tidy data is data where:
• Each variable is a column
• Each observation is a row
• Each value is a cell
Fair enough, you might say, that is a regular spreadsheet. And you are right!
However, data comes to us not tidy most of the times, and we first need to clean,
or tidy, it up. Once it’s in tidy format, we can use the tools in the tidyverse
with great efficiency to analyse the data and stop worrying about which tool to
use.

2.4.1 Reading .csv data in the tidy way

We could have used the read_csv() function from the readr package to read
our example dataset from the previous chapter. The readr function read_csv()
has a number of advantages over the built-in [Link]. For example, it is much
faster reading larger data. It also uses the tibble package to read the data as a
tibble. A tibble is simply a data frame that prints with sanity. Notice
in the output below that we are given additional information such as dimension
and variable type.
library(readr) # you need `[Link]("readr")` once!
path = [Link](package="ScPoEconometrics","datasets","[Link]")
example_data_from_disk = read_csv(path)

2.4.2 Tidy [Link] are tibbles

Let’s grab some data from the ggplot2 package:

data(mpg,package = "ggplot2") # load dataset `mpg` from `ggplot2` package
head(mpg, n = 10)

#OUT> # A tibble: 10 x 11
#OUT> manufacturer model displ year cyl trans drv cty hwy fl class
#OUT> <chr> <chr> <dbl> <int> <int> <chr> <chr> <int> <int> <chr> <chr>
2.4. THE TIDYVERSE 53

#OUT> 1 audi a4 1.8 1999 4 auto~ f 18 29 p comp~

#OUT> 2 audi a4 1.8 1999 4 manu~ f 21 29 p comp~
#OUT> 3 audi a4 2 2008 4 manu~ f 20 31 p comp~
#OUT> 4 audi a4 2 2008 4 auto~ f 21 30 p comp~
#OUT> 5 audi a4 2.8 1999 6 auto~ f 16 26 p comp~
#OUT> 6 audi a4 2.8 1999 6 manu~ f 18 26 p comp~
#OUT> 7 audi a4 3.1 2008 6 auto~ f 18 27 p comp~
#OUT> 8 audi a4 q~ 1.8 1999 4 manu~ 4 18 26 p comp~
#OUT> 9 audi a4 q~ 1.8 1999 4 auto~ 4 16 25 p comp~
#OUT> 10 audi a4 q~ 2 2008 4 manu~ 4 20 28 p comp~

The function head() will display the first n observations of the data frame, as
we have seen. The head() function was more useful before tibbles. Notice
that mpg is a tibble already, so the output from head() indicates there are only
10 observations. Note that this applies to head(mpg, n = 10) and not mpg
itself. Also note that tibbles print a limited number of rows and columns by
default. The last line of the printed output indicates with rows and columns
were omitted.
mpg

#OUT> # A tibble: 234 x 11

#OUT> manufacturer model displ year cyl trans drv cty hwy fl class
#OUT> <chr> <chr> <dbl> <int> <int> <chr> <chr> <int> <int> <chr> <chr>
#OUT> 1 audi a4 1.8 1999 4 auto~ f 18 29 p comp~
#OUT> 2 audi a4 1.8 1999 4 manu~ f 21 29 p comp~
#OUT> 3 audi a4 2 2008 4 manu~ f 20 31 p comp~
#OUT> 4 audi a4 2 2008 4 auto~ f 21 30 p comp~
#OUT> 5 audi a4 2.8 1999 6 auto~ f 16 26 p comp~
#OUT> 6 audi a4 2.8 1999 6 manu~ f 18 26 p comp~
#OUT> 7 audi a4 3.1 2008 6 auto~ f 18 27 p comp~
#OUT> 8 audi a4 q~ 1.8 1999 4 manu~ 4 18 26 p comp~
#OUT> 9 audi a4 q~ 1.8 1999 4 auto~ 4 16 25 p comp~
#OUT> 10 audi a4 q~ 2 2008 4 manu~ 4 20 28 p comp~
#OUT> # ... with 224 more rows

Let’s look at str as well to get familiar with the content of the data:
str(mpg)

#OUT> Classes 'tbl_df', 'tbl' and '[Link]': 234 obs. of 11 variables:

#OUT> $ manufacturer: chr "audi" "audi" "audi" "audi" ...
#OUT> $ model : chr "a4" "a4" "a4" "a4" ...
#OUT> $ displ : num 1.8 1.8 2 2 2.8 2.8 3.1 1.8 1.8 2 ...
#OUT> $ year : int 1999 1999 2008 2008 1999 1999 2008 1999 1999 2008 ...
#OUT> $ cyl : int 4 4 4 4 6 6 6 4 4 4 ...
#OUT> $ trans : chr "auto(l5)" "manual(m5)" "manual(m6)" "auto(av)" ...
#OUT> $ drv : chr "f" "f" "f" "f" ...
54 CHAPTER 2. WORKING WITH DATA

#OUT> $ cty : int 18 21 20 21 16 18 18 18 16 20 ...

#OUT> $ hwy : int 29 29 31 30 26 26 27 26 25 28 ...
#OUT> $ fl : chr "p" "p" "p" "p" ...
#OUT> $ class : chr "compact" "compact" "compact" "compact" ...
In this dataset an observation is for a particular model-year of a car, and the
variables describe attributes of the car, for example its highway fuel efficiency.
To understand more about the data set, we use the ? operator to pull up the
documentation for the data.
?mpg

Working with tibbles is mostly the same as working with plain [Link]:
names(mpg)

#OUT> [1] "manufacturer" "model" "displ" "year"

#OUT> [5] "cyl" "trans" "drv" "cty"
#OUT> [9] "hwy" "fl" "class"
mpg$year

#OUT> [1] 1999 1999 2008 2008 1999 1999 2008 1999 1999 2008 2008 1999 1999 2008
#OUT> [15] 2008 1999 2008 2008 2008 2008 2008 1999 2008 1999 1999 2008 2008 2008
#OUT> [29] 2008 2008 1999 1999 1999 2008 1999 2008 2008 1999 1999 1999 1999 2008
#OUT> [43] 2008 2008 1999 1999 2008 2008 2008 2008 1999 1999 2008 2008 2008 1999
#OUT> [57] 1999 1999 2008 2008 2008 1999 2008 1999 2008 2008 2008 2008 2008 2008
#OUT> [71] 1999 1999 2008 1999 1999 1999 2008 1999 1999 1999 2008 2008 1999 1999
#OUT> [85] 1999 1999 1999 2008 1999 2008 1999 1999 2008 2008 1999 1999 2008 2008
#OUT> [99] 2008 1999 1999 1999 1999 1999 2008 2008 2008 2008 1999 1999 2008 2008
#OUT> [113] 1999 1999 2008 1999 1999 2008 2008 2008 2008 2008 2008 2008 1999 1999
#OUT> [127] 2008 2008 2008 2008 1999 2008 2008 1999 1999 1999 2008 1999 2008 2008
#OUT> [141] 1999 1999 1999 2008 2008 2008 2008 1999 1999 2008 1999 1999 2008 2008
#OUT> [155] 1999 1999 1999 2008 2008 1999 1999 2008 2008 2008 2008 1999 1999 1999
#OUT> [169] 1999 2008 2008 2008 2008 1999 1999 1999 1999 2008 2008 1999 1999 2008
#OUT> [183] 2008 1999 1999 2008 1999 1999 2008 2008 1999 1999 2008 1999 1999 1999
#OUT> [197] 2008 2008 1999 2008 1999 1999 2008 1999 1999 2008 2008 1999 1999 2008
#OUT> [211] 2008 1999 1999 1999 1999 2008 2008 2008 2008 1999 1999 1999 1999 1999
#OUT> [225] 1999 2008 2008 1999 1999 2008 2008 1999 1999 2008
mpg$hwy

#OUT> [1] 29 29 31 30 26 26 27 26 25 28 27 25 25 25 25 24 25 23 20 15 20 17 17
#OUT> [24] 26 23 26 25 24 19 14 15 17 27 30 26 29 26 24 24 22 22 24 24 17 22 21
#OUT> [47] 23 23 19 18 17 17 19 19 12 17 15 17 17 12 17 16 18 15 16 12 17 17 16
#OUT> [70] 12 15 16 17 15 17 17 18 17 19 17 19 19 17 17 17 16 16 17 15 17 26 25
#OUT> [93] 26 24 21 22 23 22 20 33 32 32 29 32 34 36 36 29 26 27 30 31 26 26 28
#OUT> [116] 26 29 28 27 24 24 24 22 19 20 17 12 19 18 14 15 18 18 15 17 16 18 17
2.4. THE TIDYVERSE 55

#OUT> [139] 19 19 17 29 27 31 32 27 26 26 25 25 17 17 20 18 26 26 27 28 25 25 24
#OUT> [162] 27 25 26 23 26 26 26 26 25 27 25 27 20 20 19 17 20 17 29 27 31 31 26
#OUT> [185] 26 28 27 29 31 31 26 26 27 30 33 35 37 35 15 18 20 20 22 17 19 18 20
#OUT> [208] 29 26 29 29 24 44 29 26 29 29 29 29 23 24 44 41 29 26 28 29 29 29 28
#OUT> [231] 29 26 26 26
Subsetting is also similar to dataframe. Here, we find fuel efficient vehicles
earning over 35 miles per gallon and only display manufacturer, model and
year.
# mpg[row condition, col condition]
mpg[mpg$hwy > 35, c("manufacturer", "model", "year")]

#OUT> # A tibble: 6 x 3
#OUT> manufacturer model year
#OUT> <chr> <chr> <int>
#OUT> 1 honda civic 2008
#OUT> 2 honda civic 2008
#OUT> 3 toyota corolla 2008
#OUT> 4 volkswagen jetta 1999
#OUT> 5 volkswagen new beetle 1999
#OUT> 6 volkswagen new beetle 1999
An alternative would be to use the subset() function, which has a much more
readable syntax.
subset(mpg, subset = hwy > 35, select = c("manufacturer", "model", "year"))

Lastly, and most tidy, we could use the filter and select functions from
the dplyr package which introduces the pipe operator f(x) %>% g(z) from the
magrittr package. This operator takes the output of the first command, for
example y = f(x), and passes it as the first argument to the next function,
i.e. we’d obtain g(y,z) here.1
library(dplyr)
mpg %>%
filter(hwy > 35) %>%
select(manufacturer, model, year)

#OUT> # A tibble: 6 x 3
#OUT> manufacturer model year
#OUT> <chr> <chr> <int>
#OUT> 1 honda civic 2008
#OUT> 2 honda civic 2008
#OUT> 3 toyota corolla 2008
1 A pipe is a concept from the Unix world, where it means to take the output of some

command, and pass it on to another command. This way, one can construct a pipeline of
commands. For additional info on the pipe operator in R, you might be interested in this
tutorial.
56 CHAPTER 2. WORKING WITH DATA

#OUT> 4 volkswagen jetta 1999

#OUT> 5 volkswagen new beetle 1999
#OUT> 6 volkswagen new beetle 1999

Note that the above syntax is equivalent to the following pipe-free command
(which is much harder to read!):
library(dplyr)
select(filter(mpg, hwy > 35), manufacturer, model, year)

#OUT> # A tibble: 6 x 3
#OUT> manufacturer model year
#OUT> <chr> <chr> <int>
#OUT> 1 honda civic 2008
#OUT> 2 honda civic 2008
#OUT> 3 toyota corolla 2008
#OUT> 4 volkswagen jetta 1999
#OUT> 5 volkswagen new beetle 1999
#OUT> 6 volkswagen new beetle 1999

All three approaches produce the same results. Which you use will be largely
based on a given situation as well as your preference.

[Link] Task 1

1. Make sure to have the mpg dataset loaded by typing data(mpg) (and
library(ggplot2) if you haven’t!). Use the table function to find out
how many cars were built by mercury?
2. What is the average year the audi’s were built in this dataset? Use the
function mean on the subset of column year that corresponds to audi. (Be
careful: subsetting a tibble returns a tibble (and not a vector)!. so get
the year column after you have subset the tibble.)
3. Use the dplyr piping syntax from above first with group_by and then
with summarise(newvar=your_expression) to find the mean year by all
manufacturers (i.e. same as previous task, but for all manufacturers. don’t
write a loop!).

2.4.3 Tidy Example: Importing Non-Tidy Excel Data

The data we will look at is from Eurostat on demography and migration. You
should download the data yourself (click on previous link, then drill down to
database by themes > Population and social conditions > Demograph and mi-
gration > Population change - Demographic balance and crude rates at national
level (demo_gind)).
2.4. THE TIDYVERSE 57

Once downloaded, we can read the data with the function read_excel from the
package readxl, again part of the tidyverse suite.
It’s important to know how the data is organized in the spreadsheet. Open the
file with Excel to see:
• There is a heading which we don’t need.
• There are 5 rows with info that we don’t need.
• There is one table per variable (total population, males, females, etc)
• Each table has one row for each country, and one column for each year.
• As such, this data is not tidy.
Now we will read the first chunk of data, from the first table: total population:
library(readxl) # load the library
# Notice that if you installed the R package of this book,
# you have the .xls data file already at
# `[Link](package="ScPoEconometrics",
# "datasets","demo_gind.xls")`
# otherwise:
# * download the file to your computer
# * change the argument `path` to where you downloaded it
# you may want to change your working directory with `setwd("your/directory")
# or in RStudio by clicking Session > Set Working Directory

# total population in raw format

tot_pop_raw = read_excel(
path = [Link](package="ScPoEconometrics",
"datasets","demo_gind.xls"),
sheet="Data", # which sheet
range="A9:K68") # which excel cell range to read
names(tot_pop_raw)[1] <- "Country" # lets rename the first column
tot_pop_raw

#OUT> # A tibble: 59 x 11
#OUT> Country `2008` `2009` `2010` `2011` `2012` `2013` `2014` `2015` `2016`
#OUT> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
#OUT> 1 Europe~ 50029~ 50209~ 50317~ 50296~ 50404~ 50516~ 50701~ 50854~ 51027~
#OUT> 2 Europe~ 43872~ 44004~ 44066~ 43994~ 44055~ 44125~ 44266~ 44366~ 44489~
#OUT> 3 Europe~ 49598~ 49778~ 49886~ 49867~ 49977~ 50090~ 50276~ 50431~ 50608~
#OUT> 4 Euro a~ 33309~ 33447~ 33526~ 33457~ 33528~ 33604~ 33754~ 33856~ 33988~
#OUT> 5 Euro a~ 32988~ 33128~ 33212~ 33152~ 33228~ 33307~ 33459~ 33563~ 33699~
#OUT> 6 Belgium 10666~ 10753~ 10839~ 11000~ 11075~ 11137~ 11180~ 11237~ 11311~
#OUT> 7 Bulgar~ 75180~ 74671~ 74217~ 73694~ 73272~ 72845~ 72456~ 72021~ 71537~
#OUT> 8 Czech ~ 10343~ 10425~ 10462~ 10486~ 10505~ 10516~ 10512~ 10538~ 10553~
#OUT> 9 Denmark 54757~ 55114~ 55347~ 55606~ 55805~ 56026~ 56272~ 56597~ 57072~
#OUT> 10 German~ 82217~ 82002~ 81802~ 80222~ 80327~ 80523~ 80767~ 81197~ 82175~
#OUT> # ... with 49 more rows, and 1 more variable: `2017` <chr>
58 CHAPTER 2. WORKING WITH DATA

This shows a tibble, which we encountered just above. The column names are
Country,2008,2009,..., and the rows are numbered 1,2,3,.... Notice, in
particular, that all columns seem to be of type <chr>, i.e. characters - a string,
not a number! We’ll have to fix that, as this is clearly numeric data.

[Link] tidyr

In the previous tibble, each year is a column name (like 2008) instead of all
years being collected in one column year. We really would like to have several
rows for each Country, one row per year. We want to gather() all years into a
new column to tidy this up - and here is how:
1. specify which columns are to be gathered: in our case, all years (note
that paste(2008:2017) produces a vector like ["2008", "2009",
"2010",...])
2. say what those columns should be gathered into, i.e. what is the key for
those values: we’ll call it year.
3. Finally, what is the name of the new resulting column, containing the
value from each cell: let’s call it counts.
library(tidyr) # for the gather function
tot_pop = gather(tot_pop_raw, paste(2008:2017),key="year", value = "counts")
tot_pop

#OUT> # A tibble: 590 x 3

#OUT> Country year counts
#OUT> <chr> <chr> <chr>
#OUT> 1 European Union (current composition) 2008 500297033
#OUT> 2 European Union (without United Kingdom) 2008 438725386
#OUT> 3 European Union (before the accession of Croatia) 2008 495985066
#OUT> 4 Euro area (19 countries) 2008 333096775
#OUT> 5 Euro area (18 countries) 2008 329884170
#OUT> 6 Belgium 2008 10666866
#OUT> 7 Bulgaria 2008 7518002
#OUT> 8 Czech Republic 2008 10343422
#OUT> 9 Denmark 2008 5475791
#OUT> 10 Germany (until 1990 former territory of the FRG) 2008 82217837
#OUT> # ... with 580 more rows
That’s better! However, counts is still chr! Let’s convert it to a number:
tot_pop$counts = [Link](tot_pop$counts)

#OUT> Warning: NAs introduced by coercion

tot_pop

#OUT> # A tibble: 590 x 3

2.4. THE TIDYVERSE 59

#OUT> Country year counts

#OUT> <chr> <chr> <int>
#OUT> 1 European Union (current composition) 2008 500297033
#OUT> 2 European Union (without United Kingdom) 2008 438725386
#OUT> 3 European Union (before the accession of Croatia) 2008 495985066
#OUT> 4 Euro area (19 countries) 2008 333096775
#OUT> 5 Euro area (18 countries) 2008 329884170
#OUT> 6 Belgium 2008 10666866
#OUT> 7 Bulgaria 2008 7518002
#OUT> 8 Czech Republic 2008 10343422
#OUT> 9 Denmark 2008 5475791
#OUT> 10 Germany (until 1990 former territory of the FRG) 2008 82217837
#OUT> # ... with 580 more rows
Now you can see that column counts is indeed int, i.e. an integer number,
and we are fine. The Warning: NAs introduced by coercion means that R
converted some values to NA, because it couldn’t convert them into numeric.
More below!

[Link] dplyr

The transform chapter of Hadley Wickham’s book is a great place

to read up more on using dplyr.
With dplyr you can do the following operations on [Link] and tibbles:
• Choose observations based on a certain value (i.e. subset): filter()
• Reorder rows: arrange()
• Select variables by name: select()
• Create new variables out of existing ones: mutate()
• Summarise variables: summarise()
All of those verbs can be used with group_by(), where we apply the respective
operation on a group of the dataframe/tibble. For example, on our tot_pop
tibble we will now
• filter
• mutate
• and plot the resulting values
Let’s get a plot of the populations of France, the UK and Italy over time, in
terms of millions of people. We will make use of the piping syntax of dplyr
which we introduced just above.
library(dplyr) # for %>%, filter, mutate, ...
# 1. take the [Link] `tot_pop`
tot_pop %>%
# 2. pipe it into the filter function
60 CHAPTER 2. WORKING WITH DATA

# filter on Country being one of "France","United Kingdom" or "Italy"

filter(Country %in% c("France","United Kingdom","Italy")) %>%
# 3. pipe the result into the mutate function
# create a new column called millions
mutate(millions = counts / 1e6) %>%
# 4. pipe the result into ggplot to make a plot
ggplot(mapping = aes(x=year,y=millions,color=Country,group=Country)) + geom_line(size

65.0

Country
millions

France

62.5 Italy
United Kingdom

60.0

2008 2009 2010 2011 2012 2013 2014 2015 2016 2017
year

Arrange a tibble

• What are the top/bottom 5 most populated areas?

top5 = tot_pop %>%
arrange(desc(counts)) %>% # arrange in descending order of col `counts`
top_n(5)

bottom5 = tot_pop %>%

arrange(desc(counts)) %>%
top_n(-5)
# let's see top 5
top5

#OUT> # A tibble: 5 x 3
#OUT> Country year counts
2.4. THE TIDYVERSE 61

#OUT> <chr> <chr> <int>

#OUT> 1 European Economic Area (EU28 - current composition, plus~ 2017 5.17e8
#OUT> 2 European Economic Area (EU28 - current composition, plus~ 2016 5.16e8
#OUT> 3 European Economic Area (EU28 - current composition, plus~ 2015 5.14e8
#OUT> 4 European Economic Area (EU27 - before the accession of C~ 2017 5.13e8
#OUT> 5 European Economic Area (EU28 - current composition, plus~ 2014 5.12e8
# and bottom 5
bottom5

#OUT> # A tibble: 5 x 3
#OUT> Country year counts
#OUT> <chr> <chr> <int>
#OUT> 1 San Marino 2015 32789
#OUT> 2 San Marino 2014 32520
#OUT> 3 San Marino 2008 32054
#OUT> 4 San Marino 2011 31863
#OUT> 5 San Marino 2009 31269
Now this is not exactly what we wanted. It’s always the same country in both
top and bottom, because there are multiple years per country. Let’s compute
average population over the last 5 years and rank according to that:
topbottom = tot_pop %>%
group_by(Country) %>%
filter(year > 2012) %>%
summarise(mean_count = mean(counts)) %>%
arrange(desc(mean_count))

top5 = topbottom %>% top_n(5)

bottom5 = topbottom %>% top_n(-5)
top5

#OUT> # A tibble: 5 x 2
#OUT> Country mean_count
#OUT> <chr> <dbl>
#OUT> 1 European Economic Area (EU28 - current composition, plus IS,~ 514029320
#OUT> 2 European Economic Area (EU27 - before the accession of Croat~ 509813491.
#OUT> 3 European Union (current composition) 508502858.
#OUT> 4 European Union (before the accession of Croatia) 504287028.
#OUT> 5 European Union (without United Kingdom) 443638309.
bottom5

#OUT> # A tibble: 5 x 2
#OUT> Country mean_count
#OUT> <chr> <dbl>
#OUT> 1 Luxembourg 563319.
#OUT> 2 Malta 440467.
62 CHAPTER 2. WORKING WITH DATA

#OUT> 3 Iceland 329501.

#OUT> 4 Liechtenstein 37353
#OUT> 5 San Marino 33014.
That’s better!

Look for NAs in a tibble

Sometimes data is missing, and R represents it with the special value NA (not
available). It is good to know where in our dataset we are going to encounter
any missing values, so the task here is: let’s produce a table that has three
columns:
1. the names of countries with missing data
2. how many years of data are missing for each of those
3. and the actual years that are missing
missings = tot_pop %>%
filter([Link](counts)) %>% # [Link](x) returns TRUE if x is NA
group_by(Country) %>%
summarise(n_missing = n(),years = paste(year,collapse = ", "))
knitr:::kable(missings) # knitr:::kable makes a nice table

Country n_missing years

Albania 2 2010, 2012
Andorra 2 2014, 2015
Armenia 1 2014
France (metropolitan) 4 2014, 2015, 2016, 2017
Georgia 1 2013
Monaco 7 2008, 2009, 2010, 2011, 2012, 2013, 2014
Russia 4 2013, 2015, 2016, 2017
San Marino 1 2010

Males and Females

Let’s look at the numbers by male and female population. They are in the
same xls file, but at different cell ranges. Also, I just realised that the special
character : indicates missing data. We can feed that to read_excel and that
will spare us the need to convert data types afterwards. Let’s see:
females_raw = read_excel(
path = [Link](package="ScPoEconometrics",
"datasets","demo_gind.xls"),
sheet="Data", # which sheet
range="A141:K200", # which excel cell range to read
na=":" ) # missing data indicator
2.4. THE TIDYVERSE 63

names(females_raw)[1] <- "Country" # lets rename the first column

females_raw

#OUT> # A tibble: 59 x 11
#OUT> Country `2008` `2009` `2010` `2011` `2012` `2013` `2014` `2015` `2016`
#OUT> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#OUT> 1 Europe~ 2.56e8 2.57e8 2.58e8 2.58e8 2.58e8 2.59e8 2.60e8 2.60e8 2.61e8
#OUT> 2 Europe~ 2.25e8 2.26e8 2.26e8 2.26e8 2.26e8 2.26e8 2.27e8 2.27e8 2.28e8
#OUT> 3 Europe~ 2.54e8 2.55e8 2.55e8 2.56e8 2.56e8 2.57e8 2.57e8 2.58e8 2.59e8
#OUT> 4 Euro a~ 1.71e8 1.71e8 1.72e8 1.72e8 1.72e8 1.72e8 1.73e8 1.73e8 1.74e8
#OUT> 5 Euro a~ 1.69e8 1.70e8 1.70e8 1.70e8 1.70e8 1.71e8 1.71e8 1.72e8 1.72e8
#OUT> 6 Belgium 5.44e6 5.48e6 5.53e6 5.60e6 5.64e6 5.67e6 5.69e6 5.71e6 5.74e6
#OUT> 7 Bulgar~ 3.86e6 3.83e6 3.81e6 3.78e6 3.76e6 3.74e6 3.72e6 3.70e6 3.68e6
#OUT> 8 Czech ~ 5.28e6 5.31e6 5.33e6 5.34e6 5.35e6 5.35e6 5.35e6 5.36e6 5.37e6
#OUT> 9 Denmark 2.76e6 2.78e6 2.79e6 2.80e6 2.81e6 2.82e6 2.83e6 2.85e6 2.87e6
#OUT> 10 German~ 4.19e7 4.18e7 4.17e7 4.11e7 4.11e7 4.11e7 4.12e7 4.14e7 4.17e7
#OUT> # ... with 49 more rows, and 1 more variable: `2017` <dbl>
You can see that R now correctly read the numbers as such, after we told it that
the : character has the special missing meaning: before, it coerced the entire
2008 column (for example) to be of type chr after it hit the first :. We had
to manually convert the column back to numeric, in the process automatically
coercing the :s into NA. Now we addressed that issue directly. Let’s also get the
male data in the same way:
males_raw = read_excel(
path = [Link](package="ScPoEconometrics",
"datasets","demo_gind.xls"),
sheet="Data", # which sheet
range="A75:K134", # which excel cell range to read
na=":" ) # missing data indicator
names(males_raw)[1] <- "Country" # lets rename the first column

Next step was to tidy up this data, just as before:

females = gather(females_raw, paste(2008:2017),key="year", value = "counts")
males = gather(males_raw, paste(2008:2017),key="year", value = "counts")

Let’s try to tweak our above plot to show the same data in two separate panels:
one for males and one for females. This is easiest to do with ggplot if we have
all the data in one single [Link] (or tibble), and marked with a group
identifier. Let’s first add this to both datasets, and then let’s just combine both
into one:
females$sex = "female"
males$sex = "male"
sexes = rbind(males,females) # "row bind" 2 [Link]
sexes
64 CHAPTER 2. WORKING WITH DATA

#OUT> # A tibble: 1,180 x 4

#OUT> Country year counts sex
#OUT> <chr> <chr> <dbl> <chr>
#OUT> 1 European Union (current composition) 2008 243990548 male
#OUT> 2 European Union (without United Kingdom) 2008 213826199 male
#OUT> 3 European Union (before the accession of Croatia) 2008 241913560 male
#OUT> 4 Euro area (19 countries) 2008 162516883 male
#OUT> 5 Euro area (18 countries) 2008 161029464 male
#OUT> 6 Belgium 2008 5224309 male
#OUT> 7 Bulgaria 2008 3660367 male
#OUT> 8 Czech Republic 2008 5065117 male
#OUT> 9 Denmark 2008 2712666 male
#OUT> 10 Germany (until 1990 former territory of the FRG) 2008 40274292 male
#OUT> # ... with 1,170 more rows
Now that we have all the data nice and tidy in a [Link], this is a very
small change to our previous plotting code:
sexes %>%
filter(Country %in% c("France","United Kingdom","Italy")) %>%
mutate(millions = counts / 1e6) %>%
ggplot(mapping = aes(x=[Link](year,format="%Y"), # convert to `Date`
y=millions,colour=Country,group=Country)) +
geom_line() +
scale_x_date(name = "year") + # rename x axis
facet_wrap(~sex) # make two panels, splitting by groups `sex`

female male

34

Country
32
millions

France
Italy
United Kingdom

30

2010 2012 2014 2016 2018 2010 2012 2014 2016 2018
year
2.4. THE TIDYVERSE 65

Always Compare to Germany :-)

How do our three countries compare with respect to the biggest country in
the EU in terms of population? What fraction of Germany does the French
population make in any given year, for example?
# remember that the pipe operator %>% takes the
# result of the previous operation and passes it
# as the *first* argument to the next function call
merge_GER <- tot_pop %>%
# 1. subset to countries of interest
filter(
Country %in%
c("France",
"United Kingdom",
"Italy")
) %>%
# 2. group data by year
group_by(year) %>%
# 3. add GER's count as new column *by year*
left_join(
# Germany only
filter(tot_pop,
Country %in% "Germany including former GDR"),
# join back in `by year`
by="year")
merge_GER

#OUT> # A tibble: 30 x 5
#OUT> # Groups: year [10]
#OUT> Country.x year counts.x Country.y counts.y
#OUT> <chr> <chr> <int> <chr> <int>
#OUT> 1 France 2008 64007193 Germany including former GDR 82217837
#OUT> 2 Italy 2008 58652875 Germany including former GDR 82217837
#OUT> 3 United Kingdom 2008 61571647 Germany including former GDR 82217837
#OUT> 4 France 2009 64350226 Germany including former GDR 82002356
#OUT> 5 Italy 2009 59000586 Germany including former GDR 82002356
#OUT> 6 United Kingdom 2009 62042343 Germany including former GDR 82002356
#OUT> 7 France 2010 64658856 Germany including former GDR 81802257
#OUT> 8 Italy 2010 59190143 Germany including former GDR 81802257
#OUT> 9 United Kingdom 2010 62510197 Germany including former GDR 81802257
#OUT> 10 France 2011 64978721 Germany including former GDR 80222065
#OUT> # ... with 20 more rows
Here you see that the merge (or join) operation labelled col.x and col.y if
both datasets contained a column called col. Now let’s continue to compute
what proportion of german population each country amounts to:
66 CHAPTER 2. WORKING WITH DATA

names(merge_GER)[1] <- "Country"

merge_GER %>%
mutate(prop_GER = 100 * counts.x / counts.y) %>%
# 5. plot
ggplot(mapping =
aes(x = year,
y = prop_GER,
color = Country,
group = Country)) +
geom_line(size=1) +
scale_y_continuous("percent of German population") +
theme_bw() # new theme for a change?

80.0
percent of German population

77.5
Country
France
Italy
United Kingdom
75.0

72.5

2008 2009 2010 2011 2012 2013 2014 2015 2016 2017
year
Chapter 3

Linear Regression

In this chapter we will learn an additional way how one can represent the rela-
tionship between outcome, or dependent variable variable y and an explanatory
or independent variable x. We will refer throughout to the graphical represen-
tation of a collection of independent observations on x and y, i.e., a dataset.

3.1 How are x and y related?

3.1.1 Data on Cars

We will look at the built-in cars dataset. Let’s get a view of this by just typing
View(cars) in Rstudio. You can see something like this:
#OUT> speed dist
#OUT> 1 4 2
#OUT> 2 4 10
#OUT> 3 7 4
#OUT> 4 7 22
#OUT> 5 8 16
#OUT> 6 9 10
We have a [Link] with two columns: speed and dist. Type help(cars)
to find out more about the dataset. There you could read that
The data give the speed of cars (mph) and the distances taken to
stop (ft).
It’s good practice to know the extent of a dataset. You could just type
dim(cars)

67
68 CHAPTER 3. LINEAR REGRESSION

#OUT> [1] 50 2
to find out that we have 50 rows and 2 columns. A central question that we
want to ask now is the following:

3.1.2 How are speed and dist related?

The simplest way to start is to plot the data. Remembering that we view each
row of a [Link] as an observation, we could just label one axis of a graph
speed, and the other one dist, and go through our table above row by row. We
just have to read off the x/y coordinates and mark them in the graph. In R:
plot(dist ~ speed, data = cars,
xlab = "Speed (in Miles Per Hour)",
ylab = "Stopping Distance (in Feet)",
main = "Stopping Distance vs Speed",
pch = 20,
cex = 2,
col = "red")

Stopping Distance vs Speed

100 120
Stopping Distance (in Feet)

80
60
40
20
0

5 10 15 20 25

Speed (in Miles Per Hour)

Here, each dot represents one observation. In this case, one particular measure-
ment speed and dist for a car. Now, again:

How are speed and dist related? How could one best summarize this
relationship?
3.1. HOW ARE X AND Y RELATED? 69

One thing we could do, is draw a straight line through this scatterplot, like so:
plot(dist ~ speed, data = cars,
xlab = "Speed (in Miles Per Hour)",
ylab = "Stopping Distance (in Feet)",
main = "Stopping Distance vs Speed",
pch = 20,
cex = 2,
col = "red")
abline(a = 60,b = 0,lw=3)

Stopping Distance vs Speed

100 120
Stopping Distance (in Feet)

80
60
40
20
0

5 10 15 20 25

Speed (in Miles Per Hour)

Now that doesn’t seem a particularly good way to summarize the relationship.
Clearly, a better line would be not be flat, but have a slope, i.e. go upwards:
70 CHAPTER 3. LINEAR REGRESSION

Stopping Distance vs Speed

100 120
Stopping Distance (in Feet)

80
60
40
20
0

5 10 15 20 25

Speed (in Miles Per Hour)

That is slightly better. However, the line seems at too high a level - the point
at which it crosses the y-axis is called the intercept; and it’s too high. We just
learned how to represent a line, i.e. with two numbers called intercept and slope.
Let’s write down a simple formula which represents a line where some outcome
z is related to a variable x:

z = b0 + b1 x (3.1)

Here b0 represents the value of the intercept (i.e. z when x = 0), and b1 is the
value of the slope. The question for us is now: How to choose the number b0
and b1 such that the result is the good line?

3.1.3 Choosing the Best Line

In order to be able to reason about good or bad line, we need to denote the
output of equation (3.1). We call the value ŷi the predicted value for obseration
i, after having chosen some particular values b0 and b1 :

ŷi = b0 + b1 xi (3.2)

In general it is likely that we won’t be able to choose b0 and b1 in such as way

as to provide a perfect prediction, i.e. one where ŷi = yi for all i. That is, we
expect to make an error in our prediction ŷi , so let’s denote this value ei . If
we acknowlegdge that we will make errors, let’s at least make them as small as
possible! Exactly this is going to be our task now.
3.1. HOW ARE X AND Y RELATED? 71

10
8
y

6
4
2

0 5 10

Figure 3.1: The best line and its errors

Suppose we have the following set of 9 observations on x and y, and we put the
best straight line into it, that we can think of. It would look like this:

Here, the red arrows indicate the distance between the prediction (i.e. the
black line) to each data point, in other words, each arrow is a particular ei . An
upward pointing arrow indicates a positive value of a particular ei , and vice
versa for downward pointing arrows. The erros are also called residuals, which
comes from the way can write the equation for this relationship between two
particular values (yi , xi ) belonging to observation i:

yi = b0 + b1 xi + ei (3.3)

You realize of course that ŷi = yi − ei , which just means that our prediction is
the observed value yi minus any error ei we make. In other words, ei is what is
left to be explained on top of the line b0 + b1 xi , hence, it’s a residual to explain
yi . Here are y, ŷ and the resulting e which are plotted in figure 3.1:
72 CHAPTER 3. LINEAR REGRESSION

x y y_hat error
0.00 2.09 2.57 -0.48
1.25 2.79 3.41 -0.62
2.50 6.49 4.25 2.24
3.75 1.71 5.10 -3.39
5.00 9.89 5.94 3.95
6.25 7.62 6.78 0.83
7.50 4.86 7.63 -2.77
8.75 7.38 8.47 -1.09
10.00 10.63 9.31 1.32
If our line was a perfect fit to the data, all ei = 0, and the column error
would display 0 for each row - there would be no errors at all. (All points in
figure 3.1 would perfectly line up on a straight line).
Now, back to our claim that this particular line is the best line. What exactly
characterizes this best line? We now come back to what we said above - how to
make the errors as small as possible? Keeping in mind that each residual ei is
yi − ŷi , we have the following minization problem to solve:

ei = yi − ŷi = yi − (b0 + b1 xi ) (3.4)

| {z }
prediction

∑
N
e21 + · · · + e2N = e2i ≡ SSR(b0 , b1 ) (3.5)
i=1
∑
N
2
(b0 , b1 ) = arg min [yi − (int + slope xi )] (3.6)
int,slope
i=1

The best line chooses b0 and b1 so as to minimize the sum of squared

residuals (SSR).

Wait a moment, why squared residuals? This is easy to understand: suppose

that instead, we wanted to just make the sum of the arrows in figure 3.1 as small
as possible (that is, no squares). Choosing our line to make this number small
would not give a particularly good representation of the data – given that errors
of opposite sign and equal magnitude offset, we could have very long arrows
(but of opposite signs), and a poor resulting line. Squaring each error avoids
this (because now negative errors get positive values!)
We illustrate this in figure 3.2. This is the same data as in figure 3.1, but instead
of arrows of length ei for each observation i, now we draw a square with side ei ,
i.e. an area of e2i . We have two apps for you at this point, one where you have
to try and find the best line by choosing b0 and b1 , only focusing on the sum of
errors (and not their square), and a second one focusing on squared errors:
3.2. ORDINARY LEAST SQUARES (OLS) ESTIMATOR 73

10
8
y

6
4
2

0 5 10

Figure 3.2: The best line and its SQUARED errors

library(ScPoEconometrics)
launchApp("reg_simple_arrows")
launchApp("reg_simple") # with squared errors
launchApp("SSR_cone") # visualize the minimzation problem from above!

Most of our apps have an associated about document, which gives extra infor-
mation and explanations. After you have looked at all three apps, we invite you
thus to have a look at the associated explainers by typing
aboutApp("reg_simple_arrows")
aboutApp("reg_simple")
aboutApp("SSR_cone")

3.2 Ordinary Least Squares (OLS) Estimator

The method to compute (or estimate) b0 and b1 we illustrated above is called
Ordinary Least Squares, or OLS. b0 and b1 are therefore also often called the
OLS coefficients. By solving problem (3.6) one can derive an explicit formula
for them:

cov(x, y)
b1 = , (3.7)
var(x)

i.e. the estimate of the slope coefficient is the covariance between x and y divided
74 CHAPTER 3. LINEAR REGRESSION

by the variance of x, both computed from our sample of data. With b1 in hand,
we can get the estimate for the intercept as

b0 = ȳ − b1 x̄. (3.8)

where z̄ denotes the sample mean of variable z. The interpretation of the OLS
slope coefficient b1 is as follows. Given a line as in y = b0 + b1 x,
dy
• b1 = dx measures the change in y resulting from a one unit change in x
• For example, if y is wage and x is years of education, b1 would measure
the effect of an additional year of education on wages.
There is an alternative representation for the OLS slope coefficient which relates
to the correlation coefficient r. Remember from section 2.3 that r = cov(x,y)
sx sy ,
where sz is the standard deviation of variable z. With this in hand, we can
derive the OLS slope coefficient as

cov(x, y)
b1 = (3.9)
var(x)
cov(x, y)
= (3.10)
sx sx
sy
=r (3.11)
sx

In other words, the slope coefficient is equal to the correlation coefficient r times
the ratio of standard deviations of y and x.

3.2.1 Linear Regression without Regressor

There are several important special cases for the linear regression introduced
above. Let’s start with the most obvious one: What is the meaning of running
a regression without any regressor, i.e. without a x? Our line becomes very
simple. Instead of (3.1), we get

y = b0 . (3.12)

This means that our minization problem in (3.6) also becomes very simple: We
only have to choose b0 ! We have

∑
N
2
b0 = arg min [yi − int] ,
int
i=1
3.2. ORDINARY LEAST SQUARES (OLS) ESTIMATOR 75

which is a quadratic equation with a unique optimum such that

1 ∑
N
b0 = yi = y.
N i=1

Least Squares without regressor x estimates the sample mean of the

outcome variable y, i.e. it produces y.

3.2.2 Regression without an Intercept

We follow the same logic here, just that we miss another bit from our initial
equation and the minimisation problem in (3.6) now becomes:

∑
N
2
b1 = arg min [yi − slope xi ] (3.13)
slope
i=1
1
∑N
xi yi x̄ȳ
7→ b1 = N
1
∑N
i=1
= (3.14)
2 x2
N i=1 xi

Least Squares without intercept (i.e. with b0 = 0) is a line that passes

through the origin.

In this case we only get to choose the slope b1 of this anchored line.1 You should
now try out both of those restrictions on our linear model by spending some
time with
launchApp("reg_constrained")

3.2.3 Centering A Regression

By centering or demeaning a regression, we mean to substract from both y and

x their respective averages to obtain ỹi = yi − ȳ and x̃i = xi − x̄. We then run
a regression without intercept as above. That is, we use x̃i , ỹi instead of xi , yi
in (3.14) to obtain our slope estimate b1 :
1 This slope is related to the angle between vectors a = (x, y), and b = (x, 0). Hence, it’s

related to the scalar projection of a on b.

76 CHAPTER 3. LINEAR REGRESSION

1
∑N
N i=1 x̃i ỹi
b1 = 1
∑N 2 (3.15)
N i=1 x̃i
∑ N
i=1 (xi − x̄)(yi −
1
N ȳ)
= ∑N (3.16)
i=1 (xi − x̄)
1 2
N
cov(x, y)
= (3.17)
var(x)

This last expression is identical to the one in (3.7)! It’s the standard OLS
estimate for the slope coefficient. We note the following:

Adding a constant to a regression produces the same result as centering

all variables and estimating without intercept. So, unless all variables
are centered, always include an intercept in the regression.

To get a better feel for what is going on here, you can try this out now by
yourself by typing:
launchApp("demeaned_reg")

3.2.4 Standardizing A Regression

Standardizing a variable z means to demean as above, but in addition to divide

the demeaned value by its own standard deviation. Similarly to what we did
above for centering, we define transformed variables y̆i = yiσ−ȳ
y
and x̆i = xiσ−x̄
x
where σz is the standard deviation of variable z. From here on, you should by
now be used to what comes next! As above, we use x̆i , y̆i instead of xi , yi in
(3.14) to this time obtain:

1
∑N
N i=1 x̆i y̆i
b1 = 1
∑N 2 (3.18)
N i=1 x̆i
1
∑N xi −x̄ yi −ȳ
N i=1 σx σy
=
1
∑N ( xi −x̄ )2 (3.19)
N i=1 σx

Cov(x, y)
= (3.20)
σx σy
= Corr(x, y) (3.21)
3.3. PREDICTIONS AND RESIDUALS 77

After we standardize both y and x, the slope coefficient b1 in the regres-

sion without intercept is equal to the correlation coefficient.

And also for this case we have a practical application for you. Just type this
and play around with the app for a little while!
launchApp("reg_standardized")

3.3 Predictions and Residuals

Now we want to ask how our residuals ei relate to the prediction ∑ yˆi . Let us first
N
think about the average of all predictions yˆi , i.e. the number N1 i=1 yˆi . Let’s
just take (3.2) and plug this into this average, so that we get

1 ∑ 1 ∑
N N
yˆi = b0 + b1 x i (3.22)
N i=1 N i=1

1 ∑
N
= b0 + b1 xi (3.23)
N i=1
= b0 + b1 x̄ (3.24)
(3.25)

But that last line is just equal to the formula for the OLS intercept (3.8), b0 =
ȳ − b1 x̄! That means of course that

1 ∑
N
yˆi = b0 + b1 x̄ = ȳ
N i=1

in other words:

The average of our predictions yˆi is identically equal to the mean of the
outcome y. This implies that the average of the residuals is equal to
zero.

Related to this result, we can show that the prediction ŷ and the residuals are
uncorrelated, something that is often called orthogonality between ŷi and ei .
We would write this as
78 CHAPTER 3. LINEAR REGRESSION

1 ∑ 1 ∑
N N
Cov(ŷ, e) = (ŷi − ȳ)(ei − ē) = (ŷi − ȳ)ei (3.26)
N i=1 N i=1

1 ∑ 1 ∑
N N
= ŷi ei − ȳ ei = 0 (3.27)
N i=1 N i=1

It’s useful to bring back the sample data which generate figure 3.1 at this point
in order to verify these claims:
#OUT> y y_hat error
#OUT> 1 2.09 2.57 -0.48
#OUT> 2 2.79 3.41 -0.62
#OUT> 3 6.49 4.25 2.24
#OUT> 4 1.71 5.10 -3.39
#OUT> 5 9.89 5.94 3.95
#OUT> 6 7.62 6.78 0.83
#OUT> 7 4.86 7.63 -2.77
#OUT> 8 7.38 8.47 -1.09
#OUT> 9 10.63 9.31 1.32
Let’s check that these claims are true in this sample of data. We want that
1. The average of ŷi to be the same as the mean of y
2. The average of the errors should be zero.
3. Prediction and errors should be uncorrelated.
# 1.
[Link](mean(ss$error), 0)

#OUT> [1] TRUE

# 2.
[Link](mean(ss$y_hat), mean(ss$y))

#OUT> [1] TRUE

# 3.
[Link](cov(ss$error,ss$y_hat), 0)

#OUT> [1] TRUE

So indeed we can confirm this result with our test dataset. Great!

3.4 Correlation, Covariance and Linearity

It is important to keep in mind that Correlation and Covariance relate to a lin-
ear relationship between x and y. Given how the regression line is estimated by
3.4. CORRELATION, COVARIANCE AND LINEARITY 79

OLS (see just above), you can see that the regression line inherits this property
from the Covariance. A famous exercise by Francis Anscombe (1973) illustrates
this by constructing 4 different datasets which all have identical linear statis-
tics: mean, variance, correlation and regression line are identical. However, the
usefulness of the statistics to describe the relationship in the data is not clear.

dataset 1 dataset 2
12

12
y1

y2
8

8
6

6
4

4
5 10 15 5 10 15

x1 x2
dataset 3 dataset 4
12

12
y3

y4
8

8
6

6
4

5 10 15 5 10 15

x3 x4

The important lesson from this example is the following:

Always visually inspect your data, and don’t rely exclusively on sum-
mary statistics like mean, variance, correlation and regression line. All
of those assume a linear relationship between the variables in your data.

The mission of Anscombe has been continued recently. As a result of this we

can have a look at the datasauRus package, which pursues Anscbombe’s idea
through a multitude of funny data sets, all with the same linear statistics. Don’t
just compute the covariance, or you might actually end up looking at a Dinosaur!
What? Type this to find out:
launchApp("datasaurus")
aboutApp("datasaurus")

3.4.1 Non-Linear Relationships in Data

Suppose our data now looks like this:

80 CHAPTER 3. LINEAR REGRESSION

30
25
y

20
15
10

50 100 150 200 250 300

Figure 3.3: Best line with non-linear data?

30
25
y

20
15
10

50 100 150 200 250 300

Putting our previous best line defined in equation (3.3) as y = b0 + b1 x + e, we

get something like this:
Somehow when looking at 3.3 one is not totally convinced that the straight line
is a good summary of this relationship. For values x ∈ [50, 120] the line seems to
low, then again too high, and it completely misses the right boundary. It’s easy
to address this shortcoming by including higher order terms of an explanatory
variable. We would modify (3.3) to read now
3.5. ANALYSING V AR(Y ) 81

30
25
y

20
15
10

50 100 150 200 250 300

Figure 3.4: Better line with non-linear data!

yi = b0 + b1 xi + b2 x2i + ei (3.28)

This is a special case of multiple regression, which we will talk about in chapter
4. You can see that there are multiple slope coefficients. For now, let’s just see
how this performs:

3.5 Analysing V ar(y)

Analysis of Variance (ANOVA) refers to a method to decompose variation in one
variable as a function of several others. We can use this idea on our outcome
y. Suppose we wanted to know the variance of y, keeping in mind that, by
definition, yi = ŷi + ei . We would write

V ar(y) = V ar(ŷ + e) (3.29)

= V ar(ŷ) + V ar(e) + 2Cov(ŷ, e) (3.30)
= V ar(ŷ) + V ar(e) (3.31)

We have seen above in 3.3 that the covariance between prediction ŷ and error e
is zero, that’s why we have Cov(ŷ, e) = 0 in (3.31). What this tells us in words
is that we can decompose the variance in the observed outcome y into a part
that relates to variance as explained by the model and a part that comes from
82 CHAPTER 3. LINEAR REGRESSION

unexplained variation. Finally, we know the definition of variance, and can thus
write down the respective formulae for each part:
∑N
• V ar(y) = N1 i=1 (yi − ȳ)2
∑N
• V ar(ŷ) = N1 i=1 (yˆi − ȳ)2 , because the mean of ŷ is ȳ as we know. Finally,
∑N
• V ar(e) = N1 i=1 e2i , because the mean of e is zero.
We can thus formulate how the total variation in outcome y is aportioned be-
tween model and unexplained variation:

The total variation in outcome y (often called SST, or total sum of

squares) is equal to the sum of explained squares (SSE) plus the sum of
residuals (SSR). We have thus SST = SSE + SSR.

3.6 Assessing the Goodness of Fit

In our setup, there exists a convenient measure for how good a particular statis-
tical model fits the data. It is called R2 (R squared), also called the coefficient
of determination. We make use of the just introduced decomposition of variance,
and write the formula as

variance explained SSE SSR

R2 = = =1− ∈ [0, 1] (3.32)
total variance SST SST

It is easy to see that a good fit is one where the sum of explained squares (SSE)
is large relativ to the total variation (SST). In such a case, we observe an R2
close to one. In the opposite case, we will see an R2 close to zero. Notice that a
small R2 does not imply that the model is useless, just that it explains a small
fraction of the observed variation.

3.7 An Example: A Log Wage Equation

Let’s consider the following example concerning wage data collected in the 1976
Current Population Survey in the USA.2 We want to investigate the relationship
between average hourly earnings, and years of education. Let’s start with a plot:
data("wage1", package = "wooldridge") # load data

# a function that returns a plot

plotfun <- function(wage1,log=FALSE,rug = TRUE){
2 This example is close to the vignette of the wooldridge package, whose author I hereby

thank for the excellent work.

3.7. AN EXAMPLE: A LOG WAGE EQUATION 83

y = wage1$wage
if (log){
y = log(wage1$wage)
}
plot(y = y,
x = wage1$educ,
col = "red", pch = 21, bg = "grey",
cex=1.25, xaxt="n", frame = FALSE, # set default x-axis to none
main = ifelse(log,"log(Wages) vs. Education, 1976","Wages vs. Education, 1976"),
xlab = "years of education",
ylab = ifelse(log,"Log Hourly wages","Hourly wages"))
axis(side = 1, at = c(0,6,12,18)) # add custom ticks to x axis
if (rug) rug(wage1$wage, side=2, col="red") # add `rug` to y axis
}

par(mfcol = c(2,1)) # set up a plot with 2 panels

# plot 1: standard scatter plot
plotfun(wage1)

# plot 2: add a panel with histogram+density

hist(wage1$wage,prob = TRUE, col = "grey", border = "red",
main = "Histogram of wages and Density",xlab = "hourly wage")
lines(density(wage1$wage), col = "black", lw = 2)

Looking at the top panel of figure 3.5, you notice two things: From the red
ticks on the y axis, you see that wages are very concentrated at around 5 USD
per hour, with fewer and fewer observations at higher rates; and second, that
it seems that the hourly wage seems to increase with higher education levels.
The bottom panel reinforces the first point, showing that the estimated pdf
(probability density function) shown as a black line has a very long right tail:
there are always fewer and fewer, but always larger and larger values of hourly
wage in the data.

You have seen this shape of a distribution in the tutorial for chapter 2
already! Do you remember the name of this particular shape of a distribu-
tion? (why not type ScPoEconometrics::runTutorial('chapter2'))
to check?

Let’s run a first regression on this data to generate some intution:

wagei = b0 + b1 educi + ei (3.33)

We use the lm function for this purpose as follows:

84 CHAPTER 3. LINEAR REGRESSION

Wages vs. Education, 1976

20
Hourly wages

10
5
0

0 6 12 18

years of education

Histogram of wages and Density

Density

0.10
0.00

0 5 10 15 20 25

hourly wage

Figure 3.5: Wages vs Education from the wooldridge dataset wage1.

3.7. AN EXAMPLE: A LOG WAGE EQUATION 85

Wages vs. Education, 1976

25
20
Hourly wages

15
10
5
0

0 6 12 18

years of education

Figure 3.6: Wages vs Education from the wooldridge dataset wage1, with re-
gression

hourly_wage <- lm(formula = wage ~ educ, data = wage1)

and we can add the resulting regression line to our above plot:
plotfun(wage1)
abline(hourly_wage, col = 'black', lw = 2) # add regression line

The hourly_wage object contains the results of this estimation. We can get a
summary of those results with the summary method:
summary(hourly_wage)

#OUT>
#OUT> Call:
#OUT> lm(formula = wage ~ educ, data = wage1)
#OUT>
#OUT> Residuals:
#OUT> Min 1Q Median 3Q Max
#OUT> -5.3396 -2.1501 -0.9674 1.1921 16.6085
#OUT>
#OUT> Coefficients:
#OUT> Estimate Std. Error t value Pr(>|t|)
#OUT> (Intercept) -0.90485 0.68497 -1.321 0.187
86 CHAPTER 3. LINEAR REGRESSION

#OUT> educ 0.54136 0.05325 10.167 <2e-16 ***

#OUT> ---
#OUT> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#OUT>
#OUT> Residual standard error: 3.378 on 524 degrees of freedom
#OUT> Multiple R-squared: 0.1648, Adjusted R-squared: 0.1632
#OUT> F-statistic: 103.4 on 1 and 524 DF, p-value: < 2.2e-16
The main interpretation of this table can be read off the column labelled Esti-
mate, reporting estimated coefficients b0 , b1 :
1. With zero year of education, the hourly wage is about -0.9 dollars per hour
(row named (Intercept))
2. Each additional year of education increase hourly wage by 54 cents. (row
named educ)
3. For example, for 15 years of education, we predict roughly -0.9 + 0.541 *
15 = 7.215 dollars/h.

3.8 Scaling Regressions

Regression estimates (b0 , b1 ) are in the scale of the data. The actual
value of the estimates will vary, if we change the scale of the data. The
overall fit of the model to the data would not change, however, so that
the R2 statistic would be constant.

Suppose we wanted to use the above estimates to report the effect of years of
education on annual wages instead of hourly ones. Let’s assume we have full-
time workers, 7h per day, 5 days per week, 45 weeks per year. Calling this factor
δ = 7 × 5 × 45 = 1575, we have that x dollars per hour imply x × δ = x × 1575
dollars per year.
What would be the effect of using ỹ = wage × 1575 instead of y = wage as
outcome variable on our regression coefficients b0 and b1 ? Well, let’s try!
% Table created by stargazer v.5.2.2 by Marek Hlavac, Harvard University. E-
mail: hlavac at [Link] % Date and time: Thu, Oct 10, 2019 - 11:40:26
Let’s call the coefficients in the column labelled (1) as b0 and b1 , and let’s call the
ones in column (2) b∗0 and b∗1 . In column (1) we see that another year increaeses
hourly wage by 0.54 dollars, as before. In column (2), the corresponding number
is 852.64, i.e. another year of education will increase annual wages by 852.64
dollars, on average. Notice however, that b0 × δ = −0.9 × 1575 = −1425.14 = b∗0
and that b1 × δ = 0.54 × 1575 = 852.64 = b∗1 , that is we just had to multiply
both coefficients by the scaling factor applied to original outcome y to obtain
our new coefficients b∗0 and b∗1 ! Also, observe that the R2 s of both regressions
3.8. SCALING REGRESSIONS 87

Table 3.1: Effect of Scaling on Coefficients

Dependent variable:
wage annual_wage
(1) (2)
educ 0.541∗∗∗ 852.641∗∗∗
(0.053) (83.866)

Constant −0.905 −1,425.141

(0.685) (1,078.824)

Observations 526 526

R2 0.165 0.165
Adjusted R2 0.163 0.163
Residual Std. Error (df = 524) 3.378 5,320.963
F Statistic (df = 1; 524) 103.363∗∗∗ 103.363∗∗∗
∗ ∗∗ ∗∗∗
Note: p<0.1; p<0.05; p<0.01

are identical! So, really, we did not have to run the regression in column (2) at
all to make this change: multiplying all coefficients through by δ is enough in
this case. We keep the identically same fit to the data.

Rescaling the regressors x is slightly different, but it’s easy to work out how
different, given the linear nature of the covariance operator, which is part of the
OLS estimator. Suppose we rescale x by the number c. Then, using the OLS
formula in (3.7), we see that we get new slope coefficient b∗1 via

Cov(cx, y)
b∗1 = (3.34)
V ar(cx)
cCov(x, y)
= 2 (3.35)
c V ar(x)
1
= b1 . (3.36)
c
88 CHAPTER 3. LINEAR REGRESSION

As for the intercept, and by using (3.8)

1 ∑
N
b∗0 = ȳ − b∗1 c · xi (3.37)
N i=1

c ∑
N
= ȳ − b∗1 xi (3.38)
N i=1
1
= ȳ − b1 c ∗ x̄ (3.39)
c
= ȳ − b1 ∗ x̄ (3.40)
= b0 (3.41)

That is, we change the slope by the inverse of the scaling factor applied to
regressor x, but the intercept is unaffected from this. You should play around
for a while with our rescaling app to get a feeling for this:
library(ScPoEconometrics)
launchApp('Rescale')

3.9 A Particular Rescaling: The log Transform

The natural logarithm is a particularly important transformation that we of-

ten encounter in economics. Why would we transform a variable with the log
function to start with?
1. Several important economic variables (like wages, city size, firm size, etc)
are approximately log-normally distributed. By transforming them with
the log, we obtain an approximately normally distributed variable, which
has desirable properties for our regression.
2. Applying the log reduces the impact of outliers.
3. The transformation allows for a convenient interpretation in terms of per-
centage changes of the outcome variable.
Let’s investigate this issue in our running example by transforming the wage
data above. Look back at the bottom panel of figure 3.5: Of course you saw
immediately that this looked a lot like a log-normal distribution, so point 1.
above applies. We modify the left hand side of equation (3.33):

log(wagei ) = b0 + b1 educi + ei (3.42)

Let’s use the update function to modify our previous regression model:
log_hourly_wage = update(hourly_wage, log(wage) ~ ., data = wage1)
3.9. A PARTICULAR RESCALING: THE LOG TRANSFORM 89

The update function takes an existing lm object, like hourly_wage here, and
updates the formula. Here the . on the right hand side means leave unchanged
(so the RHS stays unchanged). How do our pictures change?
par(mfrow = c(1,2))

plotfun(wage1,rug = FALSE)
abline(hourly_wage, col = 'black', lw = 2) # add regression line

plotfun(wage1,log = TRUE, rug = FALSE)

abline(log_hourly_wage, col = 'black', lw = 2) # add regression line

Wages vs. Education, 1976 log(Wages) vs. Education, 1976

25

3
20

Log Hourly wages

Hourly wages

2
15

1
10
5

0
0

0 6 12 18 0 6 12 18

years of education years of education

par(mfrow = c(1,1))

It looks as if the regression line has the same slope, but beware of the different
scales of the y-axis! You can clearly see that all y-values have been compressed
by the log transformation. The log case behaves differently from our scaling by
a constant number case above because it is a nonlinear function. Let’s compare
the output between both models:

% Table created by stargazer v.5.2.2 by Marek Hlavac, Harvard University. E-

mail: hlavac at [Link] % Date and time: Thu, Oct 10, 2019 - 11:40:26

The interpretation of the transformed model in column (2) is now the following:
90 CHAPTER 3. LINEAR REGRESSION

Table 3.2: Log Transformed Equation

Dependent variable:
wage NA
(1) (2)
educ 0.541∗∗∗ 0.083∗∗∗
(0.053) (0.008)

Constant −0.905 0.584∗∗∗

(0.685) (0.097)

Observations 526 526

R2 0.165 0.186
Adjusted R2 0.163 0.184
Residual Std. Error (df = 524) 3.378 0.480
F Statistic (df = 1; 524) 103.363∗∗∗ 119.582∗∗∗
∗ ∗∗ ∗∗∗
Note: p<0.1; p<0.05; p<0.01

We call a regression of the form log(y) = b0 + b1 x + u a log-level speci-

fication, because we regressed the log of a variable on the level (i.e not
the log!) of another variable. Here, the impact of increasing x by one
unit is to increase y by 100 × b1 percent. In our example: an additional
year of education will increase hourly wages by 8.3%. Notice that this is
very different from saying …increases log hourly wages by 8.3%, which is
wrong.

Notice that the R2 slightly improved, so have a better fit to the data. This is
due the fact that the log compressed large outlier values. Whether we apply
the log to left or right-hand side variables makes a difference, as outlined in this
important table:
Common Regression Specifications

Specification Outcome Var Regressor Interpretation of b1 Comment

Level-level y x ∆y = b1 ∆x Standard
Level-log y log(x) ∆y = b1 ∆x less frequent
Log-level log(y) x %∆y = (100b1 )∆ x Semi-elasticity
Log-Log log(y) log(x) %∆y = %∆ b_1 x Elasticity

You may remember from your introductory micro course what the definition
of the elasticity of y with respect to x is: This number tells us by how many
percent y will change, if we change x by one percent. Let’s look at another
example from the wooldridge package of datasets, this time concerning CEO
salaries and their relationship with company sales.
data("ceosal1", package = "wooldridge")
par(mfrow = c(1,2))
plot(salary ~ sales, data = ceosal1, main = "Sales vs Salaries",xaxt = "n",frame = FALS
axis(1, at = c(0,40000, 80000))
rug(ceosal1$salary,side = 2)
rug(ceosal1$sales,side = 1)
3.9. A PARTICULAR RESCALING: THE LOG TRANSFORM 91

15000 Sales vs Salaries Log(Sales) vs Log(Salaries)

9
10000

log(salary)

8
salary

5000

7
6
0

0 40000 80000 5 6 7 8 9 10

sales log(sales)

Figure 3.7: The effect of log-transforming highly skewed data.

Refering back at table 3.9, here we have a log-log specification. Therefore we

interpret this regression as follows:

In a log-log equation, the slope coefficient b1 is the elasticity of y with

respect to changes in x. Here: A 1% increase in sales is associated to a
0.26% increase in CEO salaries. Note, again, that there is no log in this
statement.
92 CHAPTER 3. LINEAR REGRESSION
Chapter 4

Multiple Regression

We can extend the discussion from chapter 3 to more than one explanatory
variable. For example, suppose that instead of only x we now had x1 and x2
in order to explain y. Everything we’ve learned for the single variable case
applies here as well. Instead of a regression line, we now get a regression plane,
i.e. an object representable in 3 dimenions: (x1 , x2 , y). As an example, suppose
we wanted to explain how many miles per gallon (mpg) a car can travel as a
function of its horse power (hp) and its weight (wt). In other words we want to
estimate the equation

mpgi = b0 + b1 hpi + b2 wti + ei (4.1)

on our built-in dataset of cars (mtcars):

head(subset(mtcars, select = c(mpg,hp,wt)))

#OUT> mpg hp wt
#OUT> Mazda RX4 21.0 110 2.620
#OUT> Mazda RX4 Wag 21.0 110 2.875
#OUT> Datsun 710 22.8 93 2.320
#OUT> Hornet 4 Drive 21.4 110 3.215
#OUT> Hornet Sportabout 18.7 175 3.440
#OUT> Valiant 18.1 105 3.460
How do you think hp and wt will influence how many miles per gallon of gasoline
each of those cars can travel? In other words, what do you expect the signs of
b1 and b2 to be?
With two explanatory variables as here, it is still possible to visualize the re-
gression plane, so let’s start with this as an answer. The OLS regression plane
through this dataset looks like in figure 4.1:

93
94 CHAPTER 4. MULTIPLE REGRESSION

WebGL is not
supported by
data
your browser -
visit
[Link]
for more info
Figure 4.1: Multiple Regression - a plane in 3D. The red lines indicate the
residual for each observation.

This visualization shows a couple of things: the data are shown with red points
and the grey plane is the one resulting from OLS estimation of equation (4.1).
You should realize that this is exactly the same story as told in figure 3.1 - just
in three dimensions!
Furthermore, multiple regression refers the fact that there could be more than
two regressors. In fact, you could in principle have K regressors, and our theory
developed so far would still be valid:

ŷi = b0 + b1 x1i + b2 x2i + · · · + bK xKi (4.2)

ei = yi − ŷi (4.3)

Just as before, the least squares method chooses numbers (b0 , b1 , . . . , bK ) to as

to minimize SSR, exactly as in the minimization problem for the one regressor
case seen in (3.6).

4.1 All Else Equal

We can see from the above plot that cars with more horse power and greater
weight, in general travel fewer miles per gallon of combustible. Hence, we observe
4.1. ALL ELSE EQUAL 95

a plane that is downward sloping in both the weight and horse power directions.
Suppose now we wanted to know impact of hp on mpg in isolation, so as if we
could ask

Keeping the value of wt fixed for a certain car, what would be the impact
on mpg be if we were to increase only its hp? Put differently, keeping
all else equal, what’s the impact of changing hp on mpg?

We ask this kind of question all the time in econometrics. In figure 4.1 you clearly
see that both explanatory variables have a negative impact on the outcome of
interest: as one increases either the horse power or the weight of a car, one finds
that miles per gallon decreases. What is kind of hard to read off is how negative
an impact each variable has in isolation.
As a matter of fact, the kind of question asked here is so common that it has
got its own name: we’d say “ceteris paribus, what is the impact of hp on mpg?”.
ceteris paribus is latin and means the others equal, i.e. all other variables fixed.
In terms of our model in (4.1), we want to know the following quantity:

∂mpgi
= b1 (4.4)
∂hpi

The ∂ sign denotes a partial derivative of the function describing mpg with
respect to the variable hp. It measures how the value of mpg changes, as we
change the value of hp ever so slightly. In our context, this means: keeping
all other variables fixed, what is the effect of hp on mpg?. We call the value
of coefficient b1 therefore also the partial effect of hp on mpg. In terms of our
dataset, we use R to run the following multiple regression:
#OUT>
#OUT> Call:
#OUT> lm(formula = mpg ~ wt + hp, data = mtcars)
#OUT>
#OUT> Residuals:
#OUT> Min 1Q Median 3Q Max
#OUT> -3.941 -1.600 -0.182 1.050 5.854
#OUT>
#OUT> Coefficients:
#OUT> Estimate Std. Error t value Pr(>|t|)
#OUT> (Intercept) 37.22727 1.59879 23.285 < 2e-16 ***
#OUT> wt -3.87783 0.63273 -6.129 1.12e-06 ***
#OUT> hp -0.03177 0.00903 -3.519 0.00145 **
#OUT> ---
#OUT> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#OUT>
#OUT> Residual standard error: 2.593 on 29 degrees of freedom
96 CHAPTER 4. MULTIPLE REGRESSION

#OUT> Multiple R-squared: 0.8268, Adjusted R-squared: 0.8148

#OUT> F-statistic: 69.21 on 2 and 29 DF, p-value: 9.109e-12
From this table you see that the coefficient on wt has value -3.87783. You can
interpret this as follows:

Holding all other variables fixed at their observed values - or ceteris

paribus - a one unit increase in wt implies a -3.87783 units change in
mpg. In other words, increasing the weight of a car by 1000 pounds
(lbs), will lead to 3.88 miles less travelled per gallon. Similarly, a car
with one additional horse power means that we will travel 0.03177 fewer
miles per gallon of gasoline, all else (i.e. wt) equal.

4.2 Multicolinearity
One important requirement for multiple regression is that the data be not lin-
early dependent: Each variable should provide at least some new information
for the outcome, and it cannot be replicated as a linear combination of other
variables. Suppose that in the example above, we had a variable wtplus defined
as wt + 1, and we included this new variable together with wt in our regression.
In this case, wtplus provides no new information. It’s enough to know wt, and
add 1 to it. In this sense, wt_plus is a redundant variable and should not be
included in the model. Notice that this holds only for linearly dependent vari-
ables - nonlinear transformations (like for example wt2 ) are exempt from this
rule. Here is why:

y = b0 + b1 wt + b2 wtplus + e (4.5)
= b0 + b1 wt + b2 (wt + 1) + e (4.6)
= (b0 + b2 ) + wt(b1 + b2 ) + e (4.7)

This shows that we cannot identify the regression coefficients in case of linearly
dependent data. Variation in the variable wt identifies a different coefficient,
say γ = b1 + b2 , from what we actually wanted: separate estimates for b1 , b2 .

We cannot have variables which are linearly dependent, or perfectly colin-

ear. This is known as the rank condition. In particular, the condition
dictates that we need at least N ≥ K + 1, i.e. more observations than
coefficients. The greater the degree of linear dependence amongst our
explanatory variables, the less information we can extract from them,
and our estimates becomes less precise.
4.3. LOG WAGE EQUATION 97

WebGL is not
supported by
data
your browser -
visit
[Link]
for more info
Figure 4.2: Log wages vs education and experience in 3D.

4.3 Log Wage Equation

Let’s go back to our previous example of the relationship between log wages and
education. How does this relationship change if we also think that experience in
the labor market has an impact, next to years of education? Here is a picture:

Let’s add even more variables! For instance, what’s the impact of experience in
the labor market, and time spent with the current employer? Let’s first look at
how those variables co-vary with each other:
cmat = round(cor(subset(wage1,select = c(lwage,educ,exper,tenure))),2) # correlation matrix
corrplot::corrplot(cmat,type = "upper",method = "ellipse")
98 CHAPTER 4. MULTIPLE REGRESSION

tenure
lwage

exper
educ
1

0.8
lwage
0.6

0.4
educ
0.2

−0.2
exper
−0.4

−0.6
tenure
−0.8

−1
The way to read the so-called correlation plot in figure ?? is straightforward:
each row illustrates the correlation of a certain variable with the other variables.
In this example both the shape of the ellipse in each cell as well as their color
coding tell us how strongly two variables correlate. Let us put this into a
regression model now:
educ_only <- lm(lwage ~ educ , data = wage1)
educ_exper <- lm(lwage ~ educ + exper , data = wage1)
log_wages <- lm(lwage ~ educ + exper + tenure, data = wage1)
stargazer::stargazer(educ_only, educ_exper, log_wages,type = if (knitr:::is_latex_outpu

% Table created by stargazer v.5.2.2 by Marek Hlavac, Harvard University. E-

mail: hlavac at [Link] % Date and time: Thu, Oct 10, 2019 - 11:40:34
Column (1) refers to model (3.42) from the previous chapter, where we only had
educ as a regressor: we obtain an R2 of 0.186. Column (2) is the model that
generated the plane in figure 4.2 above. (3) is the model with three regressors.
You can see that by adding more regressors, the quality of our fit increases, as
more of the variation in y is now accounted for by our model. You can also
see that the values of our estimated coefficients keeps changing as we move
4.3. LOG WAGE EQUATION 99

Table 4.1:

Dependent variable:
lwage
(1) (2) (3)
educ 0.083∗∗∗ 0.098∗∗∗ 0.092∗∗∗
(0.008) (0.008) (0.007)

exper 0.010∗∗∗ 0.004∗∗

(0.002) (0.002)

tenure 0.022∗∗∗
(0.003)

Constant 0.584∗∗∗ 0.217∗∗ 0.284∗∗∗

(0.097) (0.109) (0.104)

Observations 526 526 526

R2 0.186 0.249 0.316
Adjusted R2 0.184 0.246 0.312
Residual Std. Error 0.480 (df = 524) 0.461 (df = 523) 0.441 (df = 522)
F Statistic 119.582∗∗∗ (df = 1; 524) 86.862∗∗∗ (df = 2; 523) 80.391∗∗∗ (df = 3; 522)
∗ ∗∗ ∗∗∗
Note: p<0.1; p<0.05; p<0.01
100 CHAPTER 4. MULTIPLE REGRESSION

from left to right across the columns. Given the correlation structure shown
in figure ??, it is only natural that this is happening: We see that educ and
exper are negatively correlated, for example. So, if we omit exper from the
model in column (1), educ will reflect part of this correlation with exper by a
lower estimated value. By directly controlling for exper in column (2) we get
an estimate of the effect of educ net of whatever effect exper has in isolation
on the outcome variable. We will come back to this point later on.

4.4 How To Make Predictions

So suppose we have a model like

lwage = b0 + b1 (educ) + b2 (exper) + b3 (tenure) + ϵ

How could we use this to make a prediction of log wages, given some new data?
Remember that the OLS procedure gives us estimates for the values b0 , b1 , b2 , b3 .
With those in hand, it is straightforward to make a prediction about the condi-
tional mean of the outcome - just plug in the desired numbers for educ,exper
and tenure. Suppose you want to know what the mean of lwage is conditional
on educ = 10,exper=4 and tenure = 2. You’d do

E[lwage|educ = 10, exper = 4, tenure = 2] = b0 + b1 10 + b2 4 + b3 2 (4.8)

= 1.27. (4.9)

I computed the last line directly with

x = c(1,10,4,2) # 1 for intercept
pred = coef(log_wages) %*% x

but R has a more complete prediction interface, using the function predict. For
starters, you can predict the model on all data points which were contained in
the dataset we used for estimation, i.e. wage1 in our case:
head(predict(log_wages)) # first 6 observations of wage1 as predicted by our model

#OUT> 1 2 3 4 5 6
#OUT> 1.304921 1.523506 1.304921 1.819802 1.461690 1.970451
Often you want to add that prediction to the original dataset:
wage_prediction = cbind(wage1, prediction = predict(log_wages))
head(wage_prediction[, c("lwage","educ","exper","tenure","prediction")])

#OUT> lwage educ exper tenure prediction

#OUT> 1 1.131402 11 2 0 1.304921
4.4. HOW TO MAKE PREDICTIONS 101

#OUT> 2 1.175573 12 22 2 1.523506

#OUT> 3 1.098612 11 2 0 1.304921
#OUT> 4 1.791759 8 44 28 1.819802
#OUT> 5 1.667707 12 7 2 1.461690
#OUT> 6 2.169054 16 9 8 1.970451
You’ll remember that we called the distance in prediction and observed out-
come our residual e. Well here this is just lwage - prediction. Indeed, e is
such an important quantity that R has a convenient method to compute y − ŷ
from an lm object directly - the method resid. Let’s add another column to
wage_prediction:
wage_prediction = cbind(wage_prediction, residual = resid(log_wages))
head(wage_prediction[, c("lwage","educ","exper","tenure","prediction","residual")])

#OUT> lwage educ exper tenure prediction residual

#OUT> 1 1.131402 11 2 0 1.304921 -0.17351850
#OUT> 2 1.175573 12 22 2 1.523506 -0.34793289
#OUT> 3 1.098612 11 2 0 1.304921 -0.20630832
#OUT> 4 1.791759 8 44 28 1.819802 -0.02804286
#OUT> 5 1.667707 12 7 2 1.461690 0.20601725
#OUT> 6 2.169054 16 9 8 1.970451 0.19860271
Using the data in wage_prediction, you should now check for yourself what
we already know about ŷ and e from section 3.3:
1. What is the average of the vector residual?
2. What is the average of prediction?
3. How does this compare to the average of the outcome lwage?
4. What is the correlation between prediction and residual?
102 CHAPTER 4. MULTIPLE REGRESSION
Chapter 5

Categorial Variables

Up until now, we have encountered only examples with continuous variables x

and y, that is, x, y ∈ R, so that a typical observation could have been (yi , xi ) =
(1.5, 5.62). There are many situations where it makes sense to think about
the data in terms of categories, rather than continuous numbers. For example,
whether an observation i is male or female, whether a pixel on a screen is black
or white, and whether a good was produced in France, Germany, Italy, China
or Spain are all categorical classifications of data.
Probably the simplest type of categorical variable is the binary, boolean, or just
dummy variable. As the name suggests, it can take on only two values, 0 and 1,
or TRUE and FALSE.

5.1 The Binary Regressor Case

Even though this is an extremely parsimonious way of encoding that, it is a
very powerful tool that allows us to represent that a certain observation i is a
member of a certain category j. For example, let’s imagine we have income
data on males and females, and we would create a variable called [Link] that
is TRUE whenever i is male, FALSE otherwise, and similarly for women. For
example, to encode whether subject i is male, one could do this:

{
1 if i is male
[Link] = ,
0 if i is not male.

and similarly for females, we’d have

103
104 CHAPTER 5. CATEGORIAL VARIABLES

{
1 if i is female
[Link] =
0 if i is not female.

By definition, we have just introduced a linear dependence into our dataset.

It will always be true that [Link] + [Link] = 1. This is because dummy
variables are based on data being mutually exclusively categorized - here, you
are either male or female.1 This should immediately remind you of section 4.2
where we introduced multicolinearity. A regression of income on both of our
variables like this

yi = b0 + b1 [Link] + b2 [Link] + ei
would be invalid because of perfect colinearity between [Link] and [Link] .
The solution to this is pragmatic and simple:

In dummy variable regressions, we remove one category from the regres-

sion (for example here: [Link]) and call it the reference category. The
effect of being male is absorbed in the intercept. The coefficient on
the remaining categories measures the difference in mean outcome with
respect to the reference category.

Now let’s try this out. We start by creating the female indicator as above,
{
1 if i is female
[Link] =
0 if i is not female.
and let’s suppose that yi is a measure of i’s annual labor income. Our model is

yi = b0 + b1 [Link] + ei (5.1)

and here is how we estimate this in R:

# x = sample(x = c(0, 1), size = n, replace = T)
dta$[Link] = factor(x) # convert x to factor
dummy_reg = lm(y~[Link],dta)
summary(dummy_reg)

#OUT>
#OUT> Call:
#OUT> lm(formula = y ~ [Link], data = dta)
1 There are transgender individuals where this example will not apply.
5.1. THE BINARY REGRESSOR CASE 105

#OUT>
#OUT> Residuals:
#OUT> Min 1Q Median 3Q Max
#OUT> -2.3901 -0.6565 0.1612 0.6846 2.7850
#OUT>
#OUT> Coefficients:
#OUT> Estimate Std. Error t value Pr(>|t|)
#OUT> (Intercept) 2.0847 0.2396 8.700 1.97e-11 ***
#OUT> is.female1 -3.0631 0.3202 -9.566 1.06e-12 ***
#OUT> ---
#OUT> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#OUT>
#OUT> Residual standard error: 1.124 on 48 degrees of freedom
#OUT> Multiple R-squared: 0.6559, Adjusted R-squared: 0.6488
#OUT> F-statistic: 91.5 on 1 and 48 DF, p-value: 1.061e-12

Notice that R displays the level of the factor to which coefficient b1 belongs here,
i.e. is.female1 means this coefficient is on level [Link] = 1 - the reference
level is [Link] = 0, and it has no separate coefficient. Also interesting is
that b1 is equal to the difference in conditional means between male and female

b1 = E[y|[Link] = 1] − E[y|[Link] = 0] = −3.0631.

A dummy variable measures the difference or the offset in the mean of the
response variable, E[y], conditional on x belonging to some category -
relative to a baseline category. In our artificial example, the coefficient
b1 informs us that women earn on average 3.756 units less than men.

It is instructive to reconsider this example graphically:

In figure 5.1 we see that this regression simplifies to the straight line connecting
the mean, or the expected value of y when [Link] = 0, i.e. E[y|[Link] =
0], to the mean when [Link] = 1, i.e. E[y|[Link] = 1]. It is useful to
remember that the unconditional mean of y, i.e. E[y], is going to be the result
of regressing y only on an intercept, illustrated by the blue line. This line will
always lie in between both conditional means. As indicated by the red arrow,
the estimate of the coefficient on the dummy, b1 , is equal to the difference in
conditional means for both groups. You should look at our app now to deepen
your understanding of what’s going on here:
library(ScPoEconometrics)
launchApp("reg_dummy")
106 CHAPTER 5. CATEGORIAL VARIABLES

4
E[Y | [Link] = 0]
2 b1= −3.06
y

E[Y | [Link] = 1]
−2

−1 0 1 2

Figure 5.1: regressing y ∈ R on [Link] ∈ {0, 1}. The blue line is E[y], the
red arrow is the size of b1 . Which is the same as the slope of the regression line
in this case and the difference in conditional means!

5.2 Dummy and Continuous Variables

What happens if there are more predictors than just the dummy variable in a
regression? For example, what if instead we had

yi = b0 + b1 [Link] + b2 experi + ei (5.2)

where experi would measure years of experience in the labor market? As above,
the dummy variable acts as an intercept shifter. We have

{
b0 + b1 + b2 × experi + ei if [Link]=1
yi = (5.3)
b0 + + b2 × experi + ei if [Link]=0

so that the intercept is b0 + b1 for women but b0 for men. We will see this in
the real-world example below, but for now let’s see the effect of switching the
dummy on and off in this app:
library(ScPoEconometrics)
launchApp("reg_dummy_example")
5.3. CATEGORICAL VARIABLES IN R: FACTOR 107

5.3 Categorical Variables in R: factor

R has extensive support for categorical variables built-in. The relevant data
type representing a categorical variable is called factor. We encountered them
as basic data types in section 1.8 already, but it is worth repeating this here.
We have seen that a factor categorizes a usually small number of numeric values
by labels, as in this example which is similar to what I used to create regressor
[Link] for the above regression:
[Link] = factor(x = c(0,1,1,0), labels = c(FALSE,TRUE))
[Link]

#OUT> [1] FALSE TRUE TRUE FALSE

#OUT> Levels: FALSE TRUE
You can see the result is a vector object of type factor with 4 entries, whereby 0
is represented as FALSE and 1 as TRUE. An other example could be if we wanted
to record a variable sex instead, and we could do
sex = factor(x = c(0,1,1,0), labels = c("male","female"))
sex

#OUT> [1] male female female male

#OUT> Levels: male female
You can see that this is almost identical, just the labels are different.

5.3.1 More Levels

We can go beyond binary categorical variables such as TRUE vs FALSE. For exam-
ple, suppose that x measures educational attainment, i.e. it is now something
like xi ∈ {high school,some college,BA,MSc}. In R parlance, high school, some
college, BA, MSc are the levels of factor x. A straightforward extension of the
above would dictate to create one dummy variable for each category (or level),
like

[Link] = 1[xi == high school]

[Link] = 1[xi == some college]
[Link] = 1[xi == BA]
[Link] = 1[xi == MSc]

but you can see that this is cumbersome. There is a better solution for us
available:
108 CHAPTER 5. CATEGORIAL VARIABLES

factor(x = c(1,1,2,4,3,4),labels = c("HS","someCol","BA","MSc"))

#OUT> [1] HS HS someCol MSc BA MSc

#OUT> Levels: HS someCol BA MSc

Notice here that R will apply the labels in increasing order the way you supplied
it (i.e. a numerical value 4 will correspond to “MSc”, no matter the ordering in
x.)

5.3.2 Log Wages and Dummies

The above developed factor terminology fits neatly into R’s linear model fitting
framework. Let us illustrate the simplest use by way of example.

Going back to our wage example, let’s say that a worker’s wage depends on
their education as well as their sex:

ln wi = b0 + b1 educi + b2 f emalei + ei (5.4)

data("wage1", package = "wooldridge")

wage1$female = [Link](wage1$female) # convert 0-1 to factor
lm_w = lm(lwage ~ educ, data = wage1)
lm_w_sex = lm(lwage ~ educ + female, data = wage1)
stargazer::stargazer(lm_w,lm_w_sex,type = if (knitr:::is_latex_output()) "latex" else "

% Table created by stargazer v.5.2.2 by Marek Hlavac, Harvard University. E-

mail: hlavac at [Link] % Date and time: Thu, Oct 10, 2019 - 11:40:34

We know the results from column (1) very well by now. How does the relation-
ship change if we include the female indicator? Remember from above that
female is a factor with two levels, 0 and 1, where 1 means that’s a female.
We see in the above output that R included a regressor called female1. This is
a combination of the variable name female and the level which was included
in the regression. In other words, R chooses a reference category (by default
the first of all levels by order of appearance), which is excluded - here this is
female==0. The interpretation is that b2 measures the effect of being female rel-
ative to being male. R automatically creates a dummy variable for each potential
level, excluding the first category.

Figure 5.2 illustrates this. The left panel is our previous model. The right panel
adds the female dummy. You can see that both male and female have the same
upward sloping regression line. But you can also see that there is a parallel
downward shift from male to female line. The estimate of b2 = −0.36 is the size
of the downward shift.
5.3. CATEGORICAL VARIABLES IN R: FACTOR 109

Table 5.1:

Dependent variable:
lwage
(1) (2)
∗∗∗
educ 0.083 0.077∗∗∗
(0.008) (0.007)

female1 −0.361∗∗∗
(0.039)

Constant 0.584∗∗∗ 0.826∗∗∗

(0.097) (0.094)

Observations 526 526

R2 0.186 0.300
Adjusted R2 0.184 0.298
Residual Std. Error 0.480 (df = 524) 0.445 (df = 523)
F Statistic 119.582∗∗∗ (df = 1; 524) 112.189∗∗∗ (df = 2; 523)
∗ ∗∗ ∗∗∗
Note: p<0.1; p<0.05; p<0.01

3 3

2 2
female
lwage

0
1 1 1

0 0

0 5 10 15 0 5 10 15
educ educ

Figure 5.2: log wage vs educ. Right panel with female dummy.
110 CHAPTER 5. CATEGORIAL VARIABLES

5.4 Interactions

Sometimes it is useful to let the slope of a certain variable to be dependent on

the value of another regressor. For example consider a model for the sales prices
of houses, where area is the livable surface of the property, and age is its age:

log(price) = b0 + b1 area + b2 age + b3 (area × age) + e (5.5)

In that model, the partial effect of area on log(price), keeping all other vari-
ables fixed, is

∂ log(price)
= b1 + b3 (age) (5.6)
∂area

If we find that b3 > 0 in a regression, we conclude that the size of a house values
more in older houses. We call b3 the interaction effect between area and age.
Let’s look at that regression model now.
data(hprice3, package = "wooldridge")
summary(lm(lprice ~ area*age, data = hprice3))

#OUT>
#OUT> Call:
#OUT> lm(formula = lprice ~ area * age, data = hprice3)
#OUT>
#OUT> Residuals:
#OUT> Min 1Q Median 3Q Max
#OUT> -1.27226 -0.16538 -0.00298 0.20673 0.83985
#OUT>
#OUT> Coefficients:
#OUT> Estimate Std. Error t value Pr(>|t|)
#OUT> (Intercept) 1.071e+01 6.633e-02 161.448 < 2e-16 ***
#OUT> area 3.647e-04 2.875e-05 12.686 < 2e-16 ***
#OUT> age -7.377e-03 1.358e-03 -5.434 1.1e-07 ***
#OUT> area:age 9.168e-07 4.898e-07 1.872 0.0622 .
#OUT> ---
#OUT> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#OUT>
#OUT> Residual standard error: 0.2925 on 317 degrees of freedom
#OUT> Multiple R-squared: 0.5586, Adjusted R-squared: 0.5545
#OUT> F-statistic: 133.7 on 3 and 317 DF, p-value: < 2.2e-16
In this instance, we see that indeed there is a small positive interaction between
area and age on the sales price: even though age in isolation decreases the sales
value, bigger houses command a small premium if they are older.
5.4. INTERACTIONS 111

5.4.1 Interactions with Dummies: Differential Slopes

It is straightforward to extend the interactions logic to allow not only for differ-
ent intercepts, but also different slopes for each subgroup in a dataset. Let’s go
back to our dataset of wages from section 5.3.2 above. Now that we know how
to create and interaction between two variables, we can easily modify equation
(5.4) like this:

ln w = b0 + b1 female + b2 educ + b3 (female × educ) + e (5.7)

The only peculiarity here is that female is a factor with levels 0 and 1: i.e. the
interaction term b3 will be zero for all men. Similarly to above, we can test
whether there are indeed different returns to education or men and women by
looking at the estimated value b3 :
lm_w_interact <- lm(lwage ~ educ * female , data = wage1) # R expands to full interactions model
summary(lm_w_interact)

#OUT>
#OUT> Call:
#OUT> lm(formula = lwage ~ educ * female, data = wage1)
#OUT>
#OUT> Residuals:
#OUT> Min 1Q Median 3Q Max
#OUT> -2.02673 -0.27468 -0.03721 0.26221 1.34740
#OUT>
#OUT> Coefficients:
#OUT> Estimate Std. Error t value Pr(>|t|)
#OUT> (Intercept) 8.260e-01 1.181e-01 6.997 8.08e-12 ***
#OUT> educ 7.723e-02 8.988e-03 8.593 < 2e-16 ***
#OUT> female1 -3.601e-01 1.854e-01 -1.942 0.0527 .
#OUT> educ:female1 -6.408e-05 1.450e-02 -0.004 0.9965
#OUT> ---
#OUT> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#OUT>
#OUT> Residual standard error: 0.4459 on 522 degrees of freedom
#OUT> Multiple R-squared: 0.3002, Adjusted R-squared: 0.2962
#OUT> F-statistic: 74.65 on 3 and 522 DF, p-value: < 2.2e-16

We will in the next chapter learn that the estimate for b3 on the interaction
educ:female1 is difficult for us to distinguish from zero in a statistical sense;
Hence for now we conclude that there are no significantly different returns in
education for men and women in this data. This is easy to verify visually in
this plot, where we are unable to detect a difference in slopes in the right panel.
112 CHAPTER 5. CATEGORIAL VARIABLES

Impose Parallel Slopes Allow Different Slopes

3 3

2 2
log wage

1 1

0 0

0 5 10 15 0 5 10 15
educ educ

Figure 5.3: log wage vs educ. Right panel allows slopes to be different - turns
out they are not!

5.5 (Unobserved) Individual Heterogeneity

Finally, dummary variables are sometimes very important to account for spuri-
ous relationships in that data. Consider the following (artificial example):

1. Suppose we collected data on hourly wage data together with a the number
of hours worked for a set of individuals.
2. We plot want to investigate labour supply behaviour of those individuals,
hence we run regression hours_worked ~ wage.
3. We expect to get a positive coefficient on wage: the higher the wage, the
more hours worked.
4. You know that individuals are members of either group g=0 or g=1.
5.5. (UNOBSERVED) INDIVIDUAL HETEROGENEITY 113

10

5
hours

−5

0.0 0.5 1.0 1.5 2.0

wage

Here we observe a slightly negative relationship: higher wages are associated

with fewer hours worked? Maybe. But what is this, there is a group identifier
in this data! Let’s use this and include g as a dummy in the regression - suppose
g encodes male and female.
This is an artificial example; yet it shows that you can be severly misled if you
don’t account for group-specific effects in your data. The problem is particularly
accute if we don’t know group membership - we can then resort to advanced
methods that are beyond the scope of this course to estimate which group each
individual belongs to. If we do know group membership, however, it is good
practice to include a group dummy so as to control for group effects.
114 CHAPTER 5. CATEGORIAL VARIABLES

Controlling for Group

10
10

5
5 g
hours

hours

1
2

0
0

−5 −5

0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0
wage wage

Figure 5.4: Left and right panel exhibit the same data. The right panel controls
for group composition.
Chapter 6

Standard Errors

In this chapter we want to investigate uncertainty in regression estimates. We

want to understand what the precise meaning of the Std. Error column in a
typical regression table is telling us. In terms of a picture, we want to understand
better the meaning of the shaded area as in this one here:
In order to fully understand this, we need to go back and make sure we have a
good grasp of sampling. Let’s do this first.

6.1 Sampling

In class we were confronted with a jar of Tricolore Fusilli pasta as picture in

figure 6.2.1 We asked ourselves a question which, secretly, many of you had
asked themselves at one point in their lives, namely:

What is the proportion of green Fusilli in a pack of Tricolore Fusilli?

Well, it’s time to find out.

Let’s call the fusilly in this jar our study population, i.e. the set of units about
which we want to learn something. There are several approaches to address the
question of how big a proportion in the population the green Fusilli make up.
One obvious solution is to enumerate all Fusilli according to their color, and
compute their proportion in the entire population. It works perfectly well as a
solution, but is a long and arduous process, see figures 6.3 and 6.4.
1 This part is largely based on moderndive, to which I am giving full credit hereby. Thanks

for this great idea.

115
116 CHAPTER 6. STANDARD ERRORS

Log Wages vs Education

2
log hourly wage

10 15
years of education

Figure 6.1: Confidence bands around a regression line.

Figure 6.2: A glass jar filled with Fusilli pasta in three different colors.
6.1. SAMPLING 117

Figure 6.3: Manually separating Fusilli by their color is very costly in terms of
effort and cost.

Additionally, you may draw worried looks from the people around you, while
you are doing it. Maybe this is not the right way to approach this task?2

6.1.1 Taking One Sample From the Population

We started by randomly grabbing a handful of Fusilli from the jar and by letting
drop exactly N = 20 into a paper coffee cup, pictured in 6.5. We call N the
sample size. The count and corresponding proportions of each color in this first
sample are shown in the following table:

Color Count Proportion

Red 7 0.35
Green 5 0.25
White 8 0.4

So far, so good. We have our first estimate of the population proportion of green
Fusilli in the overall population: 0.25. Notice that taking a sample of N = 20
was much quicker and much less painful than performing the full count (i.e. the
census) of Fusilli performed above.
Then, we put my sample back into the jar, and we reshuffled the Fusilli. Had
we taken another sample, again of N = 20, would we again have gotten 7
Red, 5 Green, and 8 White, just as in the first sample? Maybe, but maybe not.
Suppose we had carried on for several times drawing samples of 20 and counting
the colors: Would we also have observed 5 green Fusilli? Definitely not. We
would have noted some degree of variability in the proportions computed from
our samples. The sample proportions in this case are an example of a sample
statistic.

Sampling Variation refers to the fact that if we randomly take sam-

118 CHAPTER 6. STANDARD ERRORS

Figure 6.4: Heaps of Fusilli pasta ready to be counted.

Figure 6.5: Taking one sample of 20 Fusilli from the jar.

6.2. TAKING ELEVEN SAMPLES FROM THE POPULATION 119

Figure 6.6: Taking eleven samples of 20 Fusilli each from the jar, and plotting
the histogram of obtained sample proportions of Green Fusilli.

6.2 Taking Eleven Samples From The Popula-

tion
We formed teams of two students in class who would each in turn take samples
from the jar (the population) of size N = 20, as before. Each team computed
the proportion of green Fusilli they had in their sample, and we wrote this data
down in a table on the board. Then, we drew a histogram which showed how
many samples had fallen into which bins.
We looked at the histogram in figure 6.6 and we noted several things:
1. The largest proportions where 0.3 green
2. The smallest proportion was 0.15 green.
3. Most samples found a proportion of 0.25 green fusilli.
4. We did think that this looked suspiciouly like a normal distribution.
We collected the sample data into a [Link]:
pasta_samples <- [Link](group = 1:11, replicate = 1:11, prop_green = c(0.3,0.25,0.25,0.3,0.15
pasta_samples

#OUT> group replicate prop_green

#OUT> 1 1 1 0.30
120 CHAPTER 6. STANDARD ERRORS

#OUT> 2 2 2 0.25
#OUT> 3 3 3 0.25
#OUT> 4 4 4 0.30
#OUT> 5 5 5 0.15
#OUT> 6 6 6 0.30
#OUT> 7 7 7 0.25
#OUT> 8 8 8 0.25
#OUT> 9 9 9 0.20
#OUT> 10 10 10 0.25
#OUT> 11 11 11 0.20
This produces an associated histogram which looks very much like the one we
draws onto the board:

Histogram of 11 Pasta Samples

5
4
Frequency

3
2
1
0

0.15 0.20 0.25 0.30

Proportion of Green Fusilli

6.2.1 Recap

Let’s recaptiulate what we just did. We wanted to know what proportion of

Fusilli in the glass jar in figure 6.2 are green. We acknowledged that an ex-
clusive count, or a census, is a costly and cumbersome exercise, which in most
circumstances we will try to avoid. In order to make some progress nonetheless,
we took a random sample from the full population in the jar: we randomly se-
lected 20 Fusilli, and looked at the proportion of green ones in there. We found
a proportion of 0.25.
After replacing the Fusilli from the first sample in the jar, we asked ourselves if,
6.3. HANDOVER TO MODERNDIVE 121

Figure 6.7: The Moderndive package used red and white balls instead of fusilli
pasta.

upon drawing a new sample of 20 Fusilli, we should expect to see the same out-
come - and we concluded: maybe, but maybe not. In short, we discovered some
random variation from sample to sample. We called this sampling variation.
The purpose of this little activity was three-fold:
1. To understand that random samples differ and that there is sampling
variation.
2. To understand that bigger samples will yield smaller sampling variation.
3. To illustrate that the sampling distribution of any statistic (i.e. not only
the sample proportion as in our case) computed from a random sample
converges to a normal distribution as the sample size increases.

The value of this exercise consisted in making you perform the sampling
activity yourself. We will now hand over to the brilliant moderndive
package, which will further develop this chapter.

6.3 Handover to Moderndive

The sampling activity in moderndive was performed with red and white balls
instead of green fusilli pasta. The rest is identical. We will now read sections
7.2 and 7.3 in their book, and return to this page later.

6.4 Inference in Theory

Imagine we were tasked by the Director of our school to provide him with
our best guess of the mean body height µ amongst all SciencesPo students in
order to assess which height the new desks should have. Of course, we are
122 CHAPTER 6. STANDARD ERRORS

Histogram of height

10
8
Frequency

6
4
2
0

140 150 160 170 180 190

height

Figure 6.8: Our ficitious sample of SciencesPo students’ body height. The small
ticks indicate the location of each measurement.

econometricians and don’t guess things: we estimate them! How would we go

about this task and estimate µ?

You may want to ask: Why bother with this estimation business at all, and
not just measure all students’ height, compute µ, and that’s it? That’s a good
question! In most cases, we cannot do this, either because we do not have access
to the entire population (think of computing the mean height of all Europeans!),
or it’s too costly to measure everyone, or it’s impractical. That’s why we take
samples from the wider population, to make inference. In our example, suppose
we’d randomly measure students coming out of the SciencesPo building at 27
Rue Saint Guillaume until we have 50 measurements on any given Monday.
Suppose further that we found a sample mean height x̄ = 168.5, and that the
sample standard deviation was s = 10. In short, we found the data summarized
in figure 6.8

What are we going to tell Monsieur le Directeur now, with those two numbers
and figure 6.8 in hand? Before we address this issue, we need to make a short
detour into test statistics.
6.4. INFERENCE IN THEORY 123

6.4.1 Test Statistics

We have encountered many statistics already: think of the sample mean, or the
standard deviation. Statistics are just functions of data. Test statistics are used
to perform statistical tests.
Many test statistics rely on some notion of standardizing the sample data so
that it becomes comparable to a theoretical distribution. We encountered this
idea already in section 3.2.4, where we talked about a standardized regression.
The most common standardization is the so-called z-score, which says that

x−µ
≡ z ∼ N (0, 1), (6.1)
σ

in other words, substracting the population mean from random variable x and
dividing by it’s population standard deviation yields a standard normally dis-
tributed random variable, commonly called z.
A very similar idea applies if we don’t know the population variance (which is
our case here!). The corresponding standardization gives rise to the t-statistic,
and it looks very similar to (6.1):

√ x̄ − µ
n ≡ T ∼ tn−1 (6.2)
s

Several things to note:

• We observe the same standardization as above: dividing by the sample
standard deviation s brings x̄ − µ to a unit free scale.
• We use x̄ and s instead
√ of x and σ
• We multiply by n because we expect x̄ − µ to be a small number: we
need to rescale it again to make it compatible with the tn−1 distribution.
• tn−1 is the Student’s T distribution with n − 1 degrees of freedom. We
don’t have n degrees of freedom because we already had to estimate one
statistic (x̄) in order to construct T .

6.4.2 Confidence Intervals

Back to our example now! We are clearly in need of some measure of confidence
about our sample statistic x̄ = 168.5 before we communicate our result. It seems
reasonable to inform the Director about x̄, but surely we also need to tell him
that there was considerable dispersion in the data: Some people were as short
as 143.98cm, while others were as tall as 189.41cm!
The way to proceed is to construct a confidence interval about the true popula-
tion mean µ, based on x̄, which will take this uncertainty into account. We will
use the t statistic from above. We want to have a symmetric interval around
124 CHAPTER 6. STANDARD ERRORS

tn−1

α
= 0.025 1 − α = 0.95 α
= 0.025
2 2

T
0 −c 0 c

Figure 6.9: Confidence Interval Construction. The blue area is called *coverage
region* which contains the true µ with probability 1 − α.

x̄ which contains the true value µ with probability 1 − α. One very popular
choice of α is 0.05, hence we cover µ with 95% probability. After computing
our statistic T as defind in (6.2), this interval is defined as follows:

Pr (−c ≤ T ≤ c) = 1 − α (6.3)

where c stands for critical value, which we need to choose. This is illustrated in
figure 6.9.

Given the symmetry of the t distribution it’s enough to find c at the upper tail:
the point above which α2 of all probability mass of the tdf distribution comes to
lie. In other words, if Tdf is the CDF of the t distribution with df degrees of
freedom, we find c as

α
Tdf (c) ≡ Pr (T < c) = 1 − = 0.975 (6.4)
2
c =Tdf−1 (Tdf (c)) = Tdf−1 (0.975) (6.5)

Here Tdf−1 stands for the quantile function, i.e. the inverse of the CDF. In our ex-
ample with df = 49, you can find thus that c = 2.01 by typing qt(0.975,df=49)
into your R session.3 Now we only have to expand the definition of the T statistic
from (6.2) inside (6.3) to obtain

3 You often will see c = 1.96, which comes from the fact that one relies on the t distribution

converging to the normal distribution with large n. Type qnorm(0.975) to confirm!

6.4. INFERENCE IN THEORY 125

0.95 = 1 − α = Pr (−c ≤ T ≤ c) (6.6)

( )
√ x̄ − µ
= Pr −2.01 ≤ n ≤ 2.01 (6.7)
s
( )
s s
= Pr x̄ − 2.01 √ ≤ µ ≤ x̄ + 2.01 √ (6.8)
n n

Finally, filling in our numbers for s etc, this implies that a 95% confidence
interval about the location of the true average height of all SciencesPo students,
µ, is given by:

CI = [165.658, 171.342] (6.9)

We would tell the director that with 95% probability, the true average height of
all students comes to lie within those two bounds.

6.4.3 Hypothesis Testing

Now know by now how the standard errors of an OLS estimate are computed,
and what they stand for. We can now briefly4 discuss a very common usage
of this information, in relation to which variables we should include in our
regression. There is a statistical proceedure called hypothesis testing which
helps us to make such decisions. In hypothesis testing, we have a baseline, or
null hypothesis H0 , which we want to confront with a competing alternative
hypthesis H1 . Continuing with our example of the mean height of SciencesPo
students (µ), one potential hypothesis could be

H0 :µ = 167 (6.10)
H1 :µ ̸= 167 (6.11)

Here we state that under the null hypthesis, µ = 167, and under the alternative,
it’s not equal to that value. This would be called a two-sided test, because it
tests deviations from H0 below as well as above. An alternative formulation
could use the one-sided test that

H0 :µ = 167 (6.12)
H1 :µ > 167. (6.13)
4 We will not go into great detail here. Please refer back to your statistics course from last

spring semester (chapters 8 and 9), or the short note I wrote while ago
126 CHAPTER 6. STANDARD ERRORS

which would mean: under the null hypothesis, the average of all ScPo students’
body height is 167cm. Under the alternative, it is larger. You can immediately
see that this is very similar to confidence interval construction.
Suppose as above that we found x̄ = 168.5, and that the sample standard
deviation is still s = 10. Would you regard this as strong or weak evidence
against H0 and in favor of H1 ?
You should now remember what you saw when you did launchApp("estimate").
Look again at this app and set the slider to a sample size of 50, just as in our
running example. You can see that the app draws one hundred (100) samples
for you, locates their sample mean on the x-axis, and estimates the red density.

The crucial thing to note here is that, given we are working with a ran-
dom sample from a population with a certain distribution of height, our
sample statistic x̄ is also a random variable. Every new set of ran-
domly drawn students would yield a different x̄, and all of them together
would follow the red density in the app. In reality we often only get to
draw one single sample, and we can use knowledge about the sampling
distribution to make inference.

Our task is now to decide if given that particular sampling distribution, given
our estimate x̄ and given an observed sample variance s2 , whether x̄ = 168.5
is far away from x̄ = 167, or not. The way to proceed is by computing a test
statistic, which is to be compared to a critical value: if the test statistic exceeds
that value, we reject H0 , otherwise we cannot. The critical value depends on
the sampling distribution, and the size of the test. We talk about this next.

6.4.4 Making Errors

There are two types of error one can make when deploying such a test:
1. We might reject H0 , when in fact it is true! Here, upon observing x̄ = 168.5
we might conclude that indeed µ > 167 and thus we’d reject. But we might
have gotten unlucky and by chance have obtained an unusually tall sample
of students. This is called type one error.
2. We might fail to reject H0 when in fact H1 is true. This is called the type
two error.
We design a test with a certain probability of type one error α in mind. In
other words, we choose with which probability α we are willing to make a type
one error. (Notice that the best tests also avoid making type two errors! The
number 1 − Pr(type 2 error) is called power, hence we prefer tests with high
power). A typical choice for α is 0.05, i.e. we are willing to make a type one
error with probability 5%. α is commonly called the level of significance or
the size of a test.
6.4. INFERENCE IN THEORY 127

f (x̄|µ = 167)

1 − α = 0.95
α = 0.05

x̄
0 µ = 167 x̄c = 171.74
√
t= n x̄−µ
s
0 tc = 1.676

Figure 6.10: Cookbook Testing Proceedure. Subscripts c indicate *critical

value*. There are two x-axis: one for values of x̄, and one for the corresponding
t statistic. The red area is the rejection area. If we observe a test statistic
such that t > tc , we feel reassured that our x̄ is *sufficiently far away* from
the hypothesized value µ, such that we feel comfortable with rejecting H0 . And
vice versa: If our test statistic falls below tc , we will not reject H0

6.4.5 Performing the Test

We can stick to the following cookbook procedure, which is illustrated in figure

6.10.
1. Set up hypothesis and significance level:
1. H0 : µ = 167
2. H1 : µ > 167
3. α = 0.05
2. Test Statistic and test distribution:
• We don’t know the true population variance σ 2 , hence we estimate
it via s2 from our sample.
• The corresponding test statistic is the t-statistic, which follows the
Student’s T distribution.
• That is, our statistic is T = s/ √ ∼ t49 , where 49 is equal to the
x̄−µ
n
degrees of freedom in this case.
3. Rejection Region: We perform a one-sided test. We said we are happy
with a 5% significance level, i.e. we are looking for the t value which
corresponds just to 1−0.05 = 0.95 mass under the pdf of the t distribution.
More precisely, we are looking for the 1 − 0.05 = 0.95 quantile of the t50
distribution.5 This implies a critical value c = 1.676, which you can verify
by typing qt(0.95,df=50) in R.
x̄−µ
√ = 168.5−167
4. Calculate our test statistic: s/ n
√
10/ 50
= 1.061
5. Decide: We find that 1.061 < 1.676. Hence, we cannot reject H0 , because
we only found weak evidence against it in our sample of data.
5 See the previous footnote for an explanation of this!
128 CHAPTER 6. STANDARD ERRORS

You can see from this that whether or not our test statistic is far way from the
critical value, or just below does not change our decision: it’s either accept or
reject. We never know if we narrowly rejected a H0 , or not. P-values are an
improvement over this stark dichotomy. The p-value is defined as the particular
level of significance α∗ , up to which all H0 ’s would be rejected. If this is a
very small number, we have overwhelming support to reject the null. If, on the
contrary, α∗ turns out to be rather large, we only found weak evidence against
H0 .

We define the p-value as the sum of rejection areas for a given test statistic T ∗ .
Notice that the symmetry of the t distribution implies that we would multiply
by two each of the two tail probabilities in the case of a two-sided test.

α∗ = Pr(t > |T ∗ |) (6.14)

6.5 Uncertainty in Regression Estimates

In the previous chapters we have seen how the OLS method can produce es-
timates about intercept and slope coefficients from data. You have seen this
method at work in R by using the lm function as well. It is now time to in-
troduce the notion that given that b0 , b1 and b2 are estimates of some unkown
population parameters, there is some degree of uncertainty about their values.
An other way to say this is that we want some indication about the precision of
those estimates. The underlying issue that the data we have at hand are usually
samples from a larger population.

How confident should we be about the estimated values b?

6.6 What is true? What are Statistical Models?

A statistical model is simply a set of assumptions about how some data have
been generated. As such, it models the data-generating process (DGP), as we
have it in mind. Once we define a DGP, we could simulate data from it and
see how this compares to the data we observe in the real world. Or, we could
change the parameters of the DGP so as to understand how the real world data
would change, could we (or some policy) change the corresponding parameters
in reality. Let us now consider one particular statistical model, which in fact
we have seen so many times already.
6.7. THE CLASSICAL REGRESSION MODEL (CRM) 129

6.7 The Classical Regression Model (CRM)

Let’s bring back our simple model (3.3) to explain this concept.

yi = β0 + β1 xi + εi (6.15)

The smallest set of assumptions used to define the classical regression model as
in (6.15) are the following:
1. The data are not linearly dependent: Each variable provides new infor-
mation for the outcome, and it cannot be replicated as a linear combina-
tion of other variables. We have seen this in section 4.2. In the particular
case of one regressor, as here, we require that x exhibit some variation in
the data, i.e. V ar(x) ̸= 0.
2. The mean of the residuals conditional on x should be zero, E[ε|x] = 0.
Notice that this also means that Cov(ε, x) = 0, i.e. that the errors and
our explanatory variable(s) should be uncorrelated. It is said that x should
be strictly exogenous to the model.
These assumptions are necessary to successfully (and correctly!) run an OLS
regression. They are often supplemented with an additional set of assumptions,
which help with certain aspects of the exposition, but are not strictly necessary:
3. The data are drawn from a random sample of size n: observation (xi , yi )
comes from the exact same distribution, and is independent of observation
(xj , yj ), for all i ̸= j.
4. The variance of the error term ε is the same for each value of x: V ar(ε|x) =
σ 2 . This property is called homoskedasticity.
5. The error is normally distributed, i.e. ε ∼ N (0, σ 2 )
Invoking assumption 5. in particular defines what is commonly called the normal
linear regression model.

6.7.1 b is not β!

Let’s talk about the small but important modifications we applied to model
(3.3) to end up at (6.15) above:
• β0 and β1 and intercept and slope parameters
• ε is the error term.
First, we assumed that (6.15) is the correct represenation of the DGP. With
that assumption in place, the values β0 and β1 are the true parameter values
which generated the data. Notice that β0 and β1 are potentially different from
b0 and b1 in (3.3) for a given sample of data - they could in practice be very
close to each other, but b0 and b1 are estimates of β0 and β1 . And, crucially,
those estimates are generated from a sample of data. Now, the fact that our
130 CHAPTER 6. STANDARD ERRORS

data {yi , xi }N
i=1 are a sample from a larger population, means that there will be
sampling variation in our estimates - exactly like in the case of the sample mean
estimating the population average as mentioned above. One particular sample
of data will generate one particular set of estimates b0 and b1 , whereas another
sample of data will generate estimates which will in general be different - by how
much those estimates differ across samples is the question in this chapter. In
general, the more observations we have the greater the precision of our estimates,
hence, the closer the estimates from different samples will lie together.

6.7.2 Violating the Assumptions of the CRM

It’s interesting to consider in which circumstances we might violate those as-

sumptions. Let’s give an example for each of them:
1. No Perfect Collinearity. We have seen that a perfect collinearity makes it
impossible to compute to OLS coefficients. Remember the example about
adding wtplus = wt + 1 to the mtcars dataset? Here it is:
library(dplyr)
mtcars %>%
mutate(wtplus = wt + 1) %>%
lm(mpg ~ wt + wtplus, data = .)

#OUT>
#OUT> Call:
#OUT> lm(formula = mpg ~ wt + wtplus, data = .)
#OUT>
#OUT> Coefficients:
#OUT> (Intercept) wt wtplus
#OUT> 37.285 -5.344 NA
That the coefficient on wtplus is NA is the result of the direct linear depen-
dence. (Notice that creating wtplus2 = (wt + 1)^2) would work, since
that is not linear!)
2. Conditional Mean of errors is zero, E[ε|x] = 0. Going back to our running
example in figure 6.1 about wages and education: Suppose that each in-
dividual i in our data something like innate ability, something we might
wish to measure with an IQ-test, however imperfecty. Let’s call it ai . It
seems reasonable to think that high ai will go together with high wages.
At the same time, people with high ai will find studying for exams and
school work much less burdensome than others, hence they might select
into obtaining more years of schooling. The problem? Well, there is no
ai in our regression equation - most of time we don’t have a good mea-
sure of it to start with. So it’s an unobserved variable, and as such, it is
part of the error term ε in our model. We will attribute to educ part of
the effect on wages that is actually caused by ability ai ! Sometimes we
6.8. STANDARD ERRORS IN THEORY 131

may be able to reason about whether our estimate on educ is too high or
too low, but we will never know it’s true value. We don’t get the ceteris
paribus effect (the true partial derivative of educ on lwage). Technically,
the assumption E[ε|x] = 0 implies that Cov(ε, x) = 0, so that’s the part
that is violated.
3. Data from Random Sample. One common concern here is that the obser-
vations in the data could have been selected in a particular fashion, which
would make it less representative of the underlying population. Suppose
we had ended up with individuals only from the richest neighborhood of
town; Our interpretation the impact of education on wages might not be
valid for other areas.
4. Homoskedasticity. For correct inference (below!), we want to know
whether the variance of ε varies with our regressor x or not. Here is a
typical example where it does:
Food Expenditure vs Log(income)
2000
1500
foodexp

1000
500

6.0 6.5 7.0 7.5 8.0 8.5

log(income)

As income increases, not all people increase their food consumption in an

equal way. So V ar(ε|x) will vary with the value of x, hence it won’t be
equal to the constant σ 2 .
5. If the distribution of ε is not normal, it is more cumbersome to derive
theoretical results about inference.

6.8 Standard Errors in Theory

The standard deviation of the OLS parameters is generally called standard error.
As such, it is just the square root of the parameter’s variance. Under assump-
132 CHAPTER 6. STANDARD ERRORS

tions 1. through 4. above we can define the formula for the variance of our
slope coefficient in the context of our single regressor model (6.15) as follows:

σ2
V ar(b1 |xi ) = ∑N (6.16)
i (xi − x̄)2

In pratice, we don’t know the theoretical variance of ε, i.e. σ 2 , but we form an

estimate about it from our sample of data. A widely used estimate uses the
already encountered SSR (sum of squared residuals), and is denoted s2 :

∑n ∑n 2
SSR i=1 (yi − b0 − b1 xi )2 e
2
s = = = i=1 i
n−p n−p n−p
where n − p are the degrees of freedom available in this estimation. p is the
number of parameters we wish to estimate (here: 1). So, the variance formula
would become

SSR
V ar(b1 |xi ) = ∑N (6.17)
(n − p) i (xi − x̄)2

We most of the time work directly with the standard error of a coefficient, hence
we define

√
√ SSR
SE(b1 ) = V ar(b1 |xi ) = ∑N (6.18)
(n − p) i (xi − x̄)2

You can clearly see that, as n increases, the denominator increases, and therefore
variance and standard error of the estimate will decrease.

6.8.1 Standard Errors in Practice

We would like to further make this point in an experiential way, i.e. we want
you to experience what is going on. We invite you to spend some time with the
following apps. In particular, make sure you have a thorough understanding of
launchApp("estimate").
library(ScPoEconometrics)
launchApp("estimate")
launchApp("sampling")
launchApp("standard_errors_simple")
launchApp("standard_errors_changeN")
6.8. STANDARD ERRORS IN THEORY 133

f (bk |βk = 0)

α
= 0.025 1 − α = 0.95 α
= 0.025
2 2

bk
0 −bck βk = 0 bck
bk −0
t= sk
tdown = −1.96 0 tup = 1.96

Figure 6.11: Testing whether coefficient bk is *statistically significantly different*

from zero. Now we have two red rejection areas. We relabel critical values with
a superscript here. If we observe a test statistic falling in either red region, we
reject, else we do not. Notice that the true value under H0 is βk = 0.

6.8.2 Testing Regression Coefficients

In Regression Analysis, we often want to test a very specific alternative hypoth-

esis: We want to have a quick way to tell us whether a certain variable xk is
relevant in our statistical model or not. In hypothesis testing language, that
would be

H0 :βk = 0 (6.19)
H1 :βk ̸= 0. (6.20)

Clearly, if in the true regression model we find βk = 0, this means that xk

has a zero partial effect on the outcome, hence it should be excluded from the
regression. Notice that we are interested in βk , not in bk , which is the estimator
that we compute from our sample (similarly to x̄, which estimates µ above).
As such, this is a two-sided test. We can again illustrate this in figure 6.11.
Notice how we now have two rejection areas.
The relevant test statistic for a regression coefficient is again the t distribution.
In fact, this particular test is so important that all statistical packages report
the t statistic corresponding to (6.20) automatically. Let’s look at an example:
#OUT>
#OUT> Call:
#OUT> lm(formula = mpg ~ wt + hp + drat, data = mtcars)
#OUT>
#OUT> Residuals:
#OUT> Min 1Q Median 3Q Max
134 CHAPTER 6. STANDARD ERRORS

#OUT> -3.3598 -1.8374 -0.5099 0.9681 5.7078

#OUT>
#OUT> Coefficients:
#OUT> Estimate Std. Error t value Pr(>|t|)
#OUT> (Intercept) 29.394934 6.156303 4.775 5.13e-05 ***
#OUT> wt -3.227954 0.796398 -4.053 0.000364 ***
#OUT> hp -0.032230 0.008925 -3.611 0.001178 **
#OUT> drat 1.615049 1.226983 1.316 0.198755
#OUT> ---
#OUT> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#OUT>
#OUT> Residual standard error: 2.561 on 28 degrees of freedom
#OUT> Multiple R-squared: 0.8369, Adjusted R-squared: 0.8194
#OUT> F-statistic: 47.88 on 3 and 28 DF, p-value: 3.768e-11
The column t value is just Estimate divided by Std. Error. That is, R reports
in the column t value the following number for us:

bk − 0
t value = (6.21)
sk

where sk is the estimated standard error as introduced in 6.8, and where we

test H0 : βk = 0. Notice that this particular t statistic √ is different from our
previous formulation in (6.2): we don’t have to scale by n! This is so because
R and other statistical software assumes the normal linear regression model (see
6.7). Normality of the regression error ε implies that the t statistic looks like in
(6.21).
We have to choose a critical value for this test. Many people automatically
choose the 0.975 quantile of the standard normal distribution, qnorm(0.975),
1.96 in this case. This is fine for sample sizes greater than 100, say. In this
regression, we only have 28 degrees of freedom, so we better choose the critical
value from the t distribution as above. We get tdown = −2.048 and tup = 2.048
as critical values. Let’s test whether the coefficient on wt is statistically different
from zero:

H0 :βwt = 0 (6.22)
H1 :βwt ̸= 0 (6.23)

We just take the t value entry, and see whether it lies above or below either
critical value: Indeed, we see that −4.053 < −2.048, and we are happy to reject
H0 .
On the other hand, when testing for statistical significance of drat that does
not seem to be the case:
6.9. WHAT’S IN MY MODEL? (AND WHAT IS NOT?) 135

H0 :βdrat = 0 (6.24)
H1 :βdrat ̸= 0 (6.25)

Here we find that 1.316 ∈ [−2.048, 2.048], hence it does not lie in any rejection
region, and we can not reject H0 . We would say that coefficient βdrat is not
statistically significant at the 5% level. As such, we should not include it in our
regression.

6.8.3 P-Values and Stars

R also reports two additional columns in its regression output. The so-called p-
value in column Pr(>|t|) and a column with stars. P-values are an improvement
over the dichotomy introduced in the standard reject/accept framework above.
We never know if we narrowly rejected a H0 , or not. The p-value is defined as
the particular level of significance α∗ , up to which all H0 ’s would be rejected. If
this is a very small number, we have overwhelming support to reject the null. If,
on the contrary, α∗ turns out to be rather large, we only found weak evidence
against H0 .

We define the p-value as the sum of rejection areas for a given test statistic T ∗ .
Notice that the symmetry of the t distribution implies that we multiply by two
each of the two tail probabilities:

α∗ = 2 Pr(t > |T ∗ |) (6.26)

The stars in the final column are a visualization of this information. They show
a quick summary of the magnitude of each p-value. Commonly, *** means
an extremely small reference significance level α∗ = 0 (almost zero), ** means
α∗ = 0.001, etc. In that case, up to a significance level of 0.1%, all H0 would be
rejected. You clearly see that all columns Std. Error, t value and Pr(>|t|)
give a different type of the same information.

6.9 What’s in my model? (And what is not?)

We want to revisit the underlying assumptions of the classical model outlined in

6.7. Right now we to talk a bit more about assumption number 2 of the above
definition in 6.7. It said this:
136 CHAPTER 6. STANDARD ERRORS

The mean of the residuals conditional on x should be zero, E[ε|x] = 0.

This means that Cov(ε, x) = 0, i.e. that the errors and our explanatory
variable(s) should be uncorrelated. We want x to be strictly exogenous
to the model.

Great. But what does this mean? How could x be correlated with something
we don’t even observe?! Good questions - let’s try with an example.
Imagine that we assume that

yi = β0 + β1 xi + εi (6.27)

represents the DGP of impact the sales price of houses (y) as a function of
number of bathrooms (x). We run OLS as

yi = b0 + b1 xi + ei
You find a positive impact of bathrooms on houses:
data(Housing, package="Ecdat")
hlm = lm(price ~ bathrms, data = Housing)
summary(hlm)

#OUT>
#OUT> Call:
#OUT> lm(formula = price ~ bathrms, data = Housing)
#OUT>
#OUT> Residuals:
#OUT> Min 1Q Median 3Q Max
#OUT> -77225 -15271 -2510 11704 102729
#OUT>
#OUT> Coefficients:
#OUT> Estimate Std. Error t value Pr(>|t|)
#OUT> (Intercept) 32794 2694 12.17 <2e-16 ***
#OUT> bathrms 27477 1952 14.08 <2e-16 ***
#OUT> ---
#OUT> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#OUT>
#OUT> Residual standard error: 22880 on 544 degrees of freedom
#OUT> Multiple R-squared: 0.267, Adjusted R-squared: 0.2657
#OUT> F-statistic: 198.2 on 1 and 544 DF, p-value: < 2.2e-16
In fact, from this you conclude that each additional bathroom increases the
sales price of a house by 27477 dollars. Let’s see if our assumption E[ε|x] = 0
is satisfied:
6.9. WHAT’S IN MY MODEL? (AND WHAT IS NOT?) 137

library(dplyr)
# add residuals to the data
Housing$resid <- resid(hlm)
Housing %>%
group_by(bathrms) %>%
summarise(mean_of_resid=mean(resid))

#OUT> # A tibble: 4 x 2
#OUT> bathrms mean_of_resid
#OUT> <dbl> <dbl>
#OUT> 1 1 -118.
#OUT> 2 2 955.
#OUT> 3 3 -11195.
#OUT> 4 4 32298.
Oh, that doesn’t look good. Even though the unconditional mean E[e] = 0 is
very close to zero (type mean(resid(hlm))!), this doesn’t seem to hold at all by
categories of x. This indicates that there is something in the error term e which
is correlated with bathrms. Going back to our discussion about ceteris paribus
in section 4.1, we stated that the interpretation of our OLS slope estimate is
that

Keeping everything else fixed at the current value, what is the impact of
x on y? Everything also includes things in ε (and, hence, e)!

It looks like our DGP in (6.27) is the wrong model. Suppose instead, that in
reality sales prices are generated like this:

y i = β0 + β1 x i + β2 z i + ε i (6.28)

This would now mean that by running our regression, informed by the wrong
DGP, what we estimate is in fact this:

yi = b0 + b1 xi + (b2 zi + ei ) = b0 + b1 xi + ui .

This is to say that by omitting variable z, we relegate it to a new error term,

here called ui = b2 zi + ei . Our assumption above states that all regressors need
to be uncorrelated with the error term - so, if Corr(x, z) ̸= 0, we have a problem.
Let’s take this idea to our running example.

6.9.1 Omitted Variable Bias

What we are discussing here is called Omitted Variable Bias. There is a variable
which we omitted from our regression, i.e. we forgot to include it. It is often
difficult to find out what that variable could be, and you can go a long way by
138 CHAPTER 6. STANDARD ERRORS

just reasoning about the data-generating process. In other words, do you think
it’s reasonable that price be determined by the number of bathrooms only? Or
could there be another variable, omitted from our model, that is important to
explain prices, and at the same time correlated with bathrms?

Let’s try with lotsize, i.e. the size of the area on which the house stands. Intu-
itively, larger lots should command a higher price; At the same time, however,
larger lots imply more space, hence, you can also have more bathrooms! Let’s
check this out:

#OUT>
#OUT> Call:
#OUT> lm(formula = price ~ bathrms + lotsize, data = Housing)
#OUT>
#OUT> Residuals:
#OUT> Min 1Q Median 3Q Max
#OUT> -60752 -12532 -1674 10514 92931
#OUT>
#OUT> Coefficients:
#OUT> Estimate Std. Error t value Pr(>|t|)
#OUT> (Intercept) 1.008e+04 2.810e+03 3.588 0.000364 ***
#OUT> bathrms 2.281e+04 1.703e+03 13.397 < 2e-16 ***
#OUT> lotsize 5.575e+00 3.944e-01 14.136 < 2e-16 ***
#OUT> ---
#OUT> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#OUT>
#OUT> Residual standard error: 19580 on 543 degrees of freedom
#OUT> Multiple R-squared: 0.4642, Adjusted R-squared: 0.4622
#OUT> F-statistic: 235.2 on 2 and 543 DF, p-value: < 2.2e-16

Here we see that the estimate for the effect of an additional bathroom decreased
from 27477 to 22811.5 by almost 5000 dollars! Well that’s the problem then.
We said above that one more bathroom is worth 27477 dollars - if nothing
else changes! But that doesn’t seem to hold, because we have seen that as we
increase bathrms from 1 to 2, the mean of the resulting residuals changes quite
a bit. So there is something in ε which does change, hence, our conclusion
that one more bathroom is worth 27477 dollars is in fact invalid!

The way in which bathrms and lotsize are correlated is important here, so
let’s investigate that:

This shows that lotsize and the number of bathrooms is indeed positively related.
Larger lot of the house, more bathrooms. This leads to a general result:
6.9. WHAT’S IN MY MODEL? (AND WHAT IS NOT?) 139

2e−04

bathrms
density

1
2
3
1e−04

0e+00

5000 10000
lotsize

Figure 6.12: Distribution of ‘lotsize‘ by ‘bathrms‘

Direction of Omitted Variable Bias

If there is an omitted variable z that is positively correlated with our
explanatory variable x, then our estimate of effect of x on y will be too
large (or, biased upwards). The correlation between x and z means that
we attribute part of the impact of z on y mistakenly to x! And, of course,
vice versa for negatively correlated omitted variables.
140 CHAPTER 6. STANDARD ERRORS
Chapter 7

Instrumental Variables

• Measurement error
• Omitted Variable Bias
• Reverse Causality / Simultaneity Bias
are all called endogeneity problems.

7.1 Simultaneity Bias

• Detroit has a large police force
• Detroit has a high crime rate
• Omaha has a small police force
• Omana has a small crime rate
Do large police forces cause high crime rates?
Absurd! Absurd? How could we use data to tell?
We have the problem that large police forces and high crime rates covary posi-
tively in the data, and for obvious reasons: Cities want to protect their citizens
and therefore respond to increased crime with increased police. Using math-
ematical symbols, we have the following system of linear equations, i.e. two
equations which are jointly determined:

crimeit = f (policeit )
policeit = g(crimeit )

We need a factor that is outside this circular system, affecting only the size
of the police force, but not the actual crime rate. Such a factor is called an
instrumental variable.

141
142 CHAPTER 7. INSTRUMENTAL VARIABLES
Chapter 8

Projects

This chapter contains several empirical projects.

8.1 Trade Exercise

• Trade exercise

143