The console

When you start R you will see a command console window, into which you type R commands. It looks like this on the Mac (and will look exactly the same inside of RStudio):

The command prompt is the > symbol.

You type R statements into the console, which is a command-line interface. R then evaluates the statements, and if it finds no errors, does something, like printing a result, producing graphics, or storing a value for later use. It then waits for you to type in another command. This is how computing was done in the olden days, and it is still the most efficient and powerful way to do computing. Something about everything old is new again, in a manner of speaking.

R commands can also be stored in a script, which is simply a text file containing basically the same commands that you would enter interactively on the command line. This is one of the most powerful features of R (R is, after all, a programming language) and enables complex sequences of commands to be stored, edited, and rerun.

R can evaluate expressions and manipulate variables, e.g.,

> x <- 2

This assigns the value of 2 to the variable x. In R-speak, this assignment operation is pronounced “x gets 2,” not “x equals 2.” Many other programming languages use an equals sign (=) for assigning values to variables (e.g., x = 2), but R is different. In R you type a less-than symbol immediately followed by a dash, making a left-pointing arrow that points towards the variable to which you want to assign a value. R will let you use an equals sign, but R purists insist that using <- is preferable. This is true both for historical reasons (explained in the previous link) and for consistency, because there is another R assignment operator that uses a double arrow, namely the <<- operator.

To print the contents of a variable, simply type the variable’s name:

> x
[1] 2

The 1 in square brackets is part of R’s way of printing lists of numbers. This convention becomes more useful when there is a longer list of numbers to print. The number in brackets is the index of the first number on that line. So if we were to generate 15 random numbers from a normal distribution, our output display might look like:

> rnorm(15)

[1] -1.20251243 -0.93915306 -0.58314349 -0.28205304 -0.72031211 1.12303378
[7] 1.60557581 1.30062736 1.06739881 -2.09506242 -0.04172139 -1.66868780
[13] 0.87027623 0.43993863 -0.07720584

Here the [7] indicates that 1.60557581 is the seventh element in the list of numbers (called a vector in R parlance).

