Lab: Introduction to R and RStudio

Single-variable plots in lattice

For one quantitative variable, you might want to produce a histogram to show the distribution.

histogram(~salary, data=Salaries)

You could also view a smoothed version of the distribution with a density plot

densityplot(~salary, data=Salaries)

If you have categorical data, a barchart is more appropriate.

barchart(tally(~rank, data=Salaries))

Notice that all these plots use the same syntax: goal ( y ~ x , data = mydata )

We call this the formula syntax, and it is the most common way we will be writing R code in this class.

In this case, we ony had an x variable, so we put nothing before the tilde, ~. (Hint: to find the tilde button, look at the upper left of your computer keyboard. It usually shares a key with the backtick.)

An aside– RMarkdown

This lab (and all the rest of them) is written in RMarkdown. That means I type text, code, and formatting into a document with the file extension .Rmd. That file is available on the course site if you would like to look at it. Some students prefer to download the source code and work through it on their own machine or on the RStudio server. Once I have things typed into my .Rmd file, I click the “Knit” button to knit a formatted HMTL version of the lab. If you are doing this on your own machine and the knit button isn’t showing, check out Prof McNamara’s RMarkdown troubleshooting page. When I write the labs, I use the option eval=FALSE, which means that the output of the code doesn’t show up in my knitted document. For your homework, you need to use the option eval=TRUE so that the output shows up in your HTML. If you are working through this lab on your own computer, go to the top of the document and change this option. Then try re-knitting. See how all the plots and output show up in the document?

For homework, you need to submit the knitted HTML document. So, it is important to make sure your document has knitted correctly. Again, the RMarkdown troubleshooting page has some hints for how to do this. I usually just open it in my web browser to check it looks right.

Single-variable plots in ggplot2

ggplot2 can create all the same plot types as lattice. A histogram,

require(ggplot2)
ggplot(Salaries) + geom_histogram(aes(x=salary))

A density plot,

ggplot(Salaries) + geom_density(aes(x=salary))

Or, for categorical data, a barchart.

ggplot(Salaries) + geom_bar(aes(x=rank))

Notice that the syntax is different with ggplot. It’s now ggplot(mydata) + goal(aes(x=x, y=y)) where your goal is probably a geom (short for geometric object), like geom_histogram(), geom_density() or geom_bar(). There are actually other ways to write ggplot2 code, but we won’t think about those for now.

Multiple-variable plots in lattice

For two quantitative variables, a scatterplot is appropriate.

xyplot(salary~yrs.since.phd, data=Salaries)

For one quantitative and one categorical variable, you can make side-by-side boxplots (boxplot(), bwplot()),

bwplot(salary~sex, data=Salaries)

For two categorical variables, a mosaic plot is most appropriate

mosaicplot(sex~rank, data=Salaries, las=1)

We’re back to that goal ( y ~ x , data = mydata ) syntax.

Multiple-variable plots in ggplot2

A scatterplot,

ggplot(Salaries) + geom_point(aes(x=yrs.since.phd, y=salary)) 

Side-by-side boxplots,

ggplot(Salaries) + geom_boxplot(aes(x=sex, y=salary)) 

Again, we’re using ggplot(mydata) + goal(aes(x=x, y=y)) where your goal is probably a geom (short for geometric object), like geom_boxplot().Mosaic plots in ggplot2 are more complicated and require an additional package, so we’ll leave those for now.

Bells and Whistles

Most R plotting function can take many arguments that will modify their behavior. Read the documentation for more information.

help(xyplot)
?geom_histogram

Summary statistics

Again, there are several possible syntaxes for summary statistics. I mostly use the mosaic package syntax, which looks like this:

mean(~salary, data=Salaries)
sd(~salary, data=Salaries)

This is the same formula syntax as lattice plotting. goal ( y ~ x , data = mydata ), and your goal can be a summary statistic like mean() or sd(). Most summary statistics you would expect are in R. median(), max(), min(), IQR(), and many more all work. There is also a function that gives you a bunch of my favorite statistics, favstats()

favstats(~salary, data=Salaries)

For categorical data, you might want to tally()

tally(~rank, data=Salaries)

Linear models

The syntax for linear models is very similar to the syntax for lattice graphics and mosaic summary statistics. goal ( y ~ x , data = mydata ) We’ll use it for modeling the data we collected about ourselves.

In order to read in the data, we need to use an assignment operator so we can store the data in our environment. Without an assignment operator, the data would all print out into the Console but not be useable again.

We’ll try to be consistent and use = as our assignment operator, but there is another way to do it (of course!) that I sometimes use that looks like this: <-.

# install.packages("googlesheets")
require(googlesheets)
url = gs_url("https://docs.google.com/spreadsheets/d/1tHg_Oo88GIS8E00-lf286-4AzR-wht7ZddUMFM1s5s0/")
ds = gs_read_csv(url)

It’s not crucial that you understand the code above, because most of the data we use in this course is in packages.

Now we can start to do some modeling! We might want to start with some EDA to explore the relationship between distance from home and number of countries visited.

xyplot(numFBFriends ~ accountLength, data=ds)

Given how the data looks, it seems reasonable to choose a simple linear regression model. We can then fit it using R.

mod = lm(numFBFriends ~ accountLength, data=ds)

Again, we’re using the assignment operator to store the results of our function into a named object. This is because we might want to do things with our model output later. If you want to skip the assignment operator, try just running lm(numFBFriends ~ accountLength, data=ds) in your console to see what happens.

Now, we want to move on to assessing and using our model. Typically, this means we want to look at the model output (the fitted coefficients, etc). If we run summary() on our model object we can look at the output.

summary(mod)

The p-values is significant, which might lead us to assess that the model was effective. We can use the model for description. To write this model out in our mathematical notation,

$$\widehat{\verb#numFBFriends#} = -5.264 + 84.164\cdot \verb#accountLength#$$

1. What would this model predict for the number of friends if this person had an account for 5 years?
2. What does the model predict for you? What was the observed value of your number of friends? What is your residual?

To assess further, we can compare it to a “null” model, where we don’t use a predictor (Instead, we just use the mean as the model for every observation). We can compare the residuals from the null model and the simple linear model.

In order to visualize this, we need a special plotting function. Copy-paste this code into your Console without changing anything:

panel.regplot = function (x, y, mod, ...) {
  panel.xyplot(x, y, type = c("p"), lwd=3, pch=19, ...)
  panel.arrows(x0=x, y0=y, x1=x, y1=fitted.values(mod), code=2, length=0.1, col="darkgray")
  panel.abline(mod)
}

Now, we can run two models and visualize them to compare! The first is the average, and the second is the least squares regression line. We can think of the latter as a the null model where $\beta_1 = 0$ and $\hat{y} = \bar{y}$. Which model do you think is better?

mod.mean = lm(numFBFriends ~ 1, data=ds)
mod.line = lm(numFBFriends ~ accountLength, data=ds)
xyplot(numFBFriends ~ accountLength, data=ds, panel=panel.regplot, mod = mod.mean)
xyplot(numFBFriends ~ accountLength, data=ds, panel=panel.regplot, mod = mod.line)

Your turn:

1. Think about relationships between variables in this data. Identify two bivariate relationships that you expect to be strongest, and two that you expect to be weakest.
2. Explore these relationships using plots. Were you right? What criteria did you use to assess?
3. Use lm() to fit a linear regression model to one pair of variables and examine the output using summary(). What conclusions can you draw?