# 8 Graphics

Now that we know a bit about how data can be stored and manipulated in data frames, we can begin to analyze data, so in this Chapter we will take a closer look at the aspect of data analysis known as data visualization. We’ll become acquainted with the grammar of graphics, a general approach to visualization, and then learn how this approach is implemented in the ggplot2 with which we have worked previously. Finally, we’ll engage in a case study.

## 8.1 The Grammar of Graphics

A graph begins with data, and the data we work with will be tidy data that comes in a data frame. Leland Wilkinson’s Grammar of Graphics (see ) posits that most quantitative graphics constructed from a data frame can be understood in terms of a few basic elements. In our quite elementary introduction to the Grammar, the elements to which we will pay the most attention are as follows:

• Glyphs: the basic units of a graph. Glyphs represent cases in the data frame. Each glyph corresponds to one or more cases, but no two glyphs correspond to the same case.
• Aesthetics: perceptual properties of glyphs that are not the same for all glyphs but instead vary depending on the values of variables for the case (or cases) that each glyph represents.
• Frame: special aesthetics that relate the position of each glyphs in the graph to values of variables for the cases that the glyph represents.
• Scales: particular choices that determine the precise relationship between aesthetic properties and data values for glyphs.
• Guides: visual aids that help the human viewer to infer data values for cases from the aesthetic properties of the glyphs that represent them.

We will clarify these abstract ideas with a series of examples. Many of our examples will be drawn from the m111survey data frame in the tigerstats package.

help(m111survey)

You will recall that the data frame records the results of a survey of 71 students at Georgetown College in Kentucky. Each case (row in the frame) corresponds to an individual student. See Table 8.1.

Table 8.1: The first few rows of the m111survey data frame. Five variables are shown.
sex fastest GPA seat weight_feel
male 119 3.56 1_front 1_underweight
male 110 2.50 2_middle 2_about_right
female 85 3.80 2_middle 2_about_right
female 100 3.50 1_front 1_underweight
male 95 3.20 3_back 1_underweight
male 100 3.10 1_front 3_overweight

### 8.1.1 Example: a Scatterplot

We begin with a simple scatter plot based on the data. A scatter plot is often a good way to investigate graphically the relationship between two numerical variables. Figure 8.2 shows a scatter plot of student GPA vs. the fastest speed at which the student has ever driven a car.

#### 8.1.1.1 The Glyphs

In this scatter plot, the glyphs are points. Each case—each student in the survey—is represented by one and only one point on the plot.

#### 8.1.1.2 The Aesthetics

In ordinary discourse the term aesthetic refers to any perceptual property of an object. For a point, the list of its perceptual properties includes its location, its shape, its size, its color, and so on. In the Grammar of Graphics, however, only some of the properties—the one that vary from glyph to glyph depending on data—count as aesthetics in the graph.

For the scatter plot, the property of size is not considered to be an aesthetic: we can see that this is so because all of the points are the same size, and so the size cannot vary with the values of some variable in the data frame. The same goes for the property of shape: all of the points in this scatter plot are circular.

On the other hand, the property of color IS an aesthetic for the glyphs in the graph, since the males and the females in the study are represented by points of different colors. You could say that color is mapped to the variable sex in the data frame:

• the reddish color goes with the value “female”;
• the turquoise color goes with the value “male”.

#### 8.1.1.3 The Frame

In our scatter plot there are two other glyph properties that count as aesthetics:

• x-location: the position of the glyph relative to the horizontal axis of the graph;
• y-location: the position of the glyph relative to the vertical axis.

We can see that these properties are aesthetics because:

• x-location is mapped to the variable fastest: the further to the right the glyph is, the greater is the value of fastest for the student represented by that glyph.
• y-location is mapped to the variable GPA: the higher up the glyph is, the greater is the value of GPA for the student represented by that glyph.

Although x and y locations are just two more aesthetics, they are so crucial to the nature of a two-dimensional graph that they are classed separately in the Grammar of Graphics as the frame for the graph.

In the graphs we consider in this Chapter, the frame will always consist of at least the x-location, and sometimes—as in the case of our scatter plot—the y-location as well.

#### 8.1.1.4 Scales

We can decide that color (for example) is to be mapped to sex, but that decision leaves open the question of how, precisely, to make the connection. The computer can make thousands of colors: which one will correspond to the value “male”, and which to “female”? To answer that question is to choose a scale.

In this example our scale was:

• reddish = female
• turquoise = male

But we might have adopted a different scale, such as:

• blue = female
• pink = male

Every aesthetic mapping involves a choice of a scale. Consider x-location: apparently a point on the extreme left of the scatter plot represent a student who drove 50 miles per hour. A point at the extreme right corresponds to a speed of 200 miles per hour, and in general the relationship between x-location and fastest is linear: for example, a point halfway across the graph goes with a speed of 125 miles per hour, halfway between 50 mph and 200 mph. In the same way, the mapping of y-location to GPA involves a linear choice of scale.

#### 8.1.1.5 Guides

How were we able to see what the scales were for each of the three aesthetic mappings in the scatter plot? We were assisted by three set of guides, one for each mapping:

• the legend to the right of the plot guided us from color to value of sex;
• the labels and tick marks on the x-axis and the thin vertical white lines guided us from x-location to value of fastest;
• the labels and tick marks on the y-axis and the thin horizontal white lines guided us from y-location to value of GPA.

