Introduction

This is not intended to be an in-depth introduction to R. Nor is it intended to be an R Reference document. The purpose of this document is meant to get you started using R quickly and with minimal fuss as you proceed through your introductory course in data analysis. With this document, you should be able to quickly apply some of the statistical concepts you learn in the course to actual data.

Do not try to read through this entire document at once. Each section in this document is intended to accompany one or more of the lectures from your introductory course. My recommendation is to begin with the lecture from the course, and once you have studied the topic and done some problems, proceed to a read-through of the related section of these notes. Once you have done that, feel free to also watch the accompanying video where I discuss the section.

As you watch the video corresponding to the appropriate section of this document, I encourage you to press pause at regular intervals, in order to try out the techniques on your own. Do this with each section as we cover it throughout the course.

Over time, I think you will be impressed at the power of R. Even if you don’t see it right away, as your skills become increasingly advanced, you will begin to appreciate the endless possibilities of using R as your statistical software of choice.

Why R?

In the past, we have had to use commercial software such as SPSS or SAS. These are outstanding applications with decades of industry use and development behind them. However, they tend to be quite expensive. In addition, because of licensing issues, using them for our courses has always come with headaches for students and faculty alike.

Benefits to R

  • R is completely free!
  • R has been a staple of academic statistical software for years.
  • R is usually very easy to install.
  • There is a very easy-to-use, and also free development environment for using R called R Studio.

Downside to R

By far, the main downside to R, is that it has a bit of a steep learning curve. For this reason, in our course, we are not going to try to make you R experts. Nor are we going to spend several hours up front getting you up to speed with the language. Instead, we are going to take a bit of a “cookbook” approach in which I show you some code with which you can experiment, followed by some closely related requirements for your assignments.

Installing R and R Studio

We will be installing two applications. The first is R itself. After installing R, we will install a wonderful (and also free) tool called R Studio.

Installing R

So to begin, install R by visiting this page. You should see a page that looks something like this:

The red arrow shows you the most recent, stable installations of R. Download the appropriate one and install it on your machine.

Installing R Studio

Next, we will install R Studio.

Visit the R Studio download site. Then choose the free version of R Studio Desktop and install that application.

Beginning your work in R

Once you have installed both R and R Studio, you are all set to begin using R. To begin working, start R Studio on your computer.

Note: Do not start R. R Studio automatically starts R for you.

Once you have started R Studio, you should see something like the following:

You will probably see 4 separate panes which I have labeled in the above image. Initially, we will be typing our commands in the “Console Pane”. Very shortly, however, we will start spending the majority of our time in the “Script Pane”. As we progress, you will start to see the other two panes in action as well.

In the image shown here, I have typed a couple of simple mathematical commands in the console pane Once you have installed R Studio, try out a few for yourself.

Using R as a Calculator

A line of R code can work just like a calculator. The usual order of operations apply (parentheses, exponents, multiplication / division, addition / subtraction), etc. Here are some examples:

3+5
## [1] 8
(3*2)+5
## [1] 11
3.14*5^2  #In R, the ^ is used to represent exponent
## [1] 78.5

Creating a variable and storing some values in it

In our lecture, we discussed a series of heights in inches. Before long, we will learn how to retrieve lengthy sets of values from a text file, spreadsheet, or similar. For now, however, we will simply “hard code” that is, manually type out the values.

In the next line of code, I will create a variable called heights and will enter into that variable a series of values. We call this series of values a vector. Because we will use this term a lot, it’s worth memorizing.

We will start by placing just 3 values in our vector, but will subsequently enter all 25 heights as we discussed in our lecture.

heights <- c(58.2, 59.5, 60.7)

Note the following:

Okay. Once we’ve created our variable, we can have R do things for us such as tell us the sum of all of the values in the vector:

sum(heights)
## [1] 178.4

Tell us the mean of all of the values in the vector:

mean(heights)
## [1] 59.46667

Tell us the maximum and minimum value among all items in the vector:

max(heights)
## [1] 60.7
min(heights)
## [1] 58.2

We can have R sort the vector:

sort(heights)
## [1] 58.2 59.5 60.7

And many other things!

In fact, let’s now put in all of the remaining values from lecture into the vector, and start experimenting with it:

heights <- c(58.2, 59.5, 60.7, 60.9, 61.9, 61.9, 62.2, 62.3, 62.4, 62.9, 63.1, 63.4, 63.9, 64.0, 64.1, 64.5, 64.8, 65.2, 65.7, 66.2, 66.7, 67.1, 67.8, 68.9, 69.6)

We can print out all of the values by issuing the command print(). Inside the parentheses, we tell R to print out the values stored inside the variable ‘heights’:

print(heights)
##  [1] 58.2 59.5 60.7 60.9 61.9 61.9 62.2 62.3 62.4 62.9 63.1 63.4 63.9 64.0 64.1
## [16] 64.5 64.8 65.2 65.7 66.2 66.7 67.1 67.8 68.9 69.6

Let’s find the mean:

mean(heights)
## [1] 63.916

Median:

median(heights)
## [1] 63.9

Now let’s look at the survival times of “Disease X” referred to in lecture. In the example shown below, I create a vector called survivalTimes.

Important note: I know I said earlier that we should always name our variables with lower case letters. However, if our variable has more than one word in it such as survivalTimes, programmers usually capitalize the first letter of each subsequent word. This is called “camel case notation”. For example, if you had a variable with three words it might look like this: someRandomVariable.

In any case, here is our vector holding all of the survival times of Disease X as discussed in lecture.

survivalTimes <- c(.6, 1.2, 1.6, 1.9, 1.5, 2.1, 2.3, 2.3, 2.5, 2.8, 2.9, 3.3, 3.4, 3.6, 3.7, 3.8, 3.9, 4.1, 4.2, 4.5, 4.7, 4.9, 5.3, 5.6, 6.1)

Let’s ask R to tell us the 5-Number summary. This is easy to do using a function called fivenum(). We simply place our vector inside the parentheses, and R will give us the 5-number summary:

fivenum(survivalTimes)
## [1] 0.6 2.3 3.4 4.2 6.1

The values shown are the familiar Min, Q1, Median, Q3, Max discussed in lecture.

That being said, some of you may notice that the 5-number summary yielded by R’s fivenum() function looks different from the one demonstrated in lecture. This is because of the way Q1 and Q3 are calculated in R. They do things just a little bit differently from the way we did it manually. The good news is that moving forward, you can let R do all the work, so you don’t have to worry about this distinction.

Modifying an Existing Variable

Suppose you have an existing vector and you wish to modify it in some way. For example, suppose you wish to remove a value from it. There are easy ways to do this and more challenging ways. For now, we will learn the easy way. The disadvantage is that it is not very elegant.

To modify an existing variable, we will simply create a new variable, but give it a similar name.

For example:

golfScores <- c(78, 95, -4, 96)

However, we know that one of these values, the -4 must be incorrect, since in golf it is not possible to have a negative score. To fix this, I will simply create a new variable called, say golfScores2 without the bad value:

