7.1 Load today’s data

Start by loading today’s data. If you need to download the data, you can do so from these links, Save each to your data directory, and you should be set.

• hotdogs.csv
  • Data on nutritional information for various types of hotdogs
• rollerCoaster
  • Data on the tallest drop and highest speed of a set of roller coasters
• smokingCancer
  • Data on the smoking and mortality rates of various professions in Great Britain.
  • Smoking column gives a ratio with 100 meaning the profession smokes exactly the same amount as the general public and higher values meaning they smoke more
  • Mortality column gives a ratio with 100 meaning the profession has exactly the same lung cancer mortality rate as the general public and higher values meaning they have a higher mortality rate

# Set working directory
# Note: copy the path from your old script to match your directories
setwd("~/path/to/class/scripts")

# Load data
hotdogs       <- read.csv("../data/hotdogs.csv")
rollerCoaster <- read.csv("../data/Roller_coasters.csv")
smokingCancer <- read.csv("../data/smokingCancer.csv")

7.2 Compare Apples and Oranges

We have all heard that you can’t compare apples and oranges[*] Then again, some have: xkcd, and even a scientific journal . But, why not? In this section, we will see how to do exactly that using the standard deviation we just learned about as a yard stick.

Now, we can’t say that apples are better than oranges. But, we can say that this apple is a better apple than that orange is an orange. Said another way, we can identify which is more extreme, given the distribution of things like it.

Let’s be a bit more specific. Because the standard deviation is a measure of the width of the spread, it is directly comparable no matter what it is measuring[*] Recall that this is the formula for sd \[ \sqrt{ \frac{\sum\left(x - \bar{x}\right)^2}{n-1} } \] . Recall that, no matter what the units or thing being meausured is, approximately 95% of the values are within 2 standard deviations of the mean (as long as the distribution is unimodal and symmetrical). This means that we can measure how many standard deviations away from the mean a value is and that will tell us how “extreme” the value is – for that distribution

As an example, imagine that we decide that size (of fruit) matters. Thus, the best apple is the biggest one and the best orange is the biggest one. If we are given an orange that is 3 standard deviations bigger than the mean and an apple that is 2.5 standard deviations bigger than the mean, we can conclude that the orange is a better (or at least bigger) orange than the apple is an apple. This approach to measuring is so common that it has a name: z-scores.

# Calculate the mean
mean(rollerCoaster$Speed)

## [1] 65.11067

# How far is 85 from the mean?
85 - mean(rollerCoaster$Speed)

## [1] 19.88933

# Calulate the sd
sd(rollerCoaster$Speed)

## [1] 11.88235

# Calculate the z- score for an 85 mph coaster
( 85 - mean(rollerCoaster$Speed) ) / sd(rollerCoaster$Speed)

## [1] 1.673856

# Calculate the z- score for all coasters in the data set
( rollerCoaster$Speed - mean(rollerCoaster$Speed) ) / sd(rollerCoaster$Speed)

# Fast way to calculate z-scores
scale(rollerCoaster$Speed)

# Make histogram of scaled data
hist(scale(rollerCoaster$Speed),
     main = "How fast are roller coasters?",
     xlab = "Max Speed (sds from mean)",
     ylab = "Number of coasters")

# How extreme is a 100 mph coaster?
(100 - mean(rollerCoaster$Speed) ) / sd(rollerCoaster$Speed)

## [1] 2.936233

# How extreme is a 250 Calorie hotdog
(250 - mean(hotdogs$Calories) ) / sd(hotdogs$Calories)

## [1] 3.558322

7.3 Scatterplot

With z-scores out of the way, the next thing we might want to see is how two variables relate to each other. That is, can we visually see if high values of one variable are related to high values of another variable? For this, we obviously won’t be comparing completely unrelated data sets to each other. Rather, we need paired sets of data that can be linked together. In this section, we will see how to plot the relationships, in the next two sections, we will see how to put numerical summaries on the outcomes.

First, we will look at the structure of just such a dataset: the smokingCancer data.frame we loaded above.

# Look at the top of the data
head(smokingCancer)

As described above, this data gives information on the smoking and lung cancer mortality rates for 25 sets of professions in Britain. Each is scaled as a ratio with a value of 100 equating to the population average. We can plot the data using the default plot() function that will guess (usually correctly) at what kind of a plot we want. Here, that is a scatter plot, but we can get it the same way we got a box plot: using the formula and data = arguments to plot().

# Make a plot of smoking data
plot(Mortality ~ Smoking, data = smokingCancer)

With very little work, we have a servicable scatter plot. If we wanted to present this (in class, in a paper, etc.), we may want to clarify the labels a bit more, but for now, it will serve it’s purposes.

We now run into an interesting point that we must consider: which numeric variable should go on the x (bottom) axis and which should go on the y (left) axis? Because of the way people view plots, it is important to be careful with what may seem a small differnce. We always want to put the variable that we think “explains” the other variable (the one that “responds”) on the x-axis.

