In this article by Viswa Viswanathan and Shanthi Viswanathan, the authors of the book R Data Analysis Cookbook, we discover how R can be used in various ways such as comparison, classification, applying different functions, and so on.
We will cover the following recipes:
- Creating charts that facilitate comparisons
- Building, plotting, and evaluating – classification trees
- Using time series objects
- Applying functions to subsets of a vector
(For more resources related to this topic, see here.)
Creating charts that facilitate comparisons
In large datasets, we often gain good insights by examining how different segments behave. The similarities and differences can reveal interesting patterns. This recipe shows how to create graphs that enable such comparisons.
Getting ready
If you have not already done so, download the code files and save the daily-bike-rentals.csv file in your R working directory. Read the data into R using the following command:
> bike <- read.csv("daily-bike-rentals.csv")
> bike$season <- factor(bike$season, levels = c(1,2,3,4), labels = c("Spring", "Summer", "Fall", "Winter"))
> attach(bike)
How to do it…
We base this recipe on the task of generating histograms to facilitate the comparison of bike rentals by season.
Using base plotting system
We first look at how to generate histograms of the count of daily bike rentals by season using R’s base plotting system:
- Set up a 2 X 2 grid for plotting histograms for the four seasons:
> par(mfrow = c(2,2))
- Extract data for the seasons:
> spring <- subset(bike, season == "Spring")$cnt > summer <- subset(bike, season == "Summer")$cnt > fall <- subset(bike, season == "Fall")$cnt > winter <- subset(bike, season == "Winter")$cnt
- Plot the histogram and density for each season:
> hist(spring, prob=TRUE, xlab = "Spring daily rentals", main = "") > lines(density(spring)) > > hist(summer, prob=TRUE, xlab = "Summer daily rentals", main = "") > lines(density(summer)) > > hist(fall, prob=TRUE, xlab = "Fall daily rentals", main = "") > lines(density(fall)) > > hist(winter, prob=TRUE, xlab = "Winter daily rentals", main = "") > lines(density(winter))
You get the following output that facilitates comparisons across the seasons:
Using ggplot2
We can achieve much of the preceding results in a single command:
> qplot(cnt, data = bike) + facet_wrap(~ season, nrow=2) + geom_histogram(fill = "blue")
You can also combine all four into a single histogram and show the seasonal differences through coloring:
> qplot(cnt, data = bike, fill = season)
How it works…
When you plot a single variable with qplot, you get a histogram by default. Adding facet enables you to generate one histogram per level of the chosen facet. By default, the four histograms will be arranged in a single row. Use facet_wrap to change this.
There’s more…
You can use ggplot2 to generate comparative boxplots as well.
Creating boxplots with ggplot2
Instead of the default histogram, you can get a boxplot with either of the following two approaches:
> qplot(season, cnt, data = bike, geom = c("boxplot"), fill = season)
>
> ggplot(bike, aes(x = season, y = cnt)) + geom_boxplot()
The preceding code produces the following output:
The second line of the preceding code produces the following plot:
Building, plotting, and evaluating – classification trees
You can use a couple of R packages to build classification trees. Under the hood, they all do the same thing.
Getting ready
If you do not already have the rpart, rpart.plot, and caret packages, install them now. Download the data files and place the banknote-authentication.csv file in your R working directory.
How to do it…
This recipe shows you how you can use the rpart package to build classification trees and the rpart.plot package to generate nice-looking tree diagrams:
- Load the rpart, rpart.plot, and caret packages:
> library(rpart) > library(rpart.plot) > library(caret)
- Read the data:
> bn <- read.csv("banknote-authentication.csv") - Create data partitions. We need two partitions—training and validation. Rather than copying the data into the partitions, we will just keep the indices of the cases that represent the training cases and subset as and when needed:
> set.seed(1000) > train.idx <- createDataPartition(bn$class, p = 0.7, list = FALSE)
- Build the tree:
> mod <- rpart(class ~ ., data = bn[train.idx, ], method = "class", control = rpart.control(minsplit = 20, cp = 0.01))
- View the text output (your result could differ if you did not set the random seed as in step 3):
> mod n= 961 node), split, n, loss, yval, (yprob) * denotes terminal node 1) root 961 423 0 (0.55983351 0.44016649) 2) variance>=0.321235 511 52 0 (0.89823875 0.10176125) 4) curtosis>=-4.3856 482 29 0 (0.93983402 0.06016598) 8) variance>=0.92009 413 10 0 (0.97578692 0.02421308) * 9) variance< 0.92009 69 19 0 (0.72463768 0.27536232) 18) entropy< -0.167685 52 6 0 (0.88461538 0.11538462) * 19) entropy>=-0.167685 17 4 1 (0.23529412 0.76470588) * 5) curtosis< -4.3856 29 6 1 (0.20689655 0.79310345) 10) variance>=2.3098 7 1 0 (0.85714286 0.14285714) * 11) variance< 2.3098 22 0 1 (0.00000000 1.00000000) * 3) variance< 0.321235 450 79 1 (0.17555556 0.82444444) 6) skew>=6.83375 76 18 0 (0.76315789 0.23684211) 12) variance>=-3.4449 57 0 0 (1.00000000 0.00000000) * 13) variance< -3.4449 19 1 1 (0.05263158 0.94736842) * 7) skew< 6.83375 374 21 1 (0.05614973 0.94385027) 14) curtosis>=6.21865 106 16 1 (0.15094340 0.84905660) 28) skew>=-3.16705 16 0 0 (1.00000000 0.00000000) * 29) skew< -3.16705 90 0 1 (0.00000000 1.00000000) * 15) curtosis< 6.21865 268 5 1 (0.01865672 0.98134328) *
- Generate a diagram of the tree (your tree might differ if you did not set the random seed as in step 3):
> prp(mod, type = 2, extra = 104, nn = TRUE, fallen.leaves = TRUE, faclen = 4, varlen = 8, shadow.col = "gray")
The following output is obtained as a result of the preceding command:
- Prune the tree:
> # First see the cptable > # !!Note!!: Your table can be different because of the > # random aspect in cross-validation > mod$cptable CP nsplit rel error xerror xstd 1 0.69030733 0 1.00000000 1.0000000 0.03637971 2 0.09456265 1 0.30969267 0.3262411 0.02570025 3 0.04018913 2 0.21513002 0.2387707 0.02247542 4 0.01891253 4 0.13475177 0.1607565 0.01879222 5 0.01182033 6 0.09692671 0.1347518 0.01731090 6 0.01063830 7 0.08510638 0.1323877 0.01716786 7 0.01000000 9 0.06382979 0.1276596 0.01687712 > # Choose CP value as the highest value whose > # xerror is not greater than minimum xerror + xstd > # With the above data that happens to be > # the fifth one, 0.01182033 > # Your values could be different because of random > # sampling > mod.pruned = prune(mod, mod$cptable[5, "CP"])
- View the pruned tree (your tree will look different):
> prp(mod.pruned, type = 2, extra = 104, nn = TRUE, fallen.leaves = TRUE, faclen = 4, varlen = 8, shadow.col = "gray")
- Use the pruned model to predict for a validation partition (note the minus sign before train.idx to consider the cases in the validation partition):
> pred.pruned <- predict(mod, bn[-train.idx,], type = "class")
- Generate the error/classification-confusion matrix:
> table(bn[-train.idx,]$class, pred.pruned, dnn = c("Actual", "Predicted")) Predicted Actual 0 1 0 213 11 1 11 176
How it works…
Steps 1 to 3 load the packages, read the data, and identify the cases in the training partition, respectively. In step 3, we set the random seed so that your results should match those that we display.
Step 4 builds the classification tree model:
> mod <- rpart(class ~ ., data = bn[train.idx, ], method = "class", control = rpart.control(minsplit = 20, cp = 0.01))
The rpart() function builds the tree model based on the following:
- Formula specifying the dependent and independent variables
- Dataset to use
- A specification through method=”class” that we want to build a classification tree (as opposed to a regression tree)
- Control parameters specified through the control = rpart.control() setting; here we have indicated that the tree should only consider nodes with at least 20 cases for splitting and use the complexity parameter value of 0.01—these two values represent the defaults and we have included these just for illustration
Step 5 produces a textual display of the results. Step 6 uses the prp() function of the rpart.plot package to produce a nice-looking plot of the tree:
> prp(mod, type = 2, extra = 104, nn = TRUE, fallen.leaves = TRUE, faclen = 4, varlen = 8, shadow.col = "gray")
- use type=2 to get a plot with every node labeled and with the split label below the node
- use extra=4 to display the probability of each class in the node (conditioned on the node and hence summing to 1); add 100 (hence extra=104) to display the number of cases in the node as a percentage of the total number of cases
- use nn = TRUE to display the node numbers; the root node is node number 1 and node n has child nodes numbered 2n and 2n+1
- use fallen.leaves=TRUE to display all leaf nodes at the bottom of the graph
- use faclen to abbreviate class names in the nodes to a specific maximum length
- use varlen to abbreviate variable names
- use shadow.col to specify the color of the shadow that each node casts
Step 7 prunes the tree to reduce the chance that the model too closely models the training data—that is, to reduce overfitting. Within this step, we first look at the complexity table generated through cross-validation. We then use the table to determine the cutoff complexity level as the largest xerror (cross-validation error) value that is not greater than one standard deviation above the minimum cross-validation error.
Steps 8 through 10 display the pruned tree; use the pruned tree to predict the class for the validation partition and then generate the error matrix for the validation partition.
There’s more…
We discuss in the following an important variation on predictions using classification trees.
Computing raw probabilities
We can generate probabilities in place of classifications by specifying type=”prob”:
> pred.pruned <- predict(mod, bn[-train.idx,], type = "prob")
Create the ROC Chart
Using the preceding raw probabilities and the class labels, we can generate a ROC chart:
> pred <- prediction(pred.pruned[,2], bn[-train.idx,"class"]) > perf <- performance(pred, "tpr", "fpr") > plot(perf)
Using time series objects
In this recipe, we look at various features to create and plot time-series objects. We will consider data with both a single and multiple time series.
Getting ready
If you have not already downloaded the data files, do it now and ensure that the files are in your R working directory.
How to do it…
- Read the data. The file has 100 rows and a single column named sales:
> s <- read.csv("ts-example.csv") - Convert the data to a simplistic time series object without any explicit notion of time:
> s.ts <- ts(s) > class(s.ts) [1] "ts"
- Plot the time series:
> plot(s.ts)
- Create a proper time series object with proper time points:
> s.ts.a <- ts(s, start = 2002) > s.ts.a Time Series: Start = 2002 End = 2101 Frequency = 1 sales [1,] 51 [2,] 56 [3,] 37 [4,] 101 [5,] 66 (output truncated) > plot(s.ts.a) > # results show that R treated this as an annual > # time series with 2002 as the starting year
The result of the preceding commands is seen in the following graph:
To create a monthly time series run the following command:
> # Create a monthly time series > s.ts.m <- ts(s, start = c(2002,1), frequency = 12) > s.ts.m Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec 2002 51 56 37 101 66 63 45 68 70 107 86 102 2003 90 102 79 95 95 101 128 109 139 119 124 116 2004 106 100 114 133 119 114 125 167 149 165 135 152 2005 155 167 169 192 170 180 175 207 164 204 180 203 2006 215 222 205 202 203 209 200 199 218 221 225 212 2007 250 219 242 241 267 249 253 242 251 279 298 260 2008 269 257 279 273 275 314 288 286 290 288 304 291 2009 314 290 312 319 334 307 315 321 339 348 323 342 2010 340 348 354 291 > plot(s.ts.m) # note x axis on plot
The following plot can be seen as a result of the preceding commands:
> # Specify frequency = 4 for quarterly data > s.ts.q <- ts(s, start = 2002, frequency = 4) > s.ts.q Qtr1 Qtr2 Qtr3 Qtr4 2002 51 56 37 101 2003 66 63 45 68 2004 70 107 86 102 2005 90 102 79 95 2006 95 101 128 109 (output truncated) > plot(s.ts.q)
- Query time series objects (we use the s.ts.m object we created in the previous step):
> # When does the series start? > start(s.ts.m) [1] 2002 1 > # When does it end? > end(s.ts.m) [1] 2010 4 > # What is the frequency? > frequency(s.ts.m) [1] 12
- Create a time series object with multiple time series. This data file contains US monthly consumer prices for white flour and unleaded gas for the years 1980 through 2014 (downloaded from the website of the US Bureau of Labor Statistics):
> prices <- read.csv("prices.csv") > prices.ts <- ts(prices, start=c(1980,1), frequency = 12) - Plot a time series object with multiple time series:
> plot(prices.ts)
The plot in two separate panels appears as follows:
> # Plot both series in one panel with suitable legend > plot(prices.ts, plot.type = "single", col = 1:2) > legend("topleft", colnames(prices.ts), col = 1:2, lty = 1)Two series plotted in one panel appear as follow:
How it works…
Step 1 reads the data.
Step 2 uses the ts function to generate a time series object based on the raw data.
Step 3 uses the plot function to generate a line plot of the time series. We see that the time axis does not provide much information. Time series objects can represent time in more friendly terms.
Step 4 shows how to create time series objects with a better notion of time. It shows how we can treat a data series as an annual, monthly, or quarterly time series. The start and frequency parameters help us to control these data series.
Although the time series we provide is just a list of sequential values, in reality our data can have an implicit notion of time attached to it. For example, the data can be annual numbers, monthly numbers, or quarterly ones (or something else, such as 10-second observations of something). Given just the raw numbers (as in our data file, ts-example.csv), the ts function cannot figure out the time aspect and by default assumes no secondary time interval at all.
We can use the frequency parameter to tell ts how to interpret the time aspect of the data. The frequency parameter controls how many secondary time intervals there are in one major time interval. If we do not explicitly specify it, by default frequency takes on a value of 1. Thus, the following code treats the data as an annual sequence, starting in 2002:
> s.ts.a <- ts(s, start = 2002)
The following code, on the other hand, treats the data as a monthly time series, starting in January 2002. If we specify the start parameter as a number, then R treats it as starting at the first subperiod, if any, of the specified start period. When we specify frequency as different from 1, then the start parameter can be a vector such as c(2002,1) to specify the series, the major period, and the subperiod where the series starts. c(2002,1) represent January 2002:
> s.ts.m <- ts(s, start = c(2002,1), frequency = 12)
Similarly, the following code treats the data as a quarterly sequence, starting in the first quarter of 2002:
> s.ts.q <- ts(s, start = 2002, frequency = 4)
The frequency values of 12 and 4 have a special meaning—they represent monthly and quarterly time sequences.
We can supply start and end, just one of them, or none. If we do not specify either, then R treats the start as 1 and figures out end based on the number of data points. If we supply one, then R figures out the other based on the number of data points.
While start and end do not play a role in computations, frequency plays a big role in determining seasonality, which captures periodic fluctuations.
If we have some other specialized time series, we can specify the frequency parameter appropriately. Here are two examples:
- With measurements taken every 10 minutes and seasonality pegged to the hour, we should specify frequency as 6
- With measurements taken every 10 minutes and seasonality pegged to the day, use frequency = 24*6 (6 measurements per hour times 24 hours per day)
Step 5 shows the use of the functions start, end, and frequency to query time series objects.
Steps 6 and 7 show that R can handle data files that contain multiple time series.
Applying functions to subsets of a vector
The tapply function applies a function to each partition of the dataset. Hence, when we need to evaluate a function over subsets of a vector defined by a factor, tapply comes in handy.
Getting ready
Download the files and store the auto-mpg.csv file in your R working directory. Read the data and create factors for the cylinders variable:
> auto <- read.csv("auto-mpg.csv", stringsAsFactors=FALSE)
> auto$cylinders <- factor(auto$cylinders, levels =
c(3,4,5,6,8), labels = c("3cyl", "4cyl", "5cyl", "6cyl",
"8cyl"))
How to do it…
To apply functions to subsets of a vector, follow these steps:
- Calculate mean mpg for each cylinder type:
> tapply(auto$mpg,auto$cylinders,mean) 3cyl 4cyl 5cyl 6cyl 8cyl 20.55000 29.28676 27.36667 19.98571 14.96311
- We can even specify multiple factors as a list. The following example shows only one factor since the out file has only one, but it serves as a template that you can adapt:
> tapply(auto$mpg,list(cyl=auto$cylinders),mean) cyl 3cyl 4cyl 5cyl 6cyl 8cyl 20.55000 29.28676 27.36667 19.98571 14.96311
How it works…
In step 1 the mean function is applied to the auto$mpg vector grouped according to the auto$cylinders vector. The grouping factor should be of the same length as the input vector so that each element of the first vector can be associated with a group.
The tapply function creates groups of the first argument based on each element’s group affiliation as defined by the second argument and passes each group to the user-specified function.
Step 2 shows that we can actually group by several factors specified as a list. In this case, tapply applies the function to each unique combination of the specified factors.
There’s more…
The by function is similar to tapply and applies the function to a group of rows in a dataset, but by passing in the entire data frame. The following examples clarify this.
Applying a function on groups from a data frame
In the following example, we find the correlation between mpg and weight for each cylinder type:
> by(auto, auto$cylinders, function(x) cor(x$mpg, x$weight)) auto$cylinders: 3cyl [1] 0.6191685 --------------------------------------------------- auto$cylinders: 4cyl [1] -0.5430774 --------------------------------------------------- auto$cylinders: 5cyl [1] -0.04750808 --------------------------------------------------- auto$cylinders: 6cyl [1] -0.4634435 --------------------------------------------------- auto$cylinders: 8cyl [1] -0.5569099
Summary
Being an extensible system, R’s functionality is divided across numerous packages with each one exposing large numbers of functions. Even experienced users cannot expect to remember all the details off the top of their head. In this article, we went through a few techniques using which R helps analyze data and visualize the results.
Resources for Article:
Further resources on this subject:
- Combining Vector and Raster Datasets [article]
- Factor variables in R [article]
- Big Data Analysis (R and Hadoop) [article]