golfScores2 <- c(78, 95, 96)

Another option would be to simply recreate the original vector with the same variable name, but without the flawed -4 value:

golfScores <- c(78, 95, 96) In this particular case, I think this is a reasonable thing to do. Sometimes, however, we wish to keep the original vector as well as the updated one. In that case, we should create a new vector with it’s own variable name. We will now have two vectors, the original one, and the updated one. If, for whatever reason, we wished to revert back to our original vector, it would still be available to us under its own variable name.

Graphing in R

R can make extremely high-quality graphs. We will not get too elaborate with things for now, however. One nice thing about R is that even though graphing can get very advanced, R does make it very easy to create basic graphs as well.

Boxplots

For example, we can create a boxplot with the function boxplot():

boxplot(survivalTimes)

Histograms

We can create a histogram with the function hist():

hist(survivalTimes)

We will learn about various additional plots later in this document and as we proceed through the course.

Normal Quantile Plots

Note: When the time comes to draw this plot, I will refer you to this section of the document. If you reach this point in your reading and we have not yet discussed it in lecture, you may skip over this particular section. You do not need to concern yourselves with it until we cover it in class.

As discussed in lecture the normal quantile plot is intended to help support our theory that the variable under analysis is normally distributed. The normal quantile plot, also known as a “Q-Q Plot” is created by converting all of the observations to z-scores. A scatterplot is then created with z-scores on the x-axis and the actual observations on the y-axis. We should see a somewhat straight line beginning in the lower-left quadrant and sloping toward the upper right quadrant.

If the graph looks different, e.g. it is curved or the data clearly do not fall into the right-leaning relatively straight line pattern, then the variable is not normally distributed, and we should not proceed with normal distribution calculations.

It is also important to note that a seemingly straight line only supports normality, it does not guarantee it! Still, it’s a good bit of support to have – especially if we are about to embark on any kind of important analytics with the data.

Creating a QQ Plot

It is actually quite easy. Suppose we have a vector that is holding a bunch of race times such as:

raceTimes <- c(10.5, 13.2, 27.5, 15.6, 16, 28.7, 18.7, 4.9, 9.5, 11.4, 24.8, 17.9, 18.2, 15.9, 10.6, 29.3, 19, -0.7, 20.6, 11.2)

We can display a normal quantile plot with the simple function qqnorm() as follows:

qqnorm(raceTimes, 
       main="Normal Quantile Plot",
       xlab="Z-Score", ylab="Race Times")

We can see that the dots do indeed appear to follow a fairly straight line as discussed. This line does not have to be perfect, and I would be fairly comfortable doing normal distribution calculations moving forward.

By comparison, here is an example in which the normal quantile plot is clearly not a straight line:

General Graphing Tips

Improving Our Graphs

As mentioned above, we will start with very simple graphs, but. However, with experience, we can do quite a bit to finesse the appearance and complexity of our graphs.

We will start by making sure that we always include a title in our graphs. We do this by adding a parameter called main to our graphing function.

Let’s demonstrate by adding a title to our histogram:

hist(survivalTimes, main="Survival Times of People with Disease X")

A couple of important things to note:

  • Separate the argument main from the variable survivalTimes with a comma. We always separate each parameter from the next with a comma.
  • We must place the text of the title inside quotation marks

We can also adjust the labels on the x-axis and y-axis. For example, note how our x-axis simply is labeled with the name of the variable. This is not always ideal.

Instead, let’s relabel the x-axis to something more clear, such as: Survival Time in Years. To do this, we add another parameter to the hist() function called xlab:

hist(survivalTimes, main="Survival Times of People with Disease X", xlab="Survival Time in Years")

Reminder: As before, note how we separate each new argument from the one that preceded it with a comma. In the above example, when we decided to add an additional parameter (xlab) for the x label, we placed a comma before it.

Now let’s add a label on the y-axis. Not surprisingly, this argument is called ylab:

hist(survivalTimes, main="Survival Times of People with Disease X", xlab="Survival Time in Years", ylab="Number of People")

And again, note that we had to include an additional comma between the xlab parameter, and the subsequent ylab parameter. Any time we add an additional parameter to a function, we must separate it from the previous parameter with a comma.

TIP: Copy/Pasting Charts from R into Other Documents

You will frequently encounter situations where you want to paste a chart that you have created in R into other documents. For example, you will at various times need to paste your graphs into your assignment document. Fortunately, doing so is very simple. Here are the steps:

  1. Create your plot/chart
  2. In the output pane in R studio (the pane where your chart is visible), click on ‘Export’
  3. Choose ‘Copy to Clipboard’.

The chart will be copied to your clipboard and you can use ‘Paste’ (Control-V in Windows / Command-V on Mac) to paste the chart into your document.

Scatterplots

Recall that a scatterplot is used to look at the relationship between two variables. As an example, let’s work with our the dataset that examines the relationship between the number of beers consumed, and the resulting blood alcohol concentration.

Here is the original dataset:

Let’s begin by creating two vectors, one representing the number of beers consumed, and the other representing the blood alcohol concentration (BAC):

beers <- c(5,2,9,7,3,3,4,5,8,3,5,5,6,7,1,4)
bac <- c(.1, .03, .19, .095, .07, .02, .07, .085, .12, .04, .06, .05, .1, .09, .01, .05)

One important thing to note is that they must match up. For example, in our dataset, the first observation shows that for 5 beers, the BAC was 0.1. For the the second observation, 2 beers corresponded to a BAC of 0.03, and so on.

When creating our vectors in R, note how the first value in the beers vector (5 beers) corresponds to the first value in the bac vector (0.1). The second value in the beers vector (2 beers) corresponds to the second value in the bac vector (0.03), and so on.

Now let’s try to create a scatterplot using the plot() function. Before doing so, however, be sure and note the following:

When creating a scatterplot using the plot() function, the first argument to plot() will be printed on the x-axis, and the second argument will be printed on the y-axis. In other words, whichever variable you want to make the explanatory variable should be typed first, and whichever variable you want to have as the response variable should be typed second. Also be sure to include a comma between the two arguments. (In R, arguments are always separated by commas).

plot(beers, bac)

So we have a scatterplot, but hopefully you agree that it is not a particularly good plot. Let’s improve on it by adding labels. At the very least, a chart should always have a title, and nearly always, should have labels for both the x-axis and y-axis:

plot(beers, bac, 
     main="Effect of Beer Consumption on Blood Alcohol Level",
     xlab="# Beers Consumed", ylab="BAC")

Much better!

As you can probably deduce:

  • The parameter main allows you to specify text that will show up as the title of the chart
  • The parameter xlab allows you to specify text that will show up near the x-axis
  • The parameter ylab allows you to specify text that will show up near the y-axis

Okay, now let’s add a regression line. We will do this using a function called abline(). This function will require some discussion down the road, but for now, please don’t sweat it too much – just focus on the following:

  • The argument to the abline()function is another function called lm().
  • The argument to the lm() function should be in the form of the response variable, followed by a tilde, followed by the explanatory variable. Note that this is the opposite of the order we used for plot(). There is a good reason for this, but not one that we will delve into now.