In this plot with a quick glance, we can see a trend in the plot (the next two sections will explore that trend). The plot implies that increased Smoking explains increased Mortality. However, it may be the case that jobs with higher Mortality explain why those people have higher Smoking rates. Always think carefully when deciding which order to put variables in. For all formulas we use in R, the order will be:

Response Variable ~ Explanatory variable

7.4 Correlation

Now that we have a picture of the relationship between Smoking and Mortality, we may want to ask how strong the relationship is. By this, we mean “how often are the large values of one variable matched with the large values of the other variable (and the small with the small)?” Now, talking about “large” and “small” values of a variable sounds rather familiar – it sounds like z-scores. Let’s convert our smoking data to z-scores and plot the scatter plot again. We can also add lines at the 0’s on each axis to show us where the means are using abline() (recall that the z-score for the mean is by definition 0).

# Plot the z-scores
plot(scale(Mortality) ~ scale(Smoking), data = smokingCancer)
abline(v = 0, h = 0)

It should become immediately obvious that most of the points are either in the top-right or the bottom-left quadrants. this suggests that jobs with above average smoking rates also have above average lung cancer mortality rates (and below average have below average). We can go a step further and ask how tightly that holds true.

The formula to calculate this is messy[*] \[ r = \frac{1}{n-1}\sum{\left(\frac{x - \bar{x}}{s_x}\right)\left(\frac{y - \bar{y}}{s_y}\right)} \] See why I hid this? , but the concept is simple. For each data point, multiply the two z-scores together. When both are above the mean, the product will be a positive value. When both are below the mean, the product will be a positive value (a negative times a negative is positive).

Thus, when we have a pattern like this, the values sum up to be a positive number. If they were in the top-left and bottom-right, the products would add up to a negative number. If they were a mix, the positives and negatives would cancel out and leave a value near 0.

By scaling this sum (again by n - 1), it turns out that it will always be between -1 (for a perfectly straight line from the top-left to bottom-right) and 1 (for a perfectly straight line from the top-right to bottom-left).

Instead of calculating this by hand, we will use the cor() function in R. Here, we will pass in two variables. The variables must be the same length (have the same number of values) because the order needs to match. That is, the first value for our “x” variable needs to be linked to our first “y” value (e.g., be from the same profession). Try it here:

# Test the correlation between Smoking and Mortality
cor(smokingCancer$Smoking,
    smokingCancer$Mortality)

## [1] 0.7162398

This value is the correlation coefficient, and is denoted by either an “r” or the Greek letter rho “ρ”.

## [1] 0.9096954

7.5 Linear Regression

To close this chapter, we will do formally what many readers probably did in their heads already: draw a line through the relationship.

First, make your plot of the smokingCancer data again (copy your command to generate it here).

# Make a plot of smoking data
plot(Mortality ~ Smoking, data = smokingCancer)

Now, imagine drawing a line through the data. You probably have an intuitive idea of where it goes, and are probably pretty close to accurate. But, how can we draw that line in R, and what does it mean?

For both questions, we will turn to the function lm() which creates a linear model of the data. Note that the arguments are exactly the same as the arguments to our call to plot(). Generally, the easiest way to write this line is to copy the plot() line, and just change plot to lm. Here, we will save the output of lm() because we are going to use it again[*] The output also contains a lot of other information, which we will use later in the semester. .

# Save a line for the plot
smokeLine <- lm(Mortality ~ Smoking,
                data = smokingCancer)
# Display the output
smokeLine

## 
## Call:
## lm(formula = Mortality ~ Smoking, data = smokingCancer)
## 
## Coefficients:
## (Intercept)      Smoking  
##      -2.885        1.088

The output gives us two values of coeffiecients: one for the intercept and one for the effect of Smoking. We can put these into an equation like this:

Mortality = -2.885 + 1.088*Smoking

The coefficient for Smoking gives us the slope of the line. That is, for every 1 unit increase in smoking, this line says that lung cancer mortality will increase by 1.088 units. The intercept just tells us where the line crosses the y axis. Imagine you put a 0 in for Smoking: all you would be left with on the right hand side would be the intercept.

This equation has a more general form of:

Response Variable = Β₀ + ExplanatoryVariable*Β₁

Where Β₀ (that is a captial Greek letter beta; this is pronounced “Beta Nought”) represents the intercept and Β₁ (“Beta 1”) is the coefficient for the explanatory variable, which tells us the slope of the line[*] This probably looks familiar to you as the y = mx + b you learned in algebra. Betas with subscripts are used in statistics because we are often going to have more than one explantory variable and using subscripts helps us keep track of things. .

Now, we can add this line to the plot with a function we already know: abline() works to plot the output of lm() as well. Simply give it our saved line with any parameters we want (e.g. color and thickness), and it will do the rest.

# Add a line to the plot
abline(smokeLine,
       col = "purple",
       lwd = 3)

7.6 Warnings and interpretation

There are a lot of caveats and warnings associated with investigating linear relationships between variables. Instead of putting them here in this already long chapter, they will be discussed separately after we have a chance to talk about where data like these come from.

7.7 Output

As with other chapters, click the notebook icon at the top of the screen to compile and HTML output. This will help keep everything together and help you identify any errors in your script.