In R, the pound symbol (#) denotes a comment. Any line that starts with # is treated as a comment and is ignored by the R interpreter.

Function calls and arguments

Much of R involves function calls, arguments, and return values. Like in your high school algebra class, when you wrote f(x) (“the f of x”), the f stands for the name of the function and the x stands for the arguments passed into the function. In algebra and in computer programming, a function is an equation or a piece of computer code that takes in one or more values (arguments), does some sort of computation with them, and then returns a result.

For example, in R we can generate a vector of 1000 random values from a normal distribution having a mean = 0 and standard deviation = 1 (the default arguments) by typing:

> y = rnorm(1000)

Here, rnorm is the name of the function and 1000 is the argument that we pass into the function. The argument tells the function that we want it to generate 1000 values. The rnorm function generates the values and returns them for storage in the variable that we have named y. If you type the name of the variable y at the R console you will get a display of those 1000 values. Here are the first and last few lines of that output:

[1] 1.3395720217 -0.0115032420 -0.6823326244 0.2613036332 0.1040746632 0.0454721524 0.3640455956 0.7205450723 0.4994317717
[10] 1.0467334635 0.0941044374 0.0470323976 1.6612519966 -1.2653165072 0.2894176277 0.2140272012 -0.1166290364 1.9260137594
[19] -0.1623506191 0.2429738277 -0.0908777982 -1.0126380527 -1.4146009241 0.6927008684 1.5923438893 0.7920350474 -0.3419451639
...
[982] 0.2657230634 0.4844520987 1.9939724168 -0.3426382437 -2.6953913082 -0.7822534469 1.0135964703 0.7795363595 0.1870213245
[991] -1.4046173709 -0.3727388529 -0.2606101406 -0.9251060911 0.1755509390 -2.5640188283 0.5750884848 0.5416196143 -0.5890929375
[1000] -1.4114665861

Instead of using the default, unspecified argument values of mean = 0 and standard deviation = 1, we can invoke rnorm like this, specifying a particular mean and sd:

> y = rnorm(1000, mean=42, sd=2.5)

If you type help(rnorm) you will get a help screen that explains the various arguments that can be used with rnorm and its relatives.

To see if we got what we wanted, we can produce a few summary statistics of the variable y:

> summary(y)

Min. 1st Qu. Median Mean 3rd Qu. Max.
34.09 40.29   41.91 41.91 43.63 50.73

We see that the mean is 41.91, which is “close enough” to 42, given that rnorm generates a different batch of random values each time it is invoked.

Some graphics

To visualize the batch of numbers stored in the variable y, we can produce a histogram:

> hist(y)

(Click image to enlarge.)

This graphic shows us the shape of the distribution of the variable y. The peak of the histogram is centered on the mean value of 42 and tapers off symmetrically on both sides of the mean. This is because the data we generated using rnorm is from a normal distribution, which by definition is in the shape of the famous bell curve.

The plot we just produced, with only one line of R code, can be easily saved to a graphics file in various formats for importing into other applications such as word processors and presentations. For example, to export the plot to a PDF file, select the graphics window and then invoke the R File|Save As… menu option to save the file.

Plots can be saved to other graphics formats as well, but you need to do it from the R console. To save as a JPEG graphics file:

> jpeg('rplot.jpg')
> hist(y)
> dev.off()

Note that the plot will be saved in the current working directory (see HINT: Command history and workspace, below) unless otherwise specified. You also may get a “null device” message after the dev.off() statement; this just means that the JPEG file (the “graphics device”) has been closed.

You can use similar commands to export R graphics to other popular graphics formats, such as bmp, tiff, and png.

As is the case with many R things, graphics management is easier in RStudio, which has a dedicated Plots panel for keeping track of graphics and exporting them in various formats, without the need for typing code into the console.

To summarize thus far, R works differently than other statistical software you may have used:

“Rather than setting up a complete analysis at once, the process is highly interactive. You run a command (say fit a model), take the results and process it through another command (say a set of diagnostic plots), take those results and process it through another command (say cross-validation), etc. The cycle may include transforming the data, and looping back through the whole process again. You stop when you feel that you have fully analyzed the data.”

Data frames

R uses an internal structure called a data frame (one of several common R data types) to store data in a row-column, spreadsheet-like format, where the rows are the objects of interest (your units of analysis) and the columns are variables measured on each object. In R, a data frame can be created by reading in raw data from an external file. A data frame can be saved by exporting it to an external file or to R’s internal data format.

To illustrate this process we will use a subset (males only) of a famous dataset from biology, the Bumpus House Sparrow data:

• Bumpus sparrow data: bumpus.dat.txt
• Bumpus sparrow metadata: bumpus.met.txt

Here is a description of the dataset, from Hermon Bumpus and House Sparrows:

… on February 1 of the present year (1898), when, after an uncommonly severe storm of snow, rain, and sleet, a number of English sparrows [= House Sparrows, Passer domesticus] were brought to the Anatomical Laboratory of Brown University [, Providence, Rhode Island]. Seventy-two of these birds revived; sixty-four perished; … ” (p. 209). “… the storm was of long duration, and the birds were picked up, not in one locality, but in several localities; …” (p. 212). This event . . . described by Hermon Bumpus (1898) . . . has served as a classic example of natural selection in action. Bumpus’ paper is of special interest since he included the measurements of these 136 birds in his paper.

For more information on Hermon Bumpus and House Sparrows, see Hermon Bumpus and Natural Selection in the House Sparrow Passer domesticus (PDF).

Creating data frames

To create an R data frame of the sparrow data:

1. Download the plain text data file bumpus.dat.txt to a readily-accessible location on your computer.
2. Issue the following command in R:

> bumpus = read.table("d:/empty/bumpus.dat.txt", sep = "", header = TRUE)

replacing the d:/empty/bumpus.txt with the path to the data file on your particular computer.

In the above read.table command, sep = "" means that there is nothing separating the data items from each other (they are only separated by white space), and header = TRUE means that the first line in the data file is a header containing the names of the variables.

Once your data are in a data frame, you can then begin to manipulate, analyze, and visualize the data. Typing names(bumpus) lists the names of the variables. Typing the name of the data frame (bumpus) will scroll all of the data to the console window. To display individual variables, use syntax like bumpus$survive. Better still, enter attach(bumpus). Once a data frame is “attached” to the console window you can just type the name of a variable by itself to display the variable’s contents. Typing fix(bumpus) will place the data frame into an editing window.

Year,Make,Model,Length
1997,Ford,E350,2.34
2000,Mercury,Cougar,2.38

age = c(25, 30, 56)
gender = c("male", "female", "male")
weight = c(160, 110, 220)
mydata = data.frame(age, gender, weight)

Scatterplot matrix

A handy graphical tool for EDA is the scatterplot matrix. This is simple to produce in R: Just give the plot function the name of a data frame:

> plot(bumpus)

(Click image to enlarge.)

A scatterplot matrix shows individual scatterplots of all pairwise combinations of variables in the data frame and is an excellent way to search for odd patterns in your data, such as groupings of data points, nonlinearities, etc.

> getOption("defaultPackages")

[1] "datasets" "utils" "grDevices" "graphics" "stats" "methods"

Summary statistics and data screening

To produce summary statistics for all variables in data frame bumpus, type:

> summary(bumpus)

    survive           length           alar           weight           lbh             lhum             lfem            ltibio          wskull           lkeel       
 Min.   :0.0000   Min.   :153.0   Min.   :39.00   Min.   :3.000   Min.   :30.00   Min.   :0.6590   Min.   :0.6580   Min.   :1.000   Min.   :0.5700   Min.   :0.7880  
 1st Qu.:0.0000   1st Qu.:158.8   1st Qu.:44.75   1st Qu.:4.875   1st Qu.:30.30   1st Qu.:0.7080   1st Qu.:0.7000   1st Qu.:1.103   1st Qu.:0.5940   1st Qu.:0.8300  
 Median :1.0000   Median :160.0   Median :47.00   Median :5.750   Median :31.00   Median :0.7360   Median :0.7090   Median :1.117   Median :0.6000   Median :0.8500  
 Mean   :0.5893   Mean   :160.0   Mean   :46.98   Mean   :5.764   Mean   :30.96   Mean   :0.7291   Mean   :0.7105   Mean   :1.123   Mean   :0.6012   Mean   :0.8506  
 3rd Qu.:1.0000   3rd Qu.:161.0   3rd Qu.:50.00   3rd Qu.:6.500   3rd Qu.:31.50   3rd Qu.:0.7460   3rd Qu.:0.7180   3rd Qu.:1.150   3rd Qu.:0.6082   3rd Qu.:0.8795  
 Max.   :1.0000   Max.   :166.0   Max.   :53.00   Max.   :8.300   Max.   :31.90   Max.   :0.7800   Max.   :0.7650   Max.   :1.197   Max.   :0.6330   Max.   :0.9160

Previously we used the summary function for a single variable. Here we give summary the name of our data frame and get summary statistics for all of the variables. survive is a binary variable, taking on only two values: 1 for survived, 0 for dead. With a binary variable the mean is not very useful: A bird can’t be 0.58 dead. However, the other variables are morphological measures that vary continuously along their scale of measurement (see bumpus.met.txt), and therefore the summary statistics can be very informative. Look at the minimum and maximum values carefully for errors, such as a misplaced decimal point. Would a House Sparrow really weigh 3000 grams (6.61 pounds)? Many researchers skip this elementary data-screening step, which is especially important if humans are typing in the data from mud-splattered field notebooks. Also, laboratory instruments can go out of calibration, batteries can become weak, power surges can fry delicate circuitry, etc., etc. Check your data!

The mean and median both measure central tendency, but in very different ways. Most people are familiar with the mean, or arithmetic average, calculated by summing the values and dividing the sum by the number of observations. The median, however, is the middle value of ranked data. In other words, if you sorted the values from low to high, the median would be the one in the middle. The mean is very sensitive to outliers (extreme values), whereas the median is not. The following R experiment will illustrate this.

> x = c(1, 2, 3, 4, 5)
> summary(x)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      1       2       3       3       4       5 

Note that the median and mean are equal. But if we add an extreme value (outlier):

> x = c(1, 2, 3, 4, 3000)
> summary(x)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      1       2       3     602       4    3000 

The median is unchanged by the huge outlier (the median is still the middle value, 3) but the mean has increased enormously. When the data are normally distributed the mean and median are equal (see Measures of Central Tendency). So another exploratory diagnostic is to compare the mean and the median. If they are very different, you may have some outliers that are skewing the distribution.

Box plots

The summary statistics 1st Qu. and 3rd Qu. are the first and third quartiles. Along with the median, the quartiles divide the dataset into four equal groups, each group comprising a quarter of the observations.

This five-number summary (minimum, maximum, median, first quartile, third quartile) is better visualised using another invention of John Tukey: the box plot:

> boxplot(length)

(Click image to enlarge.)

The line running through the central rectangle is the median. The distances from the median to the top and bottom of the rectangle span the third and first quartiles respectively. The central rectangle itself spans the distance from the first to third quartiles (the interquartile range, or IQR). Above and below the rectangle are dotted lines (called “whiskers”) that usually terminate in the minimum and maximum values. However, if outliers are present (as in this example) the whiskers terminate at 1.5 times the IQR and the outlying values are shown as dots.

A box plot is a nice way to visually assess the distribution of your data. For example, if the distribution deviates from a smooth bell-shaped curve, a box plot will make this obvious. Let’s make a skewed set of numbers using the log-normal distribution:

> q = rlnorm(1000)
> hist(q)

(Click image to enlarge.)

> boxplot(q)

(Click image to enlarge.)

The severely squished box plot instantly reveals the highly skewed nature of the data.

Stem-and-leaf display

Yet another exploratory tool popularized by Tukey is the stem-and-leaf display:

> stem(weight)

  The decimal point is at the |
  3 | 089
  4 | 0133566677899
  5 | 001456667777899
  6 | 00000013555567799
  7 | 01569
  8 | 033

The stem-and-leaf plot is best explained by looking at the variable weight as a sorted vector from low to high:

sort(weight)

[1] 3.0 3.8 3.9 4.0 4.1 4.3 4.3 4.5 4.6 4.6 4.6 4.7 4.7 4.8 4.9 4.9 5.0 5.0 5.1 5.4 5.5 5.6 5.6
[24] 5.6 5.7 5.7 5.7 5.7 5.8 5.9 5.9 6.0 6.0 6.0 6.0 6.0 6.0 6.1 6.3 6.5 6.5 6.5 6.5 6.6 6.7 6.7
[47] 6.9 6.9 7.0 7.1 7.5 7.6 7.9 8.0 8.3 8.3

The vertical line in the stem-leaf plot separates the “stems” from the “leaves,” with the stems to the left of the line. The vertical line itself is interpreted as the decimal point. So, in our sorted vector of weight, the first value is 3.0. The 3 is the stem, to the left of the vertical line, the line is the decimal point, and then the first leaf to the right of the vertical line is the zero. The next two values are 3.8 and 3.9, filling out the first row of the stem-and-leaf plot. A stem-and-leaf plot is like a histogram lying on its side, with the bars made up of the actual data values. This gives you a much more detailed view of the distribution, right down to individual values, unlike a traditional histogram in which the individual values are thrown into “bins” that are the bars of the histogram. A histogram is based on frequencies, whereas a stem-leaf plot is based on individual values. A histogram obscures the original data by summarizing it, while a stem-leaf plot displays each and every value.

Categorical variables: factors

In our bumpus data frame we have the variable survive, a categorical variable taking on only two values, a 1 indicating that the bird survived the storm and a 0 indicating that it did not survive. R calls such categorical variables factors. We can generate multiple box plots according to factor levels, such as:

> boxplot(length~survive)

(Click image to enlarge.)

The tilde symbol is used to separate the left and right sides in a formula. The formula length~survive is telling the boxplot function that we want a separate box plot of the continuous variable length for each level of the factor survive.

A variant of the box plot, and one of my favorite statistical graphics, is the notched box plot:

> boxplot(length~survive, notch=TRUE)

(Click image to enlarge.)

The notches show a 95% confidence interval around the medians. If the notches do not overlap, there is “strong evidence” that the medians differ. Our notched box plot indicates that larger birds (males in this subset of Bumpus’ data) did not survive the storm, potential evidence for natural selection. [Note, however, that the Bumpus dataset has been subjected to multiple analyses and interpretations over the years. See, for example, Hermon Bumpus and Natural Selection in the House Sparrow Passer domesticus (PDF) and Differential overnight survival by Bumpus’ House Sparrows: An alternate interpretation (PDF)].

Some confirmatory analysis

At this point we are beginning to stray into confirmatory analysis. A couple of more examples will suffice.

The classical t-test of the difference in mean body length between surviving and non-surviving male sparrows is as follows:

> t.test(length~survive)

	Welch Two Sample t-test

data:  length by survive
t = 3.2576, df = 50.852, p-value = 0.002005
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.9007984 3.7948537
sample estimates:
mean in group 0 mean in group 1 
       161.3478        159.0000 

And your basic one-way analysis of variance (ANOVA) for the same data is:

> summary(aov(length~survive))

Note that, in the above R command, we are nesting the ANOVA function (aov) inside of the summary function. This is common practice in R, and in fact in many other programming languages, where the result returned by a function is used as the input of another function. The summary function, like many functions in R, can modify its output according to what kind of input it receives.

Here are the ANOVA results:

            Df Sum Sq Mean Sq F value  Pr(>F)   
survive      1   74.7   74.71   10.16 0.00239 **
Residuals   54  397.2    7.36                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Both tests seem to indicate that there is a significant difference between the two means, with survivors being smaller than non-survivors.

However, before you believe this, read Jacob Cohen’s Things I Have Learned (So Far).