Ultimately it will look as follows: abline(lm(bac~beers))

Let’s now put everything together

beers <- c(5,2,9,7,3,3,4,5,8,3,5,5,6,7,1,4)
bac <- c(.1, .03, .19, .095, .07, .02, .07, .085, .12, .04, .06, .05, .1, .09, .01, .05)

plot(beers, bac, 
     main="Effect of Beer Consumption on Blood Alcohol Level",
     xlab="# Beers Consumed", ylab="BAC")

abline(lm(bac~beers))

Creating Scripts in R

Imagine we wish to execute a long sequence of R commands. For example, suppose you wanted to create a vector:

scores <- c(82, 95, 87, 63, 92, 98, 54, 82, 79, 93, 95, 99)

followed by a histogram of those scores: hist(scores)

followed by a series of descriptive statistics:

mean(scores)

sd(scores)

median(scores)

fivenumber(scores)

and so forth.

We could,of course, type these out one command at a time in the console pane. However, suppose we needed to execute the same series of commands again at a later date. Or suppose we have an updated value in our scores vector (e.g. an additional observation, removal of an observation, correction of an erroneous value, etc. etc.)

We would have to retype these commands all over again. It might not be a big deal in this particular case, but in data analysis, we may often find ourselves issuing 5, 10, 20, or even 200 commands!

Wouldn’t it be nice to type out all of these commands once and then, with one command, tell R to run all of those lines of code?

Enter the “Script Pane”.

This is the pane that is probably visible at the top left of R Studio when you first open the application:

The script pane works just like the console. The difference is that whenever we type a command, we do NOT automatically see the results of that command executed. Instead, when we have typed out all of our commands, we can have R execute all of them at once (i.e. one at a time in order).

In the image shown here, you can see that I’ve typed out all of the commands in the script pane.

Note a couple of important points:

Try It! Type the following lines of code in the Script Pane, and save the file. Note that R will add a .R extension to the file name. Then click on the Source button to execute the code.

numbers <- c(4,3,5, 9, 10, 22)

mean(numbers)

Hmm, you should not see any errors, but you probably didn’t see any output show up in the Console Pane either.

The reason is that when we execute a script, we do NOT automatically see the output of every command. This is by design since if there were hundreds of lines of script, we could quickly get overwhelmed by output.

If we want to see the output of a command we must wrap the function inside a print() function. For example, to print out the mean of the numbers vectors we would have to type:

print( mean(numbers) )

Try it. Modify the script so that it looks as follows:

numbers <- c(4,3,5, 9, 10, 22)

print(mean(numbers))

Now let’s pretend that we have an additional value in our vector. Append the number 33 to the vector and click on the Source button to run the script again:

numbers <- c(4,3,5, 9, 10, 22,33)

print(mean(numbers))

Hopefully this went well.

Now let’s add calculation of the median to our script:

numbers <- c(4,3,5, 9, 10, 22,33)

print(mean(numbers))

print(median(numbers))

Continue with your practice and experimentation by adding the standard deviation (sd()), and five number summary (fivenum()).

Let’s finish with an example in which we create the scatterplot of the manatees vs. powerboat licenses using the Script window. In the script window, paste the following code:

boats <- c(447, 460, 481, 498, 513, 512, 526, 559, 585, 614, 645, 675, 711, 719)

manatees <- c(13, 21, 24, 16, 24, 20, 15, 34, 33, 33, 39, 43, 50, 47)

plot(boats, manatees, main="Manatee Deaths Resulting from Powerboats", xlab="Powerboat Licenses (in 1000s)", ylab="Manatee Deaths")

abline(lm(manatees~boats))

Then click on the Source button and see what happens.

Experiment by, say, changing the title of the graph or the labels. Click on the Source button again to confirm that the changes have indeed been made.

Do make sure that you can do this, as creating and submitting R scripts will be an important part of several assignments.

Creating Linear Regression Models

We will now learn how to generate regression models using R.

  1. We will start by creating a scatterplot of two variables to see if there is does indeed appear to be a linear relationship.
  2. If it looks as though a straight line may fit through the dots, we will generate a linear regression model using R’s lm() function. This function says to create a “linear model” (hence the “lm”). We will store this model inside a variable. We can call this variable anything we like, but I often simply call it model.
  3. Once we have a model, we will evaluate that model in order to check assumptions. The main assumptions we will look at from our model will be the strength of r, and the residual plot. If r appears to be reasonably strong, and the residual plot appears to be reasonably random, then we will write out our regression model.

Note: Another important assumption that must be met in order to trust a regression model, is that the relationship must be causal. If you have not yet encountered this topic in lecture, it will certainly be covered shortly in your course.

Let’s begin with our old friend, the beer / BAC dataset:

beers <- c(5,2,9,7,3,3,4,5,8,3,5,5,6,7,1,4)
bac <- c(.1, .03, .19, .095, .07, .02, .07, .085, .12, .04, .06, .05, .1, .09, .01, .05)

plot(beers, bac, 
     main="Effect of Beer Consumption on Blood Alcohol Level",
     xlab="# Beers Consumed", ylab="BAC")

abline(lm(bac~beers))

As we can see, the relationship certainly appears to be linear. Let’s therefore generate a regression model, so that we can see the value of r and check a residual plot. If the required assumptions appear to be met, we will look again at our model to obtain the regression formula.

Generating a Model:
It turns out that we have already seen this! In R, the function lm() generates a linear model. The argument to the lm() function must include the explanatory variable, and the response variable exactly as follows: Response Variable~Explanatory Variable. Note the tilde between the two. Also note that the response variable must be indicated first.

Now the lm() function returns a whole lot of information. We must therefore store this information inside a variable. If this sounds confusing, focus on the example shown here:

model <- lm(bac~beers)

R will generate a complete regression model, and all of the information will be stored inside this variable called model. Note: As indicated earlier, we do not have to call the variable model, I just choose that name since it’s pretty easy to remember. However, if in a given script we are working with multiple models, we would need to choose a different name for each one.

We should then “print” the summary of our model so that we can see things like R2 (and therefore, r as well).

Here is the code so far:

beers <- c(5,2,9,7,3,3,4,5,8,3,5,5,6,7,1,4)
bac <- c(.1, .03, .19, .095, .07, .02, .07, .085, .12, .04, .06, .05, .1, .09, .01, .05)

plot(beers, bac, 
     main="Effect of Beer Consumption on Blood Alcohol Level",
     xlab="# Beers Consumed", ylab="BAC")

abline(lm(bac~beers))

model <- lm(bac~beers)
print(summary(model))
## 
## Call:
## lm(formula = bac ~ beers)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.027118 -0.017350  0.001773  0.008623  0.041027 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.012701   0.012638  -1.005    0.332    
## beers        0.017964   0.002402   7.480 2.97e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02044 on 14 degrees of freedom
## Multiple R-squared:  0.7998, Adjusted R-squared:  0.7855 
## F-statistic: 55.94 on 1 and 14 DF,  p-value: 2.969e-06