Most of the time, every aesthetic mapping is accompanied by a guide that gives the human viewer at least a rough idea of the scale chosen for that mapping.

#### 8.1.1.6 Summary

In summary we say that for this scatter plot:

• The glyphs are points.
• This time each glyph represents one and only one case.
• The frame is:
• x = fastest
• y = GPA
• Other aesthetics are:
• color = sex
• There are scales for the three aesthetic mappings above. (But we usually don’t say much about the x and y-location scales if they are linear, and we don’t make a big deal of the color scale unless we went to some trouble to choose it ourselves.)
• The legend and the axis labels, tick marks and hash-lines are the guides.

### 8.1.2 Example: Two Bar Graphs

Bar graphs are useful in the study of categorical variables, especially factor variables that have only a few possible values. Figure 8.3 shows the distribution of the factor variable seat in the mat111survey data frame.

In this graph:

• The glyphs are bars.
• This time each glyph an entire group of cases: there is a bar for the students who prefer the front, a bar for the students who prefer the middle, and a bar for the back-sitters.
• The frame is:
• x = seat. Note that it is possible for x-location to map to a categorical variable!
• In this graph the y-location does not count as part of the frame, since it is not really an aesthetic. Instead the height of a bar along the y-axis tells us how many students are represented by that bar. In the Grammar of Graphics we say that the y-axis represents a statistic—a value computed from the data. In the situation at hand, our statistic is a simple tally of the cases for each value of seat.
• There are no other aesthetics! The glyphs have various perceptual properties such as a shape rectangular and color, but these don’t vary with the cases: the shape is always rectangular and the color is always burlywood.
• There is a scale for the x-location, but there is nothing very interesting about it: the three values of seat are equally spaced along the axis.
• There is a guide for the x-location: labels on the x-axis tell us which bar goes with which value of the variable seat.

Bar graphs can also be used to study the relationship between two categorical variables. Figure 8.4 shows the relationship between sex and seating preference in the m111survey data.

Again the glyphs are bars and each bar represents many cases, but now there is a bar for each combination of the values of sex and seat. The frame is again specified only by the x-location, but this time it is mapped to sex. There is another aesthetic as well: the color (more technically, the fill) of the bars is mapped to the variable seat, allowing us to see the relationship between the two categorical variables sex and seat. Scales and guides work much the same way as in the previous example.

### 8.1.3 Examples: Histograms, Density Plots and Box Plots

In this section we’ll examine some glyphs that are useful in the visualization of numerical variables.

#### 8.1.3.1 Histograms

How are the fastest speeds driven distributed, for students in the m111survey data? In order to investigate such a question graphically, we might make a histogram like the one in Figure 8.5.

In this graph:

• The glyphs are rectangles. Each rectangle represents cases where the value of fastest lies within a particular range covered by the bottom left and right corners of the rectangle.
• The frame is:
• x = fastest.
• As with our bar graphs, the y-location does not count as part of the frame, but instead represents a statistic, permitting the height of a rectangle to indicate the number of cases that it represents.
• Again there are no other aesthetics. The burlywood fill of the rectangles is constant.
• The scale for x-location, maps location to fastest in the familiar linear fashion, and the x-axis has the usual guides found for numerical variables.

Figure 8.6 is a variant, containing a second type of glyph: each student is now represented along the X-axis by a rug-tick located approximately at his or her fastest speed. (The ticks are actually “jittered” randomly so as to avoid over-plotting when two or more students report the same speed.) The addition of a second set of glyphs is called layering, and is a common device to enhance the communicative power of a graph.

#### 8.1.3.2 Density Plots

One may also study the distribution of numerical variable with a density plot, as in Figure 8.7. In this figure there is only one glyph, the curve itself, and it represents all of the cases. However, its height represents crowding (density) of values of fastest for the cases: when the curve is high, many values are crowded closely together on the x-axis, and for speeds where the curve is low the viewer knows that few (if any) students drive around that speed. The y-axis is again used along with a statistic: for density plots the vertical scale is chosen in such a way that the total area under the density curve is 1, so that the area under the curve between two given speeds is approximately equal to the proportion of students who had speeds within that range. For density plots a rug, provided again by slightly jittered ticks, is a useful additional layer to indicate crowding of values.

Since they don’t take up much territory on a graph, density curves are especially useful when we want to study the relationship between a numerical and a categorical variable. For example, Figure 8.8 shows density plots of the fastest speeds for each of the three possible seating preferences. The glyphs are again density curves, but since the color aesthetic has been mapped to seat, we get a separate and differently-colored glyphs, one for each seating-preference..

Another approach to the same graphing problem is to use a type of glyph known as a violin. Look at Figure 8.9.

A violin glyph is simply two mirror-images of the same density plot, pasted together along their bases. Thus the violin is thick where many values are clustered together and thin where data values are rare. In this plot, the frame is constituted by mapping x-location to seat and y-location to the variable fastest. For additional communicative power we have layered another set of glyphs–jittered points, one for each case–onto the plot.

Suppose that one wished to incorporate a third variable, such as sex, into the graph? One possible way to to this is to divide the graphs into separate plots based on the values of one of the categorical variables in question. The separate plots are known as facets, and are illustrated in Figure 8.10, where facet-ing has been done on the basis of the variable seat.