Note the text that says: Multiple R-squared. This is our R2 value. Since it is 0.7988, then the value for r is the square root of that which is 0.89. This is quite a strong value for r, so that requirement is certainly met.

Now let’s check a residual plot. There are several ways of doing this, but the easiest one (though a bit clunky), is to simply plot the model: plot(model)

The reason this is a little strange is that R will show several plots. The one we are interested in is the first one called Residuals vs Fitted. We can ignore the remaining plots for now.

In this case, there is no obvious pattern to the residuals, so that assumption has been met as well. Incidentally, we can ignore the lines that R draws on the plot such as the red line shown in this example.

plot(model)

So let’s summarize our assumptions:

  1. Is r strong? Answer: Yes, very strong.
  2. Is the relationship linear? Answer: Yes, it looks inear to the eye, and there is no pattern in the residual plot.
  3. Is the relationship causal? Answer: Do we believe that drinking beer causes a change in BAC? Since we are confident that the answer is yes, then we can say that the relationship is causal.

Now that we have checked out key assumptions, we are reasonably confident that we can trust our model!

Let’s complete the process by converting the awkward text from the model provided by R, into a formula that anybody can use:

To do so, we require the values for b0 and b1. Both of these can be found when we print the summary of the model. Recall we do this by printing: summary(model) Look for a heading called Coefficients, and check under the Estimate column. The value labeled Intercept is b0, and the value labeled beers is our b1, i.e. the coefficient for the number of beers.

Call: lm(formula = bac ~ beers)

Residuals: Min 1Q Median 3Q Max -0.027118 -0.017350 0.001773 0.008623 0.041027

Coefficients:

         ` Estimate   Std. Error t value  Pr(>|t|)`

(Intercept) -0.012701 0.012638 -1.005 0.332

beers 0.017964 0.002402 7.480 2.97e-06

The estimated coefficient for the Intercept is -0.012701, and the estimated coefficient for beers is 0.017964.

So, our final model is:

BAC^ = -0.013 + 0.018*Number of Beers

Important: When you are asked to provide a regression model, you MUST provide it in the format shown here. That is, do not simply paste the output of the R summary. You must show that you can interpret the summary information and present it into a formula that anybody can understand.

Summary of the Code Used to Generate a Linear Regression Model

Here is the entirety of the code:

beers <- c(5,2,9,7,3,3,4,5,8,3,5,5,6,7,1,4)
bac <- c(.1, .03, .19, .095, .07, .02, .07, .085, .12, .04, .06, .05, .1, .09, .01, .05)

plot(beers, bac, 
     main="Effect of Beer Consumption on Blood Alcohol Level",
     xlab="# Beers Consumed", ylab="BAC")

abline(lm(bac~beers))

model <- lm(bac~beers)

print(summary(model))

plot(model)

Again as a reminder:

  • print(summary(model)) outputs information about the model including the coefficients (i.e b0 and b1), R2
  • plot(model) shows various graphs of the model including the residual plot

Bootstrapping

Bootstrapping is a technique used for inference, i.e. inferring information about a population through the use of sample data. It is often used to estimate certain summary statistics such as means or standard deviations, and others. One key piece of data that can be obtained via bootstrapping is the confidence interval of some statistic. In this discussion, We will focus on using the bootstrapping technique to obtain confidence intervals from a sample mean.

The key to bootstrapping lies in the ability to take a single samples, and from it, simulate many, many, many additional samples, for example, 1000. For each of these 1000 samples, we obtain the statistic of interest (e.g. the mean). As a result, we now have 1000 observations which we can use to generate a confidence interval.

This process should, hopefully seem familiar. This is because what we have just described is the creation of a sampling distribution. Recall that in the real world, we only work with one sample. However, if the sample is random, free from bias, etc., we can generally assume that this sample is likely to at least be a reasonable approximation of the population.

Once we have our sample, we will resample it. That is, we will take a sample from our sample with one key modification: We will allow the computer to sample any observation more than once! This is called “sampling with replacement”.

Let’s use a very tiny dataset (sometimes called a “toy” data set) to illustrate. Suppose we want to figure out the average age of a certain population (e.g. part-time workers at McDonald’s). From that population, we randomly sample 5 individuals and end up with:

18, 22, 19, 35, 26

To bootstrap, we would have the computer randomly select 5 observations, but “with replacement”. In other words, the computer will select the same number of observations as our original sample, but we are allowing it to select any observation more than once.

Here is our first bootstrap:

22, 19, 22, 35, 26

Note that 22 was selected twice, and the 18 was not selected at all. The mean age from this bootstrap is 24.8

Here is our second bootstrap:

26, 18, 19, 19, 19

In this case 19 was selected three times, and 22 and 35 did not get selected at all. The mean age from this bootstrap is 20.2.

Because computers do not go on strike, we can repeat this process as much as we like. We are only limited by our computer’s memory and speed. The number of bootstraps can vary significantly, but it should be somewhat large. Typical sizes are 1,000 and even 10,000.

Bootstrapping is a (relatively) new technique in statistics having been developed by Bradley Efron, a professor at Stanford, in 1979. One of the main reasons it was not developed prior to that point is that it requires a fair amount of computing memory and processor speed. For example, imagine trying to take a sample of say, 200 observations and calculating the mean. That may seem reasonable, but now imagine having to repeat that step 50,000 times! With modern computing speeds and memory, this should be a trivial operation.

But does it work?

One of the most common reactions when learning about bootstrapping is the sense that it almost seems too good to be true. How can take just one sample and then keep resampling from that same set of observations and end up with meaningful data? The biggest keys lie in:

  1. The assumption that our sample a truly random, unbiased subset of our original population.
  2. The fact that we sample “with replacement”. Because some observations will be left out each time, while others may be omitted on any given bootstrap means that we are hopefully getting a reasonable approximation of the randomness that occurred when we obtained our sample from the original population.

This is a very simplified explanation of why the bootstrap works, but there are countless resources online and in statistical texts that can explain in more detail if interested.

The efficacy of the bootstrap has been borne out over time. Since it’s inception in 1979, it has been used many, many times in a wide variety of disciplines, and studies have shown that the sampling distributions generated by a bootstrapping tend to be very good approximations of sampling distributions from the original populations.

Bootstrapping in R

It should be noted that there are various different statistical approaches for applying the bootstrap to generate confidence intervals including: standard, percentile, backwards, bias-corrected, and others. In this discussion we will limit ourselves to the “standard” approach.

Below I have created a function that simulates the idea behind boostrapping. It is not a particularly robust function for this process – for example, it only works on calculating means. However, by keeping the function relatively simple, it enables us to discuss the general idea of what the bootstrapping process does.

Note: As you are reading this section, please note that I do not expect you to fully understand the R code in it. For example, while the function below is not complicated code, but it does incorporate a number of programming constructs that are beyond the scope of this course.