#### 8.1.3.3 Box Plots

The five number summary is a convenient way to summarize the distribution of a numerical variable. The five numbers involved are:

• the minimum value
• the first quartile $$Q1$$, the $$25^{\text{th}}$$ percentile of the values
• the median, which is the $$50^{\text{th}}$$ percentile
• the third quartile $$Q3$$, the $$75^{\text{th}}$$ percentile
• the maximum value

Also of interest is the interquartile range $$IQR$$, defined as:

$IQR = Q3 - Q1.$

The interquartile range covers measures the spread in the middle 50% of the values.

In R the five number summary can be got quickly with the fivenum() function:

Galton$jitteredFather <- jitteredFather Finally, the graph: ggplot(mosaicData::Galton, aes(x = father)) + geom_density(fill = "burlywood") + geom_rug(x = jitteredFather) + labs(x = "height of the father (inches)") ## 8.3 A Case Study: US Births In Section 1.6, we made a plot of the number of births in the United States for each day of that year (see Figure 8.31). We noticed that there appear to be two clouds of points. What accounts for this phenomenon? By now we have the R-programming chops to take on this question. To begin with, look at all of the variables available in the data frame Births78: str(Births78) ## 'data.frame': 365 obs. of 8 variables: ##$ date        : Date, format: "1978-01-01" ...
##  $births : int 7701 7527 8825 8859 9043 9208 8084 7611 9172 9089 ... ##$ wday        : Ord.factor w/ 7 levels "Sun"<"Mon"<"Tue"<..: 1 2 3 4 5 6 7 1 2 3 ...
##  $year : num 1978 1978 1978 1978 1978 ... ##$ month       : num  1 1 1 1 1 1 1 1 1 1 ...
##  $day_of_year : int 1 2 3 4 5 6 7 8 9 10 ... ##$ day_of_month: int  1 2 3 4 5 6 7 8 9 10 ...
##  $day_of_week : num 1 2 3 4 5 6 7 1 2 3 ... We see that the variable wday gives the name of the day of the week, for each of the days in the year. On a hunch, we make violin plots of the births for each of the days of the week. The code appears below, and the resulting plot is shown in Figure 8.32 ggplot(Births78, aes(x = wday, y = births)) + geom_violin(fill = "burlywood") + Aha! There are considerably fewer births on the weekend-days—Saturday and Sunday. Perhaps the entire lower cloud of points is composed of weekends. Let’s check this by re-coding the days according to whether or not they are during the week or at the weekend: weekend <- with(Births78, ifelse(wday %in% c("Sat","Sun"), "weekend", "weekday")) Births78$weekend <- weekend

Note that we have added the new variable to the data frame, so that it will be easy in ggplot2 to use that variable for grouping, as in the code below. The results appear in Figure 8.33.

ggplot(Births78, aes(x = date, y = births)) + geom_point(aes(color = weekend)) +
labs(x = "Day of the Year", y = "Number of U.S. Births",
title = "Daily U.S. Birth-Numbers in 1978")

Well, a few of the points in the lower cloud are weekdays. Is there anything special about them? To find out, we subset the data frame to examine only those points:

df <- subset(Births78, weekend != "weekend" & births <= 8500)
df
##           date births wday year month day_of_year day_of_month
## 2   1978-01-02   7527  Mon 1978     1           2            2
## 149 1978-05-29   7780  Mon 1978     5         149           29
## 185 1978-07-04   8433  Tue 1978     7         185            4
## 247 1978-09-04   8481  Mon 1978     9         247            4
## 327 1978-11-23   7915  Thu 1978    11         327           23
## 359 1978-12-25   7846  Mon 1978    12         359           25
##     day_of_week weekend
## 2             2 weekday
## 149           2 weekday
## 185           3 weekday
## 247           2 weekday
## 327           5 weekday
## 359           2 weekday

If you consult a calendar for the year 1978, you will find that every one of the above days was a major holiday. Apparently doctors prefer not to deliver babies on weekend and holidays. Scheduled births—induced births or births by non-emergency Cesarean section—are not usually set for weekends or holidays. Perhaps this accounts for the two clouds we saw in the original scatter plot.

### 8.3.1 Practice Exercises

1. With mosaicData::KidsFeet, make the following graph. (Hint: Use plyr::mapvalues().)

2. With mosaicData::TenMileRace, make the following graph. (Hint: Use cut() to make age groups. The fill for the boxes is the ever-popular "burlywood".)

3. With mosaicData::Gestation, make the following graph. (Hint: Use !is.na() to select the rows of Gestation where smoke is not NA. Then use plyr::mapvalues() on smoke.)

### 8.3.2 Solutions to Practice Exercises

1. Here’s the code:

data("KidsFeet", package = "mosaicData")
tempBiggerfoot <- KidsFeet$biggerfoot biggerfoot2 <- plyr::mapvalues(tempBiggerfoot, c("L", "R"), to = c("left foot", "right foot")) KidsFeet$betterBiggerfoot <- biggerfoot2
tempDomhand <- KidsFeet$domhand domhand2 <- plyr::mapvalues(tempDomhand, from = c("L", "R"), to = c("left hand", "right hand")) KidsFeet$betterDomhand <- domhand2
ggplot(KidsFeet, aes(x = betterBiggerfoot)) +
geom_bar(aes(fill = betterDomhand), position = "dodge") +
labs(x = "the bigger foot",
title = "Foot-size and Handedness for 39 Children",
subtitle = paste0("When your right foot is bigger, ",
"are you more likely to be right-handed?"))
2. Here’s the code:

data("TenMileRace", package = "mosaicData")
tempSex <- TenMileRace$sex sex2 <- plyr::mapvalues(tempSex, from = c("F", "M"), to = c("female", "male")) TenMileRace$betterSex <- sex2
ageGroup <- cut(TenMileRace$age, breaks = c(-Inf, 20, 30, 40, 50, 60, 70, Inf), labels = c("<20", "20s", "30s", "40s", "50s", "60s", "70+")) TenMileRace$ageGroup <- ageGroup
ggplot(TenMileRace, aes(x = ageGroup, y = time)) +
geom_boxplot(fill = "burlywood") +
facet_grid(betterSex ~ .) +
labs(x = "Age Group", y = "net time to finish (sec)",
title = "Ten-Mile Race Times")
3. Here’s the code:

data("Gestation", package = "mosaicData")
Gestation <- subset(Gestation, !is.na(smoke))
tempSmoke <- Gestation$smoke smoke2 <- plyr::mapvalues(tempSmoke, from = 0:3, to = c("never", "smokes now", "until curr. preg.", "once smoked")) Gestation$betterSmoke <- smoke2
ggplot(Gestation, aes(x = factor(betterSmoke), y = wt)) +
geom_boxplot(fill = "burlywood") +
labs(x = "smoking status of mother",
y = "birth weight (ounces)",
title = "Birth Weights",
subtitle = "Children of smoking mothers have lower birth-weights.")

From time to time we will return to ggplot2 and deepen our study of this remarkable graphing system. If you are impatient to learn more right way, you can explore the package’s documentation site. The site teaches the system by way of numerous examples that you can copy and modify.

## 8.5 More in Depth

### 8.5.1 Factor Variables in Plotting

It is worthwhile to reconsider factor variables (see 7.4.3) and to consider their role in the appearance of plots. We’ll do this by way of an example.

Consider the data frame firesetting in the tigerData package:

data("firesetting", package = "tigerData")
?tigerData::firesetting

As we see from the Help file, fire-setting offers information about a large sample of high-school students. Some of them are known to be arsonists, while of course most are not. The purpose for gathering the data was to determine what characteristics of a child might be risk factors in whether he/she would develop into a fire-setter. Among the variables of interest are:

• race, the ethnicity of the student, a factor with three levels: “white”, “black”, and “other”;
• school.attitude, a scaled score on a personality inventory in which high scores indicate a poor attitude toward school;
• fires, a factor variable with levels “0” (does not set fires) and “1” (sets fires).

We might study the relationship between these variables with some box-plots via the code below. The results appear as Figure 8.34.

ggplot(firesetting, aes(x = race, y = school.attitude)) +
geom_boxplot(fill = "burlywood",
outlier.alpha = 0) +
geom_jitter(width = 0.20, size = 0.1) +
facet_wrap(~ fires, nrow = 2) +
labs(x = "Ethnicity",
y = "School-Attitude Score")

We could improve the appearance of the plot with respect to two of the variables involved:

• fires: We should display values that are more meaningful to a human viewer than 0 and 1.
• race: Let’s change the names of the values somewhat, keeping their current order along the x-axis.

For fires, we simply need to map the current values onto others that we prefer. Currently the levels are:

levels(firesetting$fires) ## NULL We can re-map as follows: betterFires <- plyr::mapvalues(firesetting$fires, from = c("0", "1"),
to = c("no fires", "sets fires"))
firesetting$betterFires <- betterFires For race, we’d like to substitute • “White” for “white”; • “AfrAm” for “black” • “Other” for “other” • “Unknown” instead of NA We’d like the order along the horizontal axis to remain as it is. We might begin by replacing the NA-values with “Unknown”, as follows: tempRace <- firesetting$race
tempRace[is.na(tempRace)] <- "Unknown"

We get an error! That’s because race is a factor with only three “possible values”, as given by its levels:

levels(tempRace)
## [1] "white" "black" "other"

“Unknown” is not considered a “possible value” for tempRace, any attempt to set that value in any elements of tempRace will be resisted.

We can overcome this resistance by coercing tempRace into a mere character-vector:

tempRace <- as.character(tempRace)

Now tempRace doesn’t come with a limit on its “possible values”. We try again:

tempRace[is.na(tempRace)] <- "Unknown"

This apparently worked. Let’s check:

unique(tempRace)
## [1] "white"   "other"   "Unknown" "black"

So far, so good. Now let’s re-map the other three values:

tempRace2 <- plyr::mapvalues(tempRace,
from = c("white", "black", "other"),
to = c("White", "AfrAm", "Other"))

Let’s check that it worked:

unique(tempRace2)
## [1] "White"   "Other"   "Unknown" "AfrAm"

Great. Now let’s add tempRace2 as a new variable into firsetting. In that data frame, it will be called betterRace. Here’s the code to do this:

firesetting$betterRace <- tempRace2 Let’s try the graph again. We will need to modify the code a bit, so that we’re using the new and improved variable betterRace, as well as the new variable betterFires. Figure 8.35 is the result. ggplot(firesetting, aes(x = betterRace, y = school.attitude)) + geom_boxplot(fill = "burlywood", outlier.alpha = 0) + geom_jitter(width = 0.20, size = 0.1) + facet_wrap(~ betterFires, nrow = 2) + labs(x = "Ethnicity", y = "School-Attitude Score") This looks better, but now the order along the x-axis isn’t as we desired. Recall that race is now a character vector only. Since it doesn’t have levels that specify a particular order, ggplot2 uses alphabetical order as the default ordering along the x-axis. Accordingly we should convert betterRace to a factor variable, and be careful to set the levels in the order that we want: tempRace <- firesetting$betterRace
betterRace2 <- factor(tempRace,
levels = c("White", "AfrAm", "Other", "Unknown"))
firesetting$evenBetterRace <- betterRace2 Now try again. Note that in the code below we switch to evenBetterRace. See Figure 8.36 for the result. ggplot(firesetting, aes(x = evenBetterRace, y = school.attitude)) + geom_boxplot(fill = "burlywood", outlier.alpha = 0) + geom_jitter(width = 0.20, size = 0.1) + facet_wrap(~ betterFires, nrow = 2) + labs(x = "Ethnicity", y = "School-Attitude Score") This works! ### 8.5.2 Example: Dental Floss If you are like me, you don’t like to run out of an item that you commonly use, so you always try to have a spare source for the item on hand. Take dental floss, for example. It’s difficult to know when the little box is almost out of floss. Hence I try to keep two boxes on hand in my bathroom drawer: the one I’m using, and a spare. The problem is that I can never tell which one is “the one I’m using” and which one is the spare. In fact on any given day I’m equally likely to pull a length of floss out of either one of the boxes.24 Let’s make the following assumptions: • A box of floss starts out with 150 feet of floss. • I start with two fresh boxes. • Every time I floss, each box has a 50% chance of being the box I pull from. • When I pull out floss, the length I choose to pull is a random number between 1 and 2 feet. If the box has less than than the amount I choose to pull, then I will just pull out whatever if left. There are several interesting questions we could ask concerning this situation. In this example, let’s focus on finding the expected value of $$X$$, the amount of floss left in the other box when I discover that one of the boxes is exhausted. I would also like to get some idea of the distribution of $$X$$—the chance for it to assume each of its possible values. Unlike many random variables that we have studied previously $$X$$ is what we call a continuous random variable: it can assume any of the values in a range of real numbers. After all, the amount left in the “other” box could be anything from 0 to 150 feet: it could be 3.5633 feet or 5.21887 feet, or any real number at all, so long as it’s between 0 and 150. Although any real number value—3.5633 feet, for instance—is possible, it seems that the chance of $$X$$ being exactly 3.5633 is vanishingly small, and the same goes for any of the possible real numbers value for $$X$$. Hence it seems that the distribution of $$X$$ cannot be described well with a table of the values that we actually get in a simulation: indeed, we are highly unlikely to get any value more than once in the course of our simulation. One way to visualize the distribution of a continuous random variable is to make a density plot of a large, representative sample of its values. Let’s practice making a density plot for a random variable that R provides for us: a so-called normal random variable. Try this code: x <- rnorm(100, mean = 70, sd = 3) We have just asked R to produce 100 numbers chosen randomly from a normal distribution. The mean parameter sets what the values should be “on average”, and the sd parameter indicates how much they should typically differ from that average value. The bigger the sd, the bigger the difference between the mean and value actually obtained is liable to be. Let’s get a look at the first ten values in the vector x: x[1:10] ## [1] 71.13092 70.90465 66.70593 66.60878 61.61040 72.16172 72.81736 ## [8] 69.31187 75.27739 70.35210 Sure enough, the numbers are typically around the mean of 70, give or take 3 units or so. We can use the ggplot2 package to produce a density curve from all 100 values in the vector x. (The resulting plot appears as Figure 8.37.) p <- ggplot(mapping = aes(x = x)) + geom_density() + geom_rug() + labs(x = "x-value") p The density plot gives us a visual sense of the distribution of the normal random variable: it seems rather likely to come out between 67 and 73 (most of the area under the curve occurs there) and rather unlikely to be greater than 76, and so on. In a similar way, we can get a sense of the distribution of the the $$X$$, the length of floss left in the remaining box, if we simulate its values many times. Since we have worked through the encapsulation-and-generalization development process a number of times by now, we will simply show a final-form implementation of a simulation function called floss_sim(). This function takes the following parameters: • reps: number of simulations (default 1000); • len: the initial length of floss in any unused box (default-value is 150); • min_pull: the minimum length I could choose to pull (default value is 1); • max_pull: the maximum length of floss that I could choose to pull (default value is 2); • graph: if TRUE, show show a density plot of the simulated values of $$X$$ (default value FALSE;) • seed: the usual option to set a random seed. The function returns an estimated of the expected value of $$X$$. floss_sim <- function(reps = 1000, len = 150, min_pull = 1, max_pull = 2, graph = FALSE, seed = NULL) { if ( !is.null(seed) ) { set.seed(seed) } # helper-function to remove floss from a box: new_length <- function(current_length) { desired_amount <- runif(1, min = min_pull, max = max_pull) max(0, current_length - desired_amount) } # vector to hold the simulated leftover-box lengths: leftover <- numeric(reps) # now the simulations. Each time we loop through, we # start with two fresh boxes and use them until # one runs out. for (i in 1:reps ) { # set up two fresh boxes: box <- c(len, len) # now floss away! repeat { # determine which box is picked picked <- sample(c(1, 2), size = 1) box[picked] <- new_length(box[picked]) if (any(box <= 0)) { remaining <- max(box) leftover[i] <- remaining break } } # end repeat-loop } # end for loop if (graph) { plot_title <- paste( "Dental-Floss Simulation with ", reps, " Repetitions", sep ="" ) sim_plot <- ggplot(mapping = aes(leftover)) + geom_density(fill = "burlywood") + labs( x = "Amount in remaining box (ft)", title = plot_title ) print(sim_plot) } # return estimate: mean(leftover) } Let’s give it a try with one thousand repetitions. floss_sim(graph = TRUE, seed = 3535) ## [1] 16.3279 A density plot of the results appears as Figure 8.38. The expected value of the amount left in the remaining box is about 16.3 feet, but you can see from the plot that there is some variability: it’s not all that unlikely to have 40 or more feet left! ### 8.5.3 Example: The Drunken Turtle Let us return to the scenario of the Drunken Turtle that was introduced in Section 5.4.5. At each step the turtle would turn by a random angle—any angle between 0 and 360 degrees with equal likelihood—and then move forward by a fixed number of units. We wondered whether the turtle would simply wander away from its starting point or perhaps eventually loop back close to its starting point. In this section we will construct a simulation that sheds light on our question. We will imagine that the turtle starts at the origin on the $$xy$$-plane and takes a large but fixed number of steps—one thousand steps, say—each of a fixed length of one unit, but each of which is preceded by a drunken turn through a random angle. We will concern ourselves with the random variable $$X$$, defined as follows: $$X$$ = the number of times in the turtle returns to within $$1/2$$ of a unit of the origin, during its first 1000 steps. (Apparently we consider a “close return” to be anything within half a unit of the origin.) In particular, we would like to know the distribution of $$X$$: what is the chance that the turtle does not make a close return during the first thousand steps? What’s the chance of making only one close return? Two close returns? And so on … . We would also like to estimate the expected value of $$X$$: the average number of close returns in the first 1000 steps. We will next consider how to code a simulation. Obviously we will need a way to record the turtle’s position at each step. In order to accomplish this we will have to apply the unit-circle trigonometry you learned in high school. Instead of measuring angles in degrees, we will work in radians: by default R does trigonometry that way. Recall that 360 degrees is $$2 \pi$$ radians; hence the way to determine the angle of a turn by the drunken turtle is with the following call to runif(): runif(1, min = 0, max = 2*pi) ## [1] 5.052785 If we want a thousand drunken turn-angles, then we can get them all in a vector, as follows: runif(1000, min = 0, max = 2*pi) So much for the turning angles. But in order to decide at each step whether the turtle has made a close return, we have to know the coordinates of the turtle after each step. This is where unit-circle trigonometry comes in handy. Recall that if a right triangle has a hypotenuse of length 1 (the length of a turtle step) and if one of the angles of the triangle is $$\alpha$$, then the side adjacent to $$\alpha$$ is $$\cos(\alpha)$$ and the side opposite is $$\sin(\alpha)$$. The situation is illustrated in Figure 8.39. From this principle it follows that if the turtle turns a random angle of $$a_1$$ radians for its first step, then after it takes the one-unit step it will be at the point $$(\cos(a_1), \sin(a_2))$$. For the second step, the turtle begins by turning another random angle, say $$a_2$$ radians. It then takes another one-unit step. This second step forms the hypotenuse of a second right triangle (see Figure 8.40). The horizontal and vertical sides of this right triangle must be $$\cos(a_2)$$ and $$\sin(a_2)$$ respectively. When the turtle complete the second step, it has added these values to its previous $$x$$ and $$y$$ coordinates, so that it is now at the point: $(\cos(a_1) + \cos(a_2), \sin(a_1) +\sin(a_2)).$ We now see the idea: in order to record the position of the turtle at the end of each step, we simply need to: • take the cosine of all the turning-angles; • take the sine of all the turning-angles; • sum the cosines to get the $$x$$-coordinate; • sum the sines to get the $$y$$-coordinate. Let’s put this idea into practice, for just the first ten steps. First, get the turning angles, the sines and the cosines: angle <- runif(10, 0 , 2*pi) xSteps <- cos(angle) ySteps <- sin(angle) Next, we need to add the xSteps together. If we do a sum, we would get the $$x$$-coordinate after the tenth step: sum(xSteps) ## [1] -1.370258 (Apparently some of the cosines were negative, meaning that some of the turning angles were between 90 and 270 degrees.) We need the $$x$$-coordinates after all of the previous steps, as well. This can be accomplished with R’s cumsum() function. You may think of cumsum as short for “cumulative sum”. Given a vector vec of numbers, cumsum(vec) returns a vector of the same length, where for every index i, cumsum(vec)[i] is sum(vec[1:i]). This is made clear in the following example: vec <- c(3, 2, -1, 4, -2) cumsum(vec) ## [1] 3 5 4 8 6 Look at the table below. vec cumulative sum cumsum(vec) 3 3 3 2 $$3+2$$ 5 -1 $$3+2+-1$$ 4 4 $$3+2+-1+4$$ 8 -2 $$3+2+-1+4+-2$$ 6 See how it works? Returning to xSteps, we see that we can get the $$x$$-coordinates of the first ten steps like this: x <- cumsum(xSteps) x ## [1] -0.6034165 -1.3905805 -2.1259607 -3.1154380 -2.4593130 ## [6] -1.5476110 -0.8591924 -1.6420445 -0.6421761 -1.3702583 Similarly we can get the first ten $$y$$-coordinates: y <- cumsum(ySteps) Using the distance formula, we can compute the distance of each point from the origin as follows: dist <- sqrt(x^2 + y^2) The following logical vector will then record, at each each step, whether or not the turtle has made a close return: closeReturn <- (dist < 0.5) The results are shown in Table 8.2. Table 8.2: Record of the turtle’s first ten steps. It has not yet made a close return. angle xSteps ySteps x y dist closeReturn 4.0646104 -0.6034165 -0.7974262 -0.6034165 -0.7974262 1.000000 FALSE 2.4769935 -0.7871641 0.6167437 -1.3905805 -0.1806826 1.402270 FALSE 3.8861615 -0.7353801 -0.6776548 -2.1259607 -0.8583374 2.292695 FALSE 2.9963954 -0.9894774 0.1446876 -3.1154380 -0.7136497 3.196131 FALSE 0.8551238 0.6561251 0.7546522 -2.4593130 0.0410024 2.459655 FALSE 0.4233886 0.9117020 0.4108522 -1.5476110 0.4518546 1.612226 FALSE 0.8114898 0.6884186 0.7253136 -0.8591924 1.1771682 1.457373 FALSE 2.4700328 -0.7828521 0.6222079 -1.6420445 1.7993761 2.435993 FALSE 0.0162276 0.9998683 0.0162269 -0.6421761 1.8156030 1.925826 FALSE 3.8968689 -0.7280822 -0.6854899 -1.3702583 1.1301131 1.776165 FALSE It remains only to package our code into a simulation function. We will allow the user to specify the fixed initial number of steps (it does not have to be 1000), and also to determine what distance counts as a “close return.” drunken_sim <- function(steps = 1000, reps = 10000, close = 0.5, seed = NULL, table = FALSE) { if ( !is.null(seed) ) { set.seed(seed) } # set up returns vector to store the number of # close returns in each repetition: returns <- numeric(reps) for (i in 1:reps) { # get all the turning angles: angle <- runif(steps, 0 , 2*pi) # all the changes in x-coordinate of turtle: xSteps <- cos(angle) # all changes in the y-coordinate: ySteps <- sin(angle) # x and y coordinates at each step of the journey: x <- cumsum(xSteps) y <- cumsum(ySteps) ## all the distances dist <- sqrt(x^2 + y^2) ## whether or not you are close to (0,0), ## at each step: is_close <- (dist < close) ## count the number of times turtle was close: close_returns <- sum(is_close) returns[i] <- close_returns } if ( table ) { cat("Here is a table of the number of close returns:\n\n") tab <- prop.table(table(returns)) print(tab) cat("\n") } mean(returns) } Let’s try it out, where the journey is goes on for one thousand steps and a return is considered to be “close” when the turtle is less than 0.5 units from its starting point. drunken_sim(steps = 1000, reps = 10000, close = 0.5, seed = 3535, table = TRUE) ## Here is a table of the number of close returns: ## ## returns ## 0 1 2 3 4 5 6 7 8 ## 0.3881 0.2290 0.1517 0.0908 0.0572 0.0335 0.0227 0.0107 0.0064 ## 9 10 11 12 13 14 15 16 17 ## 0.0048 0.0018 0.0012 0.0008 0.0007 0.0003 0.0001 0.0001 0.0001 ## [1] 1.5655 Notice that 38.81% of the time the turtle did not make a close return at all during the first thousand steps. On the other hand it is possible, though not at all likely, for the turtle to make quite a few close returns.25 It is interesting to look specifically at the distances from the origin. The following function draws a graph of the distance from the origin for the specified number of turtle steps: drunken_sim_2 <- function(steps = 1000, seed = NULL) { if ( !is.null(seed) ) { set.seed(seed) } angle <- runif(steps, 0 , 2*pi) xSteps <- cos(angle) ySteps <- sin(angle) x <- cumsum(xSteps) y <- cumsum(ySteps) dist <- sqrt(x^2 + y^2) plot_title <- paste( "First ", steps, " Steps of the Drunken Turtle", sep = "" ) dist_plot <- ggplot(mapping = aes(x = 1:steps, y=dist)) + geom_line() + labs( x = "Step", y = "Distance From Origin", title = plot_title) print(dist_plot) } Let’s try it out! The results appear in Figure 8.41 drunken_sim_2(steps = 1000, seed = 2525) It appears that although the turtle might have returned close to the origin a couple of times, for the most part it seems to be wandering further and further away. In the mathematical field known as stochastic processes, it has been shown that if the turtle is allowed to go on forever then it will certainly make a close return infinitely many times—no matter how narrowly one defines “close”! Nonetheless, as our graph suggests, the close returns tend to be spaced further and further apart as time goes on. ### 8.5.4 Practice Exercises 1. In a previous Practice Exercise (see Section 6.8.5), you wrote a function called freeThrowSim() to estimate the probability that in 20 attempts you will make more shots than Lester does. Write a new version of the function that offers the option to view a histogram of the simulated differences (number you make minus number Lester makes). ### 8.5.5 Solutions to the Practice Exercises 1. Here you go: freeThrowSim <- function(reps, seed, shots, you = 0.70, lester = 0.60, graph = FALSE) { youMake <- rbinom(reps, size = shots, prob = you) lesterMakes <- rbinom(reps, size = shots, prob = lester) youWin <- (youMake > lesterMakes) diffs <- youMake - lesterMakes df <- data.frame(diffs = diffs) if (graph) { diff_plot <- ggplot(df, aes(x = diffs)) + color = "black", fill = "skyblue", binwidth = 1 ) + labs(x = "difference (you minus Lester)") print(diff_plot) } mean(youWin) } ## The Main Ideas of This Chapter • The main features of the Grammar of Graphics: • Frame • Glyphs • Aesthetics • Scales • Guides • ggplot() makes the frame. • Add more elements to a plot with +. • Glyphs come from functions that start with geom. The following are especially common: • When you represent your data with more than one type of glyph, you are *layering”. Common layering examples: • Don’t forget to make good axis-labels with labs(). • Facet-ing is done with facet_grid() or facet_wrap(). • Often you must manipulate your data set—for example, select only certain rows from it, or add a new column to it—prior to making your graph. ## Glossary Frame The aesthetics (usually x and y position) that help locate cases on a plot. Glyph The basic graphical unit that corresponds to a case in the data table. Aesthetic A perceptible property of a glyph that varies from case to case. Scale The relationship between the value of a variable and the graphical attribute to be displayed for that value. Guide An indication, for the human viewer, of the scale being used in an aesthetic mapping. ## Exercises 1. Using the bcscr::m111survey data frame, write the ggplot2 code necessary to produce the graph in the figure below. (The points are all blue.) Hint: For the points, map the aesthetic property shape to the variable sex. On the other hand, color is a fixed property in this graph. 2. Using the mosaicData::Utilities data frame, write the ggplot2 code necessary to produce the graph in the figure below. Hint: The month variable in Utilities is given numerically: 0 for January, 1 for February, and so on. You’ll need to use plyr::mapvalues() to map the numbers from 1 to 12 to the abbreviated month names. (See section 8.5.1 for a review of this.) The special R-vector month.abb will render the re-mapping task easy: month.abb ## [1] "Jan" "Feb" "Mar" "Apr" "May" "Jun" "Jul" "Aug" "Sep" "Oct" ## [11] "Nov" "Dec" In fact, you can create a new variable with abbreviated monnths, using this code: monthAbbr <- plyr::mapvalues( Utilities$month,
from = 1:12,
to = month.abb
)