In this example, I have written a function called bstrMeans that generates a vector of size B. Every item in the vector represents the mean of one sample. So if we have a vector of size 10, it means that we have created 10 samples, and stored the mean of each of those samples in the vector.

Again, I do not expect you to understand every line of this code. I do expect you to be clear on what the code is expected to do, as described above and continued below.

bstrMeans <- function(x, B) {  #This function will generate "B" bootstrap samples
  set.seed(1234)
  vecMeans <- vector(length=B)      #create an empty vector of size "B"
  for (i in 1:B+1) {
    s <- sample(x,replace=TRUE)       #generate one sample
    cat('Current Sample: ',s,'\t')    #print the sample
    cat('Mean of the current sample: ',mean(s),'\n')
    vecMeans[i-1]<-mean(s)            #add the mean of the current sample to our vector
  }
  return(vecMeans)
}

Let’s test the above function with our toy data set made of the ages of 5 randomly selected part-time McDonald’s workers. We will run the bootstrap 20 times:

mcdonalds <- c(18, 22, 19, 35, 26)
vctMeans <- bstrMeans(mcdonalds, B=20)
## Current Sample:  35 22 26 35 18  Mean of the current sample:  27.2 
## Current Sample:  26 35 22 22 35  Mean of the current sample:  28 
## Current Sample:  35 35 26 35 19  Mean of the current sample:  30 
## Current Sample:  35 26 22 26 22  Mean of the current sample:  26.2 
## Current Sample:  19 35 35 19 18  Mean of the current sample:  25.2 
## Current Sample:  19 35 22 19 22  Mean of the current sample:  23.4 
## Current Sample:  26 18 19 18 26  Mean of the current sample:  21.4 
## Current Sample:  18 18 22 18 19  Mean of the current sample:  19 
## Current Sample:  22 19 18 19 18  Mean of the current sample:  19.2 
## Current Sample:  22 26 18 19 19  Mean of the current sample:  20.8 
## Current Sample:  22 26 22 35 35  Mean of the current sample:  28 
## Current Sample:  18 26 19 26 19  Mean of the current sample:  21.6 
## Current Sample:  35 35 18 35 35  Mean of the current sample:  31.6 
## Current Sample:  19 18 35 18 22  Mean of the current sample:  22.4 
## Current Sample:  26 26 26 18 22  Mean of the current sample:  23.6 
## Current Sample:  22 22 26 19 35  Mean of the current sample:  24.8 
## Current Sample:  35 19 22 35 26  Mean of the current sample:  27.4 
## Current Sample:  22 18 19 19 18  Mean of the current sample:  19.2 
## Current Sample:  26 18 26 18 22  Mean of the current sample:  22 
## Current Sample:  19 19 26 35 19  Mean of the current sample:  23.6

Here is the vector generated by the function:

print(vctMeans)
##  [1] 27.2 28.0 30.0 26.2 25.2 23.4 21.4 19.0 19.2 20.8 28.0 21.6 31.6 22.4 23.6
## [16] 24.8 27.4 19.2 22.0 23.6

And here is a histogram:

hist(vctMeans)

I am now going to modify the above function to allow for very, very large bootstrap quantities. Note: One small but important modification I am making is to remove the print statements inside the function that is shown just above.

The reason I remove these statements is that if we want, say, 500,000 bootstraps, the last thing we wish to see is half a million of lines of text flying off the screen!

bstrMeans <- function(x, B=5) { 
  set.seed(1234)
  vecMeans <- vector(length=B)        #create an empty vector
  for (i in 1:B+1) {
    s <- sample(x,replace=TRUE)       #generate one sample
    vecMeans[i-1]<-mean(s)            #add the mean of the current sample to our vector
  }
  return(vecMeans)
}

Now let’s test our function with our McDonald’s ages variable, but with 100,000 bootstrap samples. Note that I assign B a value of 100,000 when I invoke the bstrMeans function:

mcdonalds <- c(18, 22, 19, 35, 26)
vctMeans <- bstrMeans(mcdonalds,B=100000)
hist(vctMeans)

The above histogram is made up of 100,000 observations. Every one of those observations is the mean of one sample taken from our original tiny vector of numbers: 2, 5, 7, 9, 1

Let’s try it with a slightly larger vector, one containing a series of heights taken from a sample of 25 people:

heights <- c(58.2, 59.5, 60.7, 60.9, 61.9, 61.9, 62.2, 62.3, 62.4, 62.9, 63.1, 63.4, 63.9, 64.0, 64.1, 64.5, 64.8, 65.2, 65.7, 66.2, 66.7, 67.1, 67.8, 68.9, 69.6)
vctMeans <- bstrMeans(heights,B=100000)
hist(vctMeans)

Let us now construct a 95% confidence interval from this heights variable:

We will start by calculating the mean and standard deviation of variable vctMeans.

meanHeight <- mean(vctMeans)
sdHeight <- sd(vctMeans)

To calculate the margin of error, we multiply 1.96 (since we are using a C of 95%) by the standard deviation.

IMPORTANT NOTE: We do NOT divide by the square root of n in this case. Right now we are working with a complete set of data – not one single sample. So in this case, we are NOT working with the standard error. Instead, we are working with a standard deviation since we are working with a complete set of samples. So in this case, we do not divide by the square root of n. Therefore, our margin of error is simply: z * SD.

moe <- 1.96 * sdHeight        # Calculate the margin of error and store in the variable 'moe'
negCI <- meanHeight - moe     #low CI value = mean - moe
posCI <- meanHeight + moe     #high CI value = mean + moe

So our confidence interval statement is the mean:

print(meanHeight)
## [1] 63.91489

With 95% Confidence interval of:

print(negCI)
## [1] 62.83072
print(posCI)
## [1] 64.99906

i.e. 62.83 - 64.99

Here is a summary of the code:

bstrMeans <- function(x, B) { 
  set.seed(1234)                      #Allows us to replicate random events
  vecMeans <- vector(length=B)        #create an empty vector of size 'B'
  
  for (i in 1:B+1) { 
    s <- sample(x,replace=TRUE)       #generate one sample
    vecMeans[i-1]<-mean(s)            #add the mean of the current sample to our vector
  }
  return(vecMeans)
}


#Create a vector of heights: 
heights <- c(58.2, 59.5, 60.7, 60.9, 61.9, 61.9, 62.2, 62.3, 62.4, 62.9, 63.1, 
63.4, 63.9, 64.0, 64.1, 64.5, 64.8, 65.2, 65.7, 66.2, 66.7, 67.1, 67.8, 68.9, 69.6)

means <- bstrMeans(heights,B=10000)
hist(means)

bsMean <- mean(means)
bsSD <- sd(means)
moe <- 1.96*bsSD

lowCI <- bsMean - moe
highCI <- bsMean + moe

print('Here is the mean:')
print(bsMean)
print('Lower Bound CI:')
print(lowCI)
print('Higher Bound CI:')
print(highCI)

If you wanted to replicate this process with a different set of values, you would only need to:

  1. Remove the heights vector and replace it with the vector you wish to work with. For example: ages <- c(33,44,55,66,77)
  2. Replace the vector you passed to the bstrMeans function. In the line means <- bstrMeans(heights, B=10000), replace heights with whichever vector you want to use. For example, if you had a vector of, say, people’s ages called ages, you might type: means <- bstrMeans(ages, B=10000)

Try it! Save this code into an R script. Then experiment by first trying it out the original code as shown in the summary above. Then replace the heights vector with a different vector, and try it again.

R’s boot() function

There is a function called boot() that specializes in doing bootstrap calculations. It has considerable more functionality than the function I created above. For example, in addition to bootstrapping calculations involving means, you can also do it on other statistics such as medians, regression coefficients, correlation coefficients – basically anything that you might want to calculate a statistic for. If you move on in your statistical training to a point where you start doing bootstrapping with any regularity, then you will very likely learn how to use this function or one just like it.

Understanding Data Sets

It is one thing to be given a ready-made vector and copy/paste it into R. For example: ages <- c(33,44,55)

In the real world, however, we pretty much always obtain our data from files. Files from which we obtain data can include things like text files, spreadsheets, database tables, web documents, documents created in other statistical software such as SAS or SPSS, and countless others.

In this section, we will discuss how to import data from text files, that is, files that can be opened up in any text editor. A text editor means that the file does not require specialized software to view the data (e.g. Excel to view a spreadsheet file, or Oracle to open up a database table). Instead, we can open up the file inside any text editor such as Windows Notepad or iOS TextEdit.

There are numerous reasons why data is stored in files including:

Let’s examine a couple of text files beginning with very simple ones, and progressing to more involved ones. As we do, download the file to your computer, and open it up in a text editor (e.g. Notepad, TextEdit, etc.) to see what the file looks like. Sometimes people ask me to recommend a good text editor. I do have one I think is terrific, and it’s also free. It is called Notepad++. A quick google search will take you to the download site.

Let’s begin with a very simple data set that has only one variable, “Heights”. This data set has 25 observations. Each observation represents the height of one person. You can download this file by right-clicking here and saving it to your computer.. You should see a page that looks something like this:

Note that only the first 8 observations are shown here.

Now let’s look at a dataset that has two variables.

This data set represents 10 observations. The two variables are:

You can download this file by right-clicking here and saving it to your computer.. You should see a page that looks something like this:

Again, be sure to right-click and save your own copy of the file. Do not simply try to view it in your browser window. Instead, you want to download it to your computer and open it in a text editor.

We should also note some other things about this file as they will turn out to be relevant:

Comma Separated Value (CSV) Files

In the previous example, we saw that in some text files, each item (i.e. column) was separated from the next by a tab. However, one of the most common ways in which text files are provided involves separating each item from the next with a comma. We call this a “comma separated value” or “CSV” file. These files are often given the extension .csv instead of the more familiar .txt. However, CSV files can be opened up in any text editor.

Let’s look at our Beer / Blood Alcohol example. Here is a version of the data set in which we have 16 observations, and two variables. The key thing to note is that in each row, each column is separated from the next by a comma. So the first individual had 5 beers and had a BAC of 0.1. The second individual had 2 beers and a BAC of 0.03 and so on.

The file can be downloaded here.

Here is what it looks like:

Beer / BAC Data Set With Additional Variables

Now let’s look at a version of the same data set that incorporates more information about each observation. In this improved data set, we are also provided with the individual’s gender, and their weight. So this data set has four variables: Gender, Weight, Beers, and BAC.

The file can be downloaded here.

Here is what it looks like:

A More Elaborate Data Set

Finally, let’s take a look at a very “real world” data set with many observations and many variables, such as one that might be downloaded from a government website.

This data set lists all 3000+ counties in the United States. For each county, the data set provides the name of the county, the state, the population in 2000, the population in 2010, the population in 2017, and other variables such as a poverty index, per-capita income, whether there is a smoking ban, and various other information.

The file can be downloaded here.

Here is what the first observations look like:

Opening up this data set in a text editor is certainly helpful to give us a quick glance at what is going on in the data set, but with over 3000 observations, and many variables, it is all but impossible to do any real analysis without using our statistical software.

Importing a Data Set into R

Now that we understand what data sets are, the question is how do we import them into R?

The answer is that there are many – perhaps too many – ways of doing so. One of the benefits and drawbacks to R is that there are often many different ways of accomplishing the same thing. This can be helpful as our skills advance, but it can be confusing to non-programmers learning R for the first time. For this reason, we will stick to very basic techniques. We will also limit ourselves to plain-text files such as txt and csv documents. Down the road, you may find yourself wanting to learn how to import Excel spreadsheets, database tables, SPSS documents, and similar as needed. However, CSV and TXT files are still arguably far and away the most common means by which data sets are made available.

The function to read in a CSV file is called read.csv(). There are two key arguments that must be provided:

  1. The name of the file
  2. Whether or not there is a header row

The command: df <- read.csv(file="beer_bac.csv", header=TRUE)

will look in the current folder (more on this in a moment) for a file named beer_bac.csv and will store the entire data set in a variable called df. You do not have to name the variable df. I am doing this to reflect that we will be working with a structure called a “data frame”. More on this in a bit. Also note the second argument that indicates that this data set has a header row. It is very important to specify either TRUE or FALSE as things can start behaving very strangely if this is wrong.

Location of the data set file

The argument given to the parameter called file must match a file name on your computer. Also, it is very important that you specify the proper path to your file. For example, if you have a Windows computer and your file is in a folder on your C drive called stats and a subfolder called datasets, your path would be: c:/stats/datasets/beer_bac.csv. So in this case your command would look like this:

df <- read.csv(file="c:/stats/datasets/beer_bac.csv", header=TRUE)

Also note the very important detail of using forward slashes in Windows. Mac users would use backslashes.

If the idea of paths names and folders is confusing to anyone reading this document, an easy way around this is to place the data set file in the same folder as your R script. If the data set file and your R script are in the same folder, then you do not need to specify any path. Simply type the file name and you will be good to go.

Reminder: Do not forget that R is case sensitive. So if your file happens to be called Beer_Bac.csv, then you must maintain the same case when you import your file.

Example: Importing a File into R

For this example, I will import the file bac_detailed.csv. This is the Beer / BAC data set that includes additional variables for things like gender, and weight. One key thing to note is that this file is stored on my computer in a folder on my C drive called temp.

df <- read.csv(file="c:/temp/bac_detailed.csv", header=TRUE)

R Studio (not regular R, but the development program R Studio) has a convenient function called View() that allows us to see a data set we have imported in spreadsheet format:

View(df)

One thing to keep in mind is that some datasets are huge! In this case, you may want to do something like:

head(df)
##   Gender Weight Beers   BAC
## 1 female    132     5 0.100
## 2 female    128     2 0.030
## 3 female    110     9 0.190
## 4   male    192     8 0.120
## 5   male    172     3 0.040
## 6 female    250     7 0.095

This command shows us the first 5 lines of the data set. This view is not as easy to interpret as the view given by the View() function, but it’s a great way to get a quick peek at a data set. I use it all the time.