You also need to make sure that the months come in the right order along the x-axis of the graph. To do this consider resetting the levels of your new monthAbbr variable. One way to accomplish this is to convert monthAbbr back to a factor with levels in the order you want, as done in Section 8.5.1, like this:

monthAbbr <- factor(monthAbbr, levels = month.abb)

Then you can put your new variable into the data frame, naming it anything you like:

Utilities\$monthName <- monthAbbr
3. The next few exercises pertain to the data frame CPS85 from the package mosaicData. Learn about it with help(CPS85). We will use the ggplot2 graphing package to explore whether men were being paid more than women in 1985.

Make a density plot of the wages of the people in the study. As with all plots you make, it should have well-labelled axes (with units if possible). For a density plot you should label the horizontal axis, but you can let ggplot2 provide the label for the “density” axis. As always, provide a descriptive title. Also provide a “rug” of individual values along the horizontal axis.

4. Look at the plot you made in the previous exercise: you will notice that one person made a wage that was much higher than all the rest. In data analysis, when a value is much higher or lower than the rest of the values we call it an outlier.

Write the code needed to find the age, sex and sector of employment of the person who made this extraoridinarily high wage. Report the age, sex and sector of this person.

Create a new data frame called cpsSmall that is the same as CPS85 except that it excludes the row corresponding to the outlier-individual.

5. In order to explore the relationship between wage and sex in the CPS study, make violin plots for the wages of men and women. (In this exercise and in subsequent exercises, use the cpsSmall data frame so as to exclude the outlier.) Based on the plot, who tends to earn higher wages: men or women?

6. (*) Someone might argue that men don’t earn higher wages because of sex-discrimination in the workplace, but rather because of some other factor. For example, it could be that in 1985 women chose to work in low-wage sectors of the economy, whereas men tended to work in higher-wage sectors. Of course for this explanation to be viable, some sectors of the economy have to pay more on average than other sectors do. In order to verify whether this is the case, make a box plot of wage vs. sector of employment. Use the plot to name a couple of high-wage sectors and a couple of low-wage sectors.

7. (*) From the previous exercise you now know that some sectors of the economy pay more than other sectors. Hence in order to investigate properly whether there was wage-discrimination in the workforce based on sex, we would have to compare the wages of men and women who work in the same sector. To this end it would be nice to have eight separate box plots, one for each sector. Each plot would compare the wages of men and women in that sector. Use facet_wrap() to construct a graph that displays all eight plots at once.