The variable df now holds the entire data set! We can access any particular variable from the data set using the dollar sign operator, $. It works as follows:

numBeers <- df$Beers
bac <- df$BAC

We now have all of the items in the Beers column stored in the variable numBeers. We can print out that column, get its mean, SD, plot a histogram and so on:

print(numBeers)
##  [1] 5 2 9 8 3 7 3 5 3 5 4 6 5 7 1 4
mean(numBeers)
## [1] 4.8125

SD:

sd(numBeers)
## [1] 2.197536

Histogram:

hist(numBeers)

Let’s retrieve the blood alcohol concentration. Note: Recall that R is case sensitive. So when we want to retireve that column, the command df$bac will not work! Instead, we must type: df$BAC since that is what the column was named in the original text file.

So let’s retrieve all the values in this column and store in a variable called (lower case) bac:

bac <- df$BAC

We now have all of the values for the number of beers and the BAC. Let’s plot a scatterplot:

plot(bac,numBeers)

Let’s put all of this together:

df <- read.csv(file="c:/temp/bac_detailed.csv", header=TRUE)
numBeers <- df$Beers
bac <- df$BAC
plot(bac,numBeers)

From this point, you can now do anything you like with the bac and numBeers variables such as calculate mean, SD, print graphs such as histograms, and generate models such as linear regrssion models.

Let’s do an example using our powerboat licenses and manatee deaths data. We will:

  1. Import the dataset
  2. Retrieve the values for the number of powerboat licenses (“Powerboats”) and the number of manatee deaths (“Deaths”) and store those values into the variables boats and manatees respectively.
  3. Create a scatterplot
  4. Plot a regression line

We can find this data set in the file: boats_manatees.csv which can be downloaded here.. Be sure to save this file to your computer, and keep track of which folder you have placed it in.

Also, remember that it is nearly always a very good idea to first open the file in a text editor just to get some sense of what information is contained and what it all looks like. It can also be helpful since by scanning the document, you can determine things such as whether or not the file has a header line.

df <- read.csv(file="c:/temp/boats_manatees.csv", header=TRUE) 
#In the above line, be sure to modify your path to reflect where you have stored the file

View(df) #View the data set to get a sense of what it looks like
boats <- df$Powerboats
manatees <- df$Deaths
plot(boats, manatees,  
      main="Manatee Deaths Resulting from Powerboats", 
      xlab="Powerboat Licenses (in 1000s)", ylab="Manatee Deaths") 

abline(lm(manatees~boats))

model <- lm(manatees~boats)
print(summary(model))
## 
## Call:
## lm(formula = manatees ~ boats)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.2468 -2.0217  0.0217  2.3369  5.6328 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -41.4304     7.4122  -5.589 0.000118 ***
## boats         0.1249     0.0129   9.675 5.11e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.276 on 12 degrees of freedom
## Multiple R-squared:  0.8864, Adjusted R-squared:  0.8769 
## F-statistic: 93.61 on 1 and 12 DF,  p-value: 5.109e-07

Working with Packages

As with any widespread programming language or application (R qualifies as both), it soon becomes necessary to organize the many, many add-ons, specialized and/or esoteric uses, variants on existing tools, data sets, and so on into their own “packages”.

Most commonly, an R package is a collection of functions relating to a specific specific need or domain, along with help files to explain what the package does and how to use it. Packages may also include things like data sets intended to help demonstrate how the package works or to enable a new user to practice with it. However, none of this is set in stone. For example, some R packages are made up exclusively of data sets for people to use in their own work whether it be research, learning R, practicing building statistical models, and more.

Examples of packages include:

In order to work with a package, it must first be installed on your computer. This is exactly like installing any application on your computer.

Installing an R Package

Repositories: Installing an R package is (usually) quite easy. However, new R users are sometimes confused by the term “repositories”. A repository is simply a place where R packages are stored for people to download. The most popular one in the R world is called “CRAN” which is a group of web and FTP servers all over the world that store nearly all of the popular R packages in existence. However, a newer, or perhaps lesser known package may at times be found only on specific web sites or Git repositories. For the new R user, the key point for now is to bear in mind is that if you are asked for a “default repository”, or if you are asked which repository to try first, then CRAN is almost certainly your best bet.

To install a package, you would type the following: install.packages("PACKAGE_NAME")

Note that the name of the package should be in quotation marks. The install.packages() function defaults to looking for the package in the CRAN repository.

At this point, you will frequently see quite a lot of information scrolling down your R screen. You may even see warnings and errors, but these can usually be disregarded. Obviously, any overt error such as “package not available” would need to be investigated further. However, this should not be an issue for the majority of packages you are likely to encounter.

Loading a Package

As mentioned earlier, installing an R package is just like installing a program such as a web browser or any other application on your computer. However, as with computer applications, installing is not enough. WHen you wish to actually use the program, you must first load it. This is exactly like double clicking an icon for a program on your computer. However, in R, we must issue the library() command:

library("PACKAGE_NAME")

Technically, the library() function does not require quotation marks around the package name, but people sometimes get confused as to when quotes are and are not necessary, so for that reason I usually suggest that people include them.

If you are interested in more information about packages, there is a more detailed, but very easy to follow discussion here.

Use the following technique to install and load packages

Recall that we issue the install.packages() function only once on a given computer in order to install the package. We must then issue the library() function in every script where we plan to use that package.

One issue that comes up a lot, is that students frequently issue the install.packages() command in all of their scripts. This means that every single time they run their script, R will download and install the package anew. This can be tedious and time consuming. Imagine if you are sending your R script to a client (or a grader!) who already has those packages installed. When they run your script, they are going to have to sit back and wait while all of your packages reinstall themselves.

For this reason, please use the following command for every package you wish to use:

if  (!require("PACKAGE_NAME")){
  install.packages("PACKAGE_NAME")
  library("PACKAGE_NAME")
}

Obviously you should replace PACKAGE_NAME with your desired package on all three lines. Repeat those commands for every package you wish to use. What this code does is to first check to see if the package has been installed (this is the purpose of the require() function), and if the package is not installed, the code will first install, and then load the package.

You can repeat that block of code for every package you wish to use. For example, if you wanted to use the MASS and car packages, you would have:

if  (!require("MASS")){
  install.packages("MASS")
  library("MASS")
}

if  (!require("car")){
  install.packages("car")
  library("car")
}

That being said, if you are working in a situation in which you can confidently assume that the package has been installed (as is often the case), you can dispense with the somewhat awkward code above, and simply type:

library("MASS")
library("car")
library("ggplot2")
#etc

This technique of simply invoking library() when you wish to load a package is the preferred method. However, you should only do this once you’ve reached a point in your R work where you know that anyone who uses your code will also have these packages installed. Until that point, use the if (!require()) version from earlier.

IMPORTANT:

Be sure to place your code for installing and loading packages at the top of your R script, since a package must be loaded into your R script before any of the functions or data sets contained inside can be used.

Tip: A convenient command to see which packages are installed on your computer is: library(). If you issue this command in your console, you will see that you already have quite a few libraries installed on your computer. This is because the base installation of R comes bundled with some commonly used packages. In fact, some of them are so common, that they do not even need to be loaded.

Example

Let’s do an example. R’s generic plot() function, while powerful and versatile, does have its limitations. If you are planning on doing a lot of scatterplots, there is a function called (surprise, surprise), scatterplot() that has all kinds of useful features for people who do a lot of these including automatic regression lines, boxplots for each variable being plotted, interactive identification of observations, and more. This function is present in a package called car.

#Ensure the car package is installed and loaded
if  (!require("car")){
  install.packages("car")
  library("car")
}
## Loading required package: car
## Warning: package 'car' was built under R version 4.0.2
## Loading required package: carData
beers <- c(5,2,9,7,3,3,4,5,8,3,5,5,6,7,1,4)
bac <- c(.1, .03, .19, .095, .07, .02, .07, .085, .12, .04, .06, .05, .1, .09, .01, .05)

#Issue the scatterplot() function
#Note that the syntax is a little bit different from the plot() function,
#namely:  response_variable ~ explanatory_variable
scatterplot(bac~beers, 
     main="Effect of Beer Consumption on Blood Alcohol Level",
     xlab="# Beers Consumed", ylab="BAC")

Warnings

Recall that you may well encounter various warnings when you load a package. Typically, warnings are just that and can be safely ignored. However, if you are writing code intended for “real world” use, you would want to look at any warnings a little more closely, so as to ensure that they are not doing anything that might cause problems with your code down the road.

Getting help

Once you have loaded a package, you can get help for individual functions and data sets contained in the package. For example:

help(scatterplot)

will open up a document in the Help pane of R Studio with all kinds of information about that particular function. Don’t be intimidated by all of the details in the help function. With just a little time and additional experience, you will get better at interpreting what all of it means.

Multiple Regression

In this section, we are going to demonstrate some basic steps to enable you to get going with multiple linear regression in R. The objective here is not to go into detail about multiple regression, but rather, to understand the basics of how to create a model using R, how to evaluate the summary of the model to get a sense of its efficacy, and how to generate a model that can be interpreted by a relative layperson,

We will use the bac_detailed.csv data set. As always, it’s a good idea to study any information provided to us about the data set. If we do not have any, then we must, at the very least, open up the data set in a text editor to get a sense of what it looks like. Recall that we are looking for things like whether or not there is a header row, the number of variables, how many observations, and more. If and when you proceed with more advanced work in data analysis, you will learn to examine the data for other details such as whether or not there are missing values, erroneous values, inconsistent data entry, etc. etc.

Here is the file as seen in a text editor:

Because this is a small data set, it is quite easy to quickly get a sense of what is going on. Some things one might note include:

We begin, of course, by reading the file into into R:

df <- read.csv(file="c:/temp/bac_detailed.csv", header=TRUE)

Recall that the read.csv() function returns the entire data set as a structure called a “data frame”. This is basically a table – albeit one that can have more than two dimensions. This data frame is now stored inside the variable called df.

It’s always a good idea to take a quick look to make sure that we have done things correctly. Let’s print out the first few rows of the data frame:

head(df)
##   Gender Weight Beers   BAC
## 1 female    132     5 0.100
## 2 female    128     2 0.030
## 3 female    110     9 0.190
## 4   male    192     8 0.120
## 5   male    172     3 0.040
## 6 female    250     7 0.095

We can now see that the data set has indeed been properly read in from the file.

Just for fun and review, let’s do a few summary statistics:

print( mean(df$Weight))
## [1] 171.5625
print( sd(df$Weight ))
## [1] 48.96389
hist(df$Weight)

Let’s build a regression model using R’s lm() function. We will begin by doing a simple linear regression, i.e. having only one explanatory variable. In this case, the variable, Beers.

model <- lm(df$BAC~df$Beers)

Be sure to note that in this particular data set, the variables are called BAC and Beers, that is, make sure to respect the case of the letters. If we had typed, say,

model <- lm(df$bac~df$beers)

we would have gotten an error and our script would have immediately halted.

Now let’s take a look at our model:

summary(model)
## 
## Call:
## lm(formula = df$BAC ~ df$Beers)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.027118 -0.017350  0.001773  0.008623  0.041027 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.012701   0.012638  -1.005    0.332    
## df$Beers     0.017964   0.002402   7.480 2.97e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02044 on 14 degrees of freedom
## Multiple R-squared:  0.7998, Adjusted R-squared:  0.7855 
## F-statistic: 55.94 on 1 and 14 DF,  p-value: 2.969e-06

Here are some important details to note:

So our model is: BAC^ = -0.013 + 0.018 * Beers

However, there are a couple of additional very important details to note:

However, with just a little bit of thought, we can all hopefully recognize that the variable Weight surely also plays a significant role in predicting BAC. Let’s now include a second explanatory variable in our model. In doing so, we are not doing “multiple” linear regression since there is more than one explanatory variable in the model.

To include additional variables in a model, we simply add them to the first explanatory variable, separating each one from the next with a + sign:

model2 <- lm(df$BAC~df$Beers + df$Weight)

Note the + df$Weight that has been added.

Now let’s take a look at our newer model:

summary(model2)
## 
## Call:
## lm(formula = df$BAC ~ df$Beers + df$Weight)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0162968 -0.0067796  0.0003985  0.0085287  0.0155621 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.986e-02  1.043e-02   3.821  0.00212 ** 
## df$Beers     1.998e-02  1.263e-03  15.817 7.16e-10 ***
## df$Weight   -3.628e-04  5.668e-05  -6.401 2.34e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01041 on 13 degrees of freedom
## Multiple R-squared:  0.9518, Adjusted R-squared:  0.9444 
## F-statistic: 128.3 on 2 and 13 DF,  p-value: 2.756e-09

Note:

Had the value for R2 decreased, we would probably have elected to not include this additional variable in our model. However, because including Weight has improved our R2, we will keep it.

Therefore, our new and improved model is: BAC^ = 0.03986 + 0.01998 * Beers - 0.0003628 * Weight

Also, note the minus sign before Weight. This is because the coefficient for Weight is negative 0.003628.

Conclusion: We have barely scratched the surface of multiple regression, but hopefully, you have a sense of what it means to include additional variables in a model, and how to create and interpret models using R. For those of you who continue on with more advanced study in statistics, you will learn all about multiple regression and model building. Some topics you will cover will include:

I don’t know whether this seems exciting to you or dry – but rest assured, it can be an absolutely fascinating field! The world runs on data and model-building. In fact, this data accumulation, analysis, and model-building is the “money” behind the massively successful e-sites such as Google and Facebook. Even if you take issue with how those sites work, these techniques are also how we predict the weather, evaluate whether or not a bridge will survive an earthquake, decide whether or not you should take a certain medication for a chronic disease, predict what movies you might like to see on Netflix, suggest some products you might want to buy on Amazon, and more.

It’s a pretty amazing field to be a part of. I wish you luck.

-Yosef