Basic R-Programming

By
Lalit Mohan Joshi
20-01-2020
 [Link]
A (very) short
and do the following (assuming you work on a
introduction to R windows computer):
• click download CRAN in the left bar
• choose a download site
Paul Torfs & Claudia Brauer • choose Windows as target operation system
Hydrology and Quantitative Water Management Group • click base
Wageningen University, The Netherlands • choose Download R 3.0.3 for Windows † and
choose default answers for all questions
3 March 2014
It is also possible to run R and RStudio from
a USB stick instead of installing them. This
could be useful when you don’t have administra-
1 Introduction tor rights on your computer. See our separate note
“How to use portable versions of R and RStudio”
R is a powerful language and environment for sta- for help on this topic.
tistical computing and graphics. It is a public do-
main (a so called “GNU”) project which is similar
2.2 Install RStudio
to the commercial S language and environment
which was developed at Bell Laboratories (for- After finishing this setup, you should see an ”R”
merly AT&T, now Lucent Technologies) by John icon on you desktop. Clicking on this would start
Chambers and colleagues. R can be considered as up the standard interface. We recommend, how-
a different implementation of S, and is much used ever, to use the RStudio interface. ‡ To install
in as an educational language and research tool. RStudio, go to:
The main advantages of R are the fact that R
is freeware and that there is a lot of help available [Link]
online. It is quite similar to other programming and do the following (assuming you work on a win-
packages such as MatLab (not freeware), but more dows computer):
user-friendly than programming languages such as • click Download RStudio
C++ or Fortran. You can use R as it is, but for • click Download RStudio Desktop
educational purposes we prefer to use R in combi- • click Recommended For Your System
nation with the RStudio interface (also freeware), • download the .exe file and run it (choose default
which has an organized layout and several extra answers for all questions)
options.
This document contains explanations, exam- 2.3 RStudio layout
ples and exercises, which can also be understood
(hopefully) by people without any programming The RStudio interface consists of several windows
experience. Going through all text and exercises (see Figure 1).
takes about 1 or 2 hours. Examples of frequently
• Bottom left: console window (also called
used commands and error messages are listed on
command window). Here you can type
the last two pages of this document and can be
simple commands after the “>” prompt and
used as a reference while programming.
R will then execute your command. This is
the most important window, because this is
2 Getting started where R actually does stuff.

2.1 Install R • Top left: editor window (also called script

window). Collections of commands (scripts)
To install R on your computer (legally for free!), can be edited and saved. When you don’t get
go to the home website of R∗ : †
At the moment of writing 3.0.3 was the latest version.
∗
On the R-website you can also find this docu- Choose the most recent one.
‡
ment: [Link] There are many other (freeware) interfaces, such as Tinn-
[Link] R.

1
1
 Figure 1 The editor, workspace, console and plots windows in RStudio.

this window, you can open it with File → you ask R to open a certain file, it will look in the
New → R script working directory for this file, and when you tell
Just typing a command in the editor window R to save a data file or figure, it will save it in the
is not enough, it has to get into the command working directory.
window before R executes the command. If Before you start working, please set your work-
you want to run a line from the script window ing directory to where all your data and script files
(or the whole script), you can click Run or are or should be stored.
press CTRL+ENTER to send it to the command Type in the command window:
window. setwd("directoryname"). For example:

• Top right: workspace / history window. > setwd("M:/Hydrology/R/")

In the workspace window you can see which
data and values R has in its memory. You Make sure that the slashes are forward slashes and
can view and edit the values by clicking on that you don’t forget the apostrophes (for the rea-
them. The history window shows what has son of the apostrophes, see section 10.1). R is case
been typed before. sensitive, so make sure you write capitals where
necessary.
• Bottom right: files / plots / packages /
Within RStudio you can also go to Tools / Set
help window. Here you can open files, view
working directory.
plots (also previous plots), install and load
packages or use the help function.
2.5 Libraries
You can change the size of the windows by drag-
ging the grey bars between the windows. R can do many statistical and data analyses. They
are organized in so-called packages or libraries.
With the standard installation, most common
2.4 Working directory
packages are installed.
Your working directory is the folder on your com- To get a list of all installed packages, go to the
puter in which you are currently working. When packages window or type library() in the console

22
window. If the box in front of the package name is You can see that a appears in the workspace win-
ticked, the package is loaded (activated) and can dow, which means that R now remembers what
be used. a is.§ You can also ask R what a is (just type a
There are many more packages available on the ENTER in the command window):
R website. If you want to install and use a pack-
> a
age (for example, the package called “geometry”)
[1] 4
you should:
• Install the package: click install packages or do calculations with a:
in the packages window and type geometry
> a * 5
or type [Link]("geometry") in the
[1] 20
command window.
• Load the package: check box in front of If you specify a again, it will forget what value
geometry or type library("geometry") in the it had before. You can also assign a new value to
command window. a using the old one.
> a = a + 10
3 Some first examples of R > a
commands [1] 14

3.1 Calculator To remove all variables from R’s memory, type

> rm(list=ls())
R can be used as a calculator. You can just type
your equation in the command window after the or click “clear all” in the workspace window. You
“>”: can see that RStudio then empties the workspace
> 10^2 + 36 window. If you only want to remove the variable
a, you can type rm(a).
and R will give the answer ToDo
[1] 136 Repeat the previous ToDo, but with several
steps in between. You can give the variables
ToDo any name you want, but the name has to start
Compute the difference between 2014 and the with a letter.
year you started at this university and divide
this by the difference between 2014 and the year
you were born. Multiply this with 100 to get
the percentage of your life you have spent at
3.3 Scalars, vectors and matrices
this university. Use brackets if you need them. Like in many other programs, R organizes num-
bers in scalars (a single number – 0-dimensional),
vectors (a row of numbers, also called arrays –
If you use brackets and forget to add the closing 1-dimensional) and matrices (like a table – 2-
bracket, the “>” on the command line changes dimensional).
into a “+”. The “+” can also mean that R is still The a you defined before was a scalar. To define
busy with some heavy computation. If you want a vector with the numbers 3, 4 and 5, you need the
R to quit what it was doing and give back the “>”, function¶ c, which is short for concatenate (paste
press ESC (see the reference list on the last page). together).
b=c(3,4,5)
3.2 Workspace
Matrices and other 2-dimensional structures
You can also give numbers a name. By doing so, will be introduced in Section 6.
they become so-called variables which can be used § Some people prefer te use <- instead of = (they do the
later. For example, you can type in the command same thing). <- consists of two characters, < and -, and
window: represents an arrow pointing at the object receiving the
value of the expression.
> a = 4 ¶
See next Section for the explanation of functions.

33
3.4 Functions
> rnorm(10, mean=1.2, sd=3.4)
If you would like to compute the mean of all the
elements in the vector b from the example above, showing that the same function (rnorm) may have
you could type different interfaces and that R has so called named
> (3+4+5)/3 arguments (in this case mean and sd). By the way,
the spaces around the “,” and “=” do not matter.
But when the vector is very long, this is very bor- Comparing this example to the previous one
ing and time-consuming work. This is why things also shows that for the function rnorm only the
you do often are automated in so-called functions. first argument (the number 10) is compulsory, and
Some functions are standard in R or in one of the that R gives default values to the other so-called
packages. You can also program your own func- optional arguments.k
tions (Section 11.3). When you use a function to RStudio has a nice feature: when you type
compute a mean, you’ll type: rnorm( in the command window and press TAB,
> mean(x=b) RStudio will show the possible arguments (Fig. 2).

Within the brackets you specify the arguments. 3.5 Plots

Arguments give extra information to the function.
In this case, the argument x says of which set R can make graphs. The following is a very sim-
∗∗ example:
of numbers (vector) the mean should computed ple
(namely of b). Sometimes, the name of the argu- 1 > x = rnorm(100)
ment is not necessary: mean(b) works as well. 2 > plot(x)

ToDo
Compute the sum of 4, 5, 8 and 11 by first com- • In the first line, 100 random numbers are
bining them into a vector and then using the assigned to the variable x, which becomes a
function sum. vector by this operation.
• In the second line, all these values are plotted
in the plots window.
The function rnorm, as another example, is a
standard R function which creates random sam- ToDo
ples from a normal distribution. Hit the ENTER Plot 100 normal random numbers.
key and you will see 10 random numbers as:
1 > rnorm(10)
2 [1] -0.949 1.342 -0.474 0.403
3 [5] -0.091 -0.379 1.015 0.740 4 Help and documentation
4 [9] -0.639 0.950
There is a large amount of (free) documentation
and help available. Some help is automatically
• Line 1 contains the command: rnorm is the func-
installed. Typing in the console window the com-
tion and the 10 is an argument specifying how
mand
many random numbers you want — in this case
10 numbers (typing n=10 instead of just 10 would > help(rnorm)
also work).
• Lines 2-4 contain the results: 10 random num- gives help on the rnorm function. It gives a de-
bers organised in a vector with length 10. scription of the function, possible arguments and
Entering the same command again produces 10 the values that are used as default for optional
new random numbers. Instead of typing the same arguments. Typing
text again, you can also press the upward arrow
key (↑) to access previous commands. If you want > example(rnorm)
10 random numbers out of normal distribution k Use the help function (Sect. 4) to see which values are
with mean 1.2 and standard deviation 3.4 you can used as default.
∗∗
type See Section 7 for slightly less trivial examples.

4
4
Figure 2 RStudio shows possible arguments when you press TAB after the function name and bracket.

gives some examples of how the function can be You can store your commands in files, the so-
used. called scripts. These scripts have typically file
An HTML-based global help can be called with: names with the extension .R, e.g. foo.R. You can
open an editor window to edit these files by click-
> [Link]()
ing File and New or Open file... †† .
You can run (send to the console window)
or by going to the help window.
part of the code by selecting lines and pressing
The following links can also be very useful:
CTRL+ENTER or click Run in the editor window. If
• [Link]
you do not select anything, R will run the line
[Link] A full manual.
your cursor is on. You can always run the whole
• [Link]
script with the console command source, so e.g.
[Link] A short reference card.
for the script in the file foo.R you type:
• [Link]
html > source("foo.R")
A very rich source of examples.
• [Link] You can also click Run all in the editor window
A typical user wiki. or type CTRL+SHIFT + S to run the whole script
• [Link] at once.
Also called Quick-R. Gives very productive
direct help. Also for users coming from other ToDo
programming languages. Make a file called firstscript.R containing R-
• [Link] code that generates 100 random numbers and
Dictionary for programming languages (e.g. R for plots them, and run this script several times.
Matlab users).
• Just using Google (type e.g. “R rnorm” in the
search field) can also be very productive.
6 Data structures
ToDo
If you are unfamiliar with R, it makes sense to just
Find help for the sqrt function. retype the commands listed in this section. Maybe
you will not need all these structures in the begin-
ning, but it is always good to have at least a first
glimpse of the terminology and possible applica-
5 Scripts tions.
R is an interpreter that uses a command line based
environment. This means that you have to type 6.1 Vectors
commands, rather than use the mouse and menus.
Vectors were already introduced, but they can do
This has the advantage that you do not always
more:
have to retype all commands and are less likely to
††
get complaints of arms, neck and shoulders. Where also the options Save and Save as are available.

5
5
 ToDo
1 > vec1 = c(1,4,6,8,10) Put the numbers 31 to 60 in a vector named
2 > vec1 P and in a matrix with 6 rows and 5 columns
3 [1] 1 4 6 8 10 named Q. Tip: use the function seq. Look at
4 > vec1[5] the different ways scalars, vectors and matrices
5 [1] 10 are denoted in the workspace window.
6 > vec1[3] = 12
7 > vec1
8 [1] 1 4 12 8 10 Matrix-operations are similar to vector opera-
9 > vec2 = seq(from=0, to=1, by=0.25) tions:
10 > vec2
11 [1] 0.00 0.25 0.50 0.75 1.00 1 > mat[1,2]
12 > sum(vec1) 2 [1] 3
13 [1] 35 3 > mat[2,]
14 > vec1 + vec2 4 [1] 2 4 6
15 [1] 1.00 4.25 12.50 8.75 11.00 5 > mean(mat)
6 [1] 4.8333

• In line 1, a vector vec1 is explicitly constructed • Elements of a matrix can be addressed in the
by the concatenation function c(), which was in- usual way: [row,column] (line 1).
troduced before. Elements in vectors can be ad- • Line 3: When you want to select a whole row,
dressed by standard [i] indexing, as shown in you leave the spot for the column number empty
lines 4-5. (the other way around for columns of course).
• In line 6, one of the elements is replaced with a • Line 5 shows that many functions also work
new number. The result is shown in line 8. with matrices as argument.
• Line 9 demonstrates another useful way of con-
structing a vector: the seq() (sequence) function.
• Lines 10-15 show some typical vector oriented 6.3 Data frames
calculations. If you add up two vectors of the
same length, the first elements of both vectors are Time series are often ordered in data frames. A
summed, and the second elements, etc., leading to data frame is a matrix with names above the
a new vector of length 5 (just like in regular vector columns. This is nice, because you can call and
calculus). Note that the function sum sums up the use one of the columns without knowing in which
elements within a vector, leading to one number position it is.
(a scalar). 1 > t = [Link](x = c(11,12,14),
2 y = c(19,20,21), z = c(10,9,7))
3 > t
6.2 Matrices 4 x y z
5 1 11 19 10
Matrices are nothing more than 2-dimensional
6 2 12 20 9
vectors. To define a matrix, use the function
7 3 14 21 7
matrix:
8 > mean(t$z)
1 mat=matrix(data=c(9,2,3,4,5,6),ncol=3) 9 [1] 8.666667
2 > mat 10 > mean(t[["z"]])
3 [,1] [,2] [,3] 11 [1] 8.666667
4 [1,] 9 3 5
• In lines 1-2 a typical data frame called t is
5 [2,] 2 4 6
constructed. The columns have the names x, y
and z.
The argument data specifies which numbers • Line 8-11 show two ways of how you can select
should be in the matrix. Use either ncol to spec- the column called z from the data frame called t.
ify the number of columns or nrow to specify the
number of rows.

6
6
 ToDo Hundred random numbers are plotted by connect-
Make a script file which constructs three ran- ing the points by lines (the symbol between quotes
dom normal vectors of length 100. Call these after the type=, is the letter l, not the number 1)
vectors x1, x2 and x3. Make a data frame called in a gold color.
t with three columns (called a, b and c) con- Another very simple example is the classical sta-
taining respectively x1, x1+x2 and x1+x2+x3. tistical histogram plot, generated by the simple
Call the following functions for this data frame: command
plot(t) and sd(t). Can you understand the
results? Rerun this script a few times. > hist(rnorm(100))

which generates the plot in Figure 3.

6.4 Lists
Histogram of rnorm(100)
Another basic structure in R is a list. The main
advantage of lists is that the “columns” (they’re

20
not really ordered in columns any more, but are

15
more a collection of vectors) don’t have to be of
the same length, unlike matrices and data frames. Frequency
10

1 > L = list(one=1, two=c(1,2),

5

2 five=seq(0, 1, length=5))
0

3 > L
−3 −2 −1 0 1 2
4 $one
rnorm(100)
5 [1] 1
6 $two
7 [1] 1 2 Figure 3 A simple histogram plot.
8 $five
9 [1] 0.00 0.25 0.50 0.75 1.00 The following few lines create a plot using the data
10 > names(L) frame t constructed in the previous ToDo:
11 [1] "one" "two" "five"
12 > L$five + 10 1 plot(t$a, type="l", ylim=range(t),
13 [1] 10.00 10.25 10.50 10.75 11.00 2 lwd=3, col=rgb(1,0,0,0.3))
3 lines(t$b, type="s", lwd=2,
• Lines 1-2 construct a list by giving names and 4 col=rgb(0.3,0.4,0.3,0.9))
values. The list also appears in the workspace 5 points(t$c, pch=20, cex=4,
window. 6 col=rgb(0,0,1,0.3))
• Lines 3-9 show a typical printing (after pressing
L ENTER).
ToDo
• Line 10 illustrates how to find out what’s in the
list. Add these lines to the script file of the previous
• Line 12 shows how to use the numbers. section. Try to find out, either by experiment-
ing or by using the help, what the meaning is of
rgb, the last argument of rgb, lwd, pch, cex.

7 Graphics
Plotting is an important statistical activity. So it To learn more about formatting plots, search
should not come as a surprise that R has many for par in the R help. Google “R color chart” for
plotting facilities. The following lines show a sim- a pdf file with a wealth of color options.
ple plot: To copy your plot to a document, go to the plots
window, click the “Export” button, choose the
> plot(rnorm(100), type="l", col="gold") nicest width and height and click Copy or Save.

7
7
8 Reading and writing data files
There are many ways to write data from within the
R environment to files, and to read data from files.
We will illustrate one way here. The following
lines illustrate the essential:

1 > d = [Link](a = c(3,4,5), Figure 4 The files [Link] of section 8 (left)

2 b = c(12,43,54)) and [Link] from the ToDo below (right)
3 > d opened in two text editors.
4 a b
5 1 3 12 9 Not available data
6 2 4 43
7 3 5 54 ToDo
8 > [Link](d, file="[Link]", Compute the mean of the square root of a vec-
9 [Link]=FALSE) tor of 100 random numbers. What happens?
10 > d2 = [Link](file="[Link]",
11 header=TRUE)
12 > d2 When you work with real data, you will en-
13 a b counter missing values because instrumentation
14 1 3 12 failed or because you didn’t want to measure in
15 2 4 43 the weekend. When a data point is not available,
16 3 5 54 you write NA instead of a number.

• In lines 1-2, a simple example data frame is > j = c(1,2,NA)

constructed and stored in the variable d.
• Lines 3-7 show the content of this data frame: Computing statistics of incomplete data sets
two columns (called a and b), each containing is strictly speaking not possible. Maybe the
three numbers. largest value occurred during the weekend when
• Line 8 writes this data frame to a text file, you didn’t measure. Therefore, R will say that it
called [Link] The argument [Link]=FALSE doesn’t know what the largest value of j is:
prevents that row names are written to the file.
Because nothing is specified about [Link], > max(j)
the default option [Link]=TRUE is chosen and [1] NA
column names are written to the file. Figure 4
shows the resulting file (opened in an editor, such If you don’t mind about the missing data and
as Notepad), with the column names (a and b) in want to compute the statistics anyway, you can
the first line. add the argument [Link]=TRUE (Should I remove
• Lines 10-11 illustrate how to read a file into the NAs? Yes!).
a data frame. Note that the column names are
also read. The data frame also appears in the > max(j, [Link]=TRUE)
workspace window. [1] 2

ToDo 10 Classes
Make a file called [Link] in Notepad from
the example in Figure 4 and store it in your The exercises you did before were nearly all with
working directory. Write a script to read it, to numbers. Sometimes you want to specify some-
multiply the column called g by 5 and to store thing which is not a number, for example the name
it as [Link]. of a measurement station or data file. In that case
you want the variable to be a character string in-
stead of a number.

8
8
 An object in R can have several so-called ToDo
classes. The most important three are numeric, Make a graph with on the x-axis: today, Sin-
character and POSIX (date-time combinations). terklaas 2014 and your next birthday and on
You can ask R what class a certain variable is by the y-axis the number of presents you expect on
typing class(...). each of these days. Tip: make two vectors first.

10.1 Characters
To tell R that something is a character string, you 11 Programming tools
should type the text between apostrophes, other-
wise R will start looking for a defined variable with When you are building a larger program than in
the same name: the examples above or if you’re using someone
else’s scripts, you may encounter some program-
> m = "apples"
ming statements. In this Section we describe a
> m
few tips and tricks.
[1] "apples"
> n = pears
Error: object ‘pears’ not found 11.1 If-statement
The if-statement is used when certain computa-
Of course, you cannot do computations with tions should only be done when a certain condi-
character strings: tion is met (and maybe something else should be
done when the condition is not met). An example:
> m + 2
Error in m + 2 : non-numeric argument to
1 > w = 3
binary operator
2 > if( w < 5 )

3 {
10.2 Dates 4 d=2
5 }else{
Dates and times are complicated. R has to know
6 d=10
that 3 o’clock comes after 2:59 and that February
7 }
has 29 days in some years. The easiest way to tell
8 > d
R that something is a date-time combination is
9 2
with the function strptime:

1 > date1=strptime( c("20100225230000", • In line 2 a condition is specified: w should be

2 "20100226000000", "20100226010000"), less than 5.
3 format="%Y%m%d%H%M%S") • If the condition is met, R will execute what is
4 > date1 between the first brackets in line 4.
5 [1] "2010-02-25 23:00:00" • If the condition is not met, R will execute what
6 [2] "2010-02-26 00:00:00" is between the second brackets, after the else in
7 [3] "2010-02-26 01:00:00" line 6. You can leave the else{...}-part out if
you don’t need it.
• In lines 1-2 you create a vector with c(...). • In this case, the condition is met and d has been
The numbers in the vectors are between apostro- assigned the value 2 (lines 8-9).
phes because the function strptime needs char- To get a subset of points in a vector for which
acter strings as input. a certain condition holds, you can use a shorter
• In line 3 the argument format specifies how the method:
character string should be read. In this case the
year is denoted first (%Y), then the month (%m), 1 > a = c(1,2,3,4)
day (%d), hour (%H), minute (%M) and second 2 > b = c(5,6,7,8)
(%S). You don’t have to specify all of them, as 3 > f = a[b==5 | b==8]
long as the format corresponds to the character 4 > f
string. 5 [1] 1 4

9
9
• In line 1 and 2 two vectors are made. 11.3 Writing your own functions
• In line 3 you say that f is composed of those
Functions you program yourself work in the same
elements of vector a for which b equals 5 or b
way as pre-programmed R functions.
equals 8.
1 > fun1 = function(arg1, arg2 )
Note the double = in the condition. Other con- 2 {
ditions (also called logical or Boolean operators) 3 w = arg1 ^ 2
are <, >, != (6=), <= (≤) and >= (≥). To test more 4 return(arg2 + w)
than one condition in one if-statement, use & if 5 }
both conditions have to be met (“and”) and | if 6 > fun1(arg1 = 3, arg2 = 5)
one of the conditions has to be met (“or”). 7 [1] 14
8
11.2 For-loop
If you want to model a time series, you usually do • In line 1 the function name (fun1) and its argu-
the computations for one time step and then for ments (arg1 and arg2) are defined.
the next and the next, etc. Because nobody wants • Lines 2-5 specify what the function should do if
to type the same commands over and over again, it is called. The return value (arg2+w) is shown
these computations are automated in for-loops. on the screen.
In a for-loop you specify what has to be done • In line 6 the function is called with arguments 3
and how many times. To tell “how many times”, and 5.
you specify a so-called counter. An example: ToDo
1 > h = seq(from=1, to=8) Write a function for the previous ToDo, so
2 > s = c() that you can feed it any vector you like
3 > for(i in 2:10) (as argument). Use a for-loop in the func-
4 { tion to do the computation with each ele-
5 s[i] = h[i] * 10 ment. Use the standard R function length
6 } in the specification of the counter. a)

7 > s
8 [1] NA 20 30 40 50 60 70 80 NA NA a
Actually, people often use more for-loops than nec-
essary. The ToDo above can be done more easily
• First the vector h is made. and quickly without a for-loop but with regular vector-
• In line 2 an empty vector ( s) is created. This is computations.
necessary because when you introduce a variable
within the for-loop, R will not remember it when
it has gotten out of the for-loop.
• In line 3 the for-loop starts. In this case, i is
the counter and runs from 2 to 10.
• Everything between the curly brackets (line 5) is
processed 9 times. The first time i=2, the second
element of h is multiplied with 10 and placed in
the second position of the vector s. The second
time i=3, etc. In the last two runs, the 9th and
10th elements of h are requested, which do not
exist. Note that these statements are evaluated
without any explicit error messages.
ToDo
Make a vector from 1 to 100. Make a for-loop
which runs through the whole vector. Multiply
the elements which are smaller than 5 and larger
than 90 with 10 and the other elements with 0.1.

10
10
12 Some useful references • max or min: largest or smallest element
• rowSums (or rowMeans, colSums and colMeans):
12.1 Functions sums (or means) of all numbers in each row (or
column) of a matrix. The result is a vector.
This is a subset of the functions explained in the • quantile(x,c(0.1,0.5)): sample the 0.1 and
R reference card. 0.5th quantiles of vector x

Data creation Data processing

• [Link]: read a table from file. Arguments: • seq: create a vector with equal steps between
header=TRUE: read first line as titles of the the numbers
columns; sep=",": numbers are separated by • rnorm: create a vector with random numbers
commas; skip=n: don’t read the first n lines. with normal distribution (other distributions are
• [Link]: write a table to file also available)
• c: paste numbers together to create a vector • sort: sort elements in increasing order
• array: create a vector, Arguments: dim: length • t: transpose a matrix
• matrix: create a matrix, Arguments: ncol • aggregate(x,by=ls(y),FUN="mean"): split
and/or nrow: number of rows/columns data set x into subsets (defined by y) and com-
• [Link]: create a data frame putes means of the subsets. Result: a new list.
• list: create a list • [Link]: interpolate (in zoo package). Ar-
• rbind and cbind: combine vectors into a gument: vector with NAs. Result: vector without
matrix by row or column NAs.
• cumsum: cumulative sum. Result is a vector.
Extracting data • rollmean: moving average (in the zoo package)
• x[n]: the nth element of a vector • paste: paste character strings together
• x[m:n]: the mth to nth element • substr: extract part of a character string
• x[c(k,m,n)]: specific elements
• x[x>m & x<n]: elements between m and n Fitting
• x$n: element of list or data frame named n • lm(v1∼v2): linear fit (regression line) between
• x[["n"]]: idem vector v1 on the y-axis and v2 on the x-axis
• [i,j]: element at ith row and jth column • nls(v1∼a+b*v2, start=ls(a=1,b=0)): non-
• [i,]: row i in a matrix linear fit. Should contain equation with variables
(here v1 and v2 and parameters (here a and b)
Information on variables with starting values
• length: length of a vector • coef: returns coefficients from a fit
• ncol or nrow: number of columns or rows in a • summary: returns all results from a fit
matrix
• class: class of a variable Plotting
• names: names of objects in a list • plot(x): plot x (y-axis) versus index number
• print: show variable or character string on the (x-axis) in a new window
screen (used in scripts or for-loops) • plot(x,y): plot y (y-axis) versus x (x-axis) in
• return: show variable on the screen (used in a new window
functions) • image(x,y,z): plot z (color scale) versus x
• [Link]: test if variable is NA (x-axis) and y (y-axis) in a new window
• [Link] or [Link]: change class to • lines or points: add lines or points to a
number or character string previous plot
• strptime: change class from character to • hist: plot histogram of the numbers in a vector
date-time (POSIX) • barplot: bar plot of vector or data frame
• contour(x,y,z): contour plot
Statistics • abline: draw line (segment). Arguments: a,b
• sum: sum of a vector (or matrix) for intercept a and slope b; or h=y for horizontal
• mean: mean of a vector line at y; or v=x for vertical line at x.
• sd: standard deviation of a vector • curve: add function to plot. Needs to have an

11
11
x in the expression. Example: curve(x^2) paste
• legend: add legend with given symbols (lty • ALT+TAB: change to another program window
or pch and col) and text (legend) at location • ↑, ↓, ← or →: move cursor
(x="topright") • HOME or END: move cursor to begin or end of line
• axis: add axis. Arguments: side – 1=bottom, • Page Up or Page Down: move cursor one page
2=left, 3=top, 4=right up or down
• mtext: add text on axis. Arguments: text • SHIFT+↑/↓/←/→/HOME/END/PgUp/PgDn: select
(character string) and side
• grid: add grid
• par: plotting parameters to be specified before 12.3 Error messages
the plots. Arguments: e.g. mfrow=c(1,3)):
number of figures per page (1 row, 3 columns); • No such file or directory or Cannot
new=TRUE: draw plot over previous plot. change working directory
Make sure the working directory and file names
Plotting parameters are correct.
These can be added as arguments to plot, lines, • Object ‘x’ not found
image, etc. For help see par. The variable x has not been defined yet. Define
• type: "l"=lines, "p"=points, etc. x or write apostrophes if x should be a character
• col: color – "blue", "red", etc string.
• lty: line type – 1=solid, 2=dashed, etc. • Argument ‘x’ is missing without default
• pch: point type – 1=circle, 2=triangle, etc. You didn’t specify the compulsory argument x.
• main: title - character string •+
• xlab and ylab: axis labels – character string R is still busy with something or you forgot
• xlim and ylim: range of axes – e.g. c(1,10) closing brackets. Wait, type } or ) or press ESC.
• log: logarithmic axis – "x", "y" or "xy" • Unexpected ’)’ in ")" or Unexpected ’}’
in "}"
The opposite of the previous. You try to close
Programming
something which hasn’t been opened yet. Add
• function(arglist){expr}: function defini-
opening brackets.
tion: do expr with list of arguments arglist
• Unexpected ‘else’ in "else"
• if(cond){expr1}else{expr2}: if-statement:
Put the else of an if-statement on the same line
if cond is true, then expr1, else expr2
as the last bracket of the “then”-part: }else{.
• for(var in vec) {expr}: for-loop: the
• Missing value where TRUE/FALSE needed
counter var runs through the vector vec and does
Something goes wrong in the condition-part
expr each run
(if(x==1)) of an if-statement. Is x NA?
• while(cond){expr}: while-loop: while cond is
• The condition has length > 1 and only
true, do expr each run
the first element will be used
In the condition-part (if(x==1)) of an if-
statement, a vector is compared with a scalar. Is
12.2 Keyboard shortcuts x a vector? Did you mean x[i]?
• Non-numeric argument to binary operator
There are several useful keyboard shortcuts for
You are trying to do computations with something
RStudio (see Help → Keyboard Shortcuts):
which is not a number. Use class(...) to find
• CRL+ENTER: send commands from script window
out what went wrong or use [Link](...) to
to command window
transform the variable to a number.
• ↑ or ↓ in command window: previous or next
• Argument is of length zero or Replacement
command
is of length zero
• CTRL+1, CTRL+2, etc.: change between the
The variable in question is NULL, which means
windows
that it is empty, for example created by c().
Check the definition of the variable.
Not R-specific, but very useful keyboard short-
cuts:
• CTRL+C, CTRL+X and CTRL+V: copy, cut and

12
12
2 Basic Arithmetic and Objects

R has a command line interface, and will accept simple commands to it. This
is marked by a > symbol, called the prompt. If you type a command and
press return, R will evaluate it and print the result for you.

> 6 + 9

[1] 15

> x <- 15
> x - 1

[1] 14

The expression x <- 15 creates a variable called x and gives it the value 15.
This is called assignment; the variable on the left is assigned to the value
on the right. The left hand side must contain only contain a single variable.

> x + 4 <- 15 # doesn't work

Assignment can also be done with = (or ->).

> x = 5
> 5*x -> x
> x

[1] 25

The operators = and <- are identical, but many people prefer <- because it
is not used in any other context, but = is, so there is less room for confusion.

2.1 Vectors

The key feature which makes R very useful for statistics is that it is vector-
ized. This means that many operations can be performed point-wise on a
vector. The function c() is used to create vectors:

13
> x <- c(1, -1, 3.5, 2)
> x

[1] 1.0 -1.0 3.5 2.0

Then if we want to add 2 to everything in this vector, or to square each

entry:

> x + 2

[1] 3.0 1.0 5.5 4.0

> x^2

[1] 1.00 1.00 12.25 4.00

This is very useful in statistics:

> sum((x - mean(x))^2)

[1] 10.69

Exercise 2.1. The weights of five people before and after a diet programme
are given in the table.
Before 78 72 78 79 105
After 67 65 79 70 93
Read the ‘before’ and ‘after’ values into two different vectors called before
and after. Use R to evaluate the amount of weight lost for each participant.
What is the average amount of weight lost?

*Exercise 2.2. How would you write a function equivalent to sum((x - mean(x))^2)
in a language like C or Java?

Some useful vectors can be created quickly with R. The colon operator is
used to generate integer sequences

> 1:10

[1] 1 2 3 4 5 6 7 8 9 10

147

14
> -3:4

[1] -3 -2 -1 0 1 2 3 4

> 9:5

[1] 9 8 7 6 5

More generally, the function seq() can generate any arithmetic progression.

> seq(from=2, to=6, by=0.4)

[1] 2.0 2.4 2.8 3.2 3.6 4.0 4.4 4.8 5.2 5.6 6.0

> seq(from=-1, to=1, length=6)

[1] -1.0 -0.6 -0.2 0.2 0.6 1.0

Sometimes it’s necessary to have repeated values, for which we use rep()

> rep(5,3)

[1] 5 5 5

> rep(2:5,each=3)

[1] 2 2 2 3 3 3 4 4 4 5 5 5

> rep(-1:3, [Link]=10)

[1] -1 0 1 2 3 -1 0 1 2 3

We can also use R’s vectorization to create more interesting sequences:

> 2^(0:10)

[1] 1 2 4 8 16 32 64 128 256 512 1024

158

15
> 1:3 + rep(seq(from=0,by=10,to=30), each=3)

[1] 1 2 3 11 12 13 21 22 23 31 32 33

The last example demonstrates recycling, which is also an important part of

vectorization. If we perform a binary operation (such as +) on two vectors
of different lengths, the shorter one is used over and over again until the
operation has been applied to every entry in the longer one. If the longer
length is not a multiple of the shorter length, a warning is given.

> 1:10 * c(-1,1)

[1] -1 2 -3 4 -5 6 -7 8 -9 10

> 1:7 * 1:2

Warning: longer object length is not a multiple of shorter object

length

[1] 1 4 3 8 5 12 7

Exercise 2.3. Create the following vectors in R using seq() and rep().

(i) 1, 1.5, 2, 2.5, . . . , 12

(ii) 1, 8, 27, 64, . . . , 1000.
(iii) 1, − 21 , 31 , − 14 , . . . , − 100
1
.
(iv) 1, 0, 3, 0, 5, 0, 7, . . . , 0, 49.
(v) 1, 3, 6, 10, 15, . . . , ni=1 i, . . . , 210 [look up ?cumsum].
P

(vi) ∗ 1, 2, 2, 3, 3, 3, 4, . . . , 9, 10, . . . , 10. [Hint: type ?seq, and read about

| {z }
10 times
the times argument.]
Exercise 2.4. The ith term in the Taylor expansion of log(1+x) is (−1)i+1 xi /i.
Create a vector containing the first 100 terms for x = 0.5. [Write out the
first few entries by hand if that helps.]
Let
n
X (−1)i+1 xi
rn (x) = log(1 + x) − .
i
i=1

Evaluate rn (1) for n = 10, 100, 1000, . . . , 106 .

16
2.2 Subsetting

It’s frequently necessary to extract some of the elements of a larger vector.

In R you can use square brackets to select an individual element or group of
elements:

> x <- c(5,9,2,14,-4)

> x[3]

[1] 2

> # note indexing starts from 1

> x[c(2,3,5)]

[1] 9 2 -4

> x[1:3]

[1] 5 9 2

> x[3:length(x)]

[1] 2 14 -4

There are two other methods for getting subvectors. The first is using a
logical vector (i.e. containing TRUE and FALSE) of the same length:

> x > 4

[1] TRUE TRUE FALSE TRUE FALSE

> x[x > 4]

[1] 5 9 14

or using negative indices to specify which elements should not be selected:

10

17
> x[-1]

[1] 9 2 14 -4

> x[-c(1,4)]

[1] 9 2 -4

(Note that this is rather different to what other languages such as C or

Python would interpret negative indices to mean.)

Exercise 2.5. The built-in vector LETTERS contains the uppercase letters
of the alphabet. Produce a vector of (i) the first 12 letters; (ii) the odd
‘numbered’ letters; (iii) the (English) consonants.

2.3 Logical Operators

As we see above, the comparison operator > returns a logical vector indi-
cating whether or not the left hand side is greater than the right hand side.
Here we demonstrate the other comparison operators:

> x <= 2 # less than or equal to

[1] FALSE FALSE TRUE FALSE TRUE

> x == 2 # equal to

[1] FALSE FALSE TRUE FALSE FALSE

> x != 2 # not equal to

[1] TRUE TRUE FALSE TRUE TRUE

Note the double equals sign ==, to distinguish between assignment and com-
parison.
We may also wish to combine logical vectors. If we want the elements of x
within a range, we can use the following:

11

18
> (x > 0) & (x < 10) # 'and'

[1] TRUE TRUE TRUE FALSE FALSE

The & operator does a pointwise ‘and’ comparison between the two sides.
Similarly, the vertical bar | does pointwise ‘or’, and the unary ! operator
performs negation.

> (x == 5) | (x > 10)

[1] TRUE FALSE FALSE TRUE FALSE

> !(x > 5)

[1] TRUE FALSE TRUE FALSE TRUE

Exercise 2.6. The function rnorm() generates normal random variables.

For instance, rnorm(10) gives a vector of 10 i.i.d. standard normals. Gen-
erate 20 standard normals, and store them as x. Then obtain subvectors
of

(i) the entries in x which are less than 1;

(ii) the entries between − 12 and 1;

(iii) the entries whose absolute value is larger than 1.5.

2.4 Character Vectors

As you might have noticed in the exercise above, vectors don’t have to
contain numbers. We can equally create a character vector, in which
each entry is a string of text. Strings in R are contained within double
quotes ":

> x <- c("Hello", "how do you do", "lovely to meet you", 42)
> x

[1] "Hello" "how do you do" "lovely to meet you"

[4] "42"

12

19
Notice that you cannot mix numbers with strings: if you try to do so the
number will be converted into a string. Otherwise character vectors are
much like their numerical counterparts.

> x[2:3]

[1] "how do you do" "lovely to meet you"

> x[-4]

[1] "Hello" "how do you do" "lovely to meet you"

> c(x[1:2], "goodbye")

[1] "Hello" "how do you do" "goodbye"

2.5 Matrices

Matrices are much used in statistics, and so play an important role in R. To

create a matrix use the function matrix(), specifying elements by column
first:

> matrix(1:12, nrow=3, ncol=4)

[,1] [,2] [,3] [,4]

[1,] 1 4 7 10
[2,] 2 5 8 11
[3,] 3 6 9 12

This is called column-major order. Of course, we need only give one of

the dimensions:

> matrix(1:12, nrow=3)

unless we want vector recycling to help us:

> matrix(1:3, nrow=3, ncol=4)

[,1] [,2] [,3] [,4]

13

20
[1,] 1 1 1 1
[2,] 2 2 2 2
[3,] 3 3 3 3

Sometimes it’s useful to specify the elements by row first

> matrix(1:12, nrow=3, byrow=TRUE)

There are special functions for constructing certain matrices:

> diag(3)

[,1] [,2] [,3]

[1,] 1 0 0
[2,] 0 1 0
[3,] 0 0 1

> diag(1:3)

[,1] [,2] [,3]

[1,] 1 0 0
[2,] 0 2 0
[3,] 0 0 3

> 1:5 %o% 1:5

[,1] [,2] [,3] [,4] [,5]

[1,] 1 2 3 4 5
[2,] 2 4 6 8 10
[3,] 3 6 9 12 15
[4,] 4 8 12 16 20
[5,] 5 10 15 20 25

The last operator performs an outer product, so it creates a matrix with

(i, j)-th entry xi yj . The function outer() generalizes this to any function
f on two arguments, to create a matrix with entries f (xi , yj ). (More on
functions later.)

> outer(1:3, 1:4, "+")

[,1] [,2] [,3] [,4]

14

21
[1,] 2 3 4 5
[2,] 3 4 5 6
[3,] 4 5 6 7

Matrix multiplication is performed using the operator %*%, which is quite

distinct from scalar multiplication *.

> A <- matrix(c(1:8,10), 3, 3)

> x <- c(1,2,3)
> A %*% x # matrix multiplication

[,1]
[1,] 30
[2,] 36
[3,] 45

> A*x # NOT matrix multiplication

[,1] [,2] [,3]

[1,] 1 4 7
[2,] 4 10 16
[3,] 9 18 30

Standard functions exist for common mathematical operations on matrices.

> t(A) # transpose

[,1] [,2] [,3]

[1,] 1 2 3
[2,] 4 5 6
[3,] 7 8 10

> det(A) # determinant

[1] -3

> diag(A) # diagonal

[1] 1 5 10

15

22
> solve(A) # inverse

[,1] [,2] [,3]

[1,] -0.6667 -0.6667 1
[2,] -1.3333 3.6667 -2
[3,] 1.0000 -2.0000 1

Exercise 2.7. Construct the matrix

 
1 2 3
B= 4 2 6
−3 −1 −3

Show that B × B × B is a scalar multiple of the identity matrix, and find

the scalar.

Matrices can be subsetted much the same way as vectors, although of course
they have two indices. Row number comes first:

> A[2,1]

[1] 2

> A[2,2:ncol(A)]

[1] 5 8

> A[,1:2] # blank indices give everything

[,1] [,2]
[1,] 1 4
[2,] 2 5
[3,] 3 6

> A[c(),1:2] # empty indices give nothing!

[,1] [,2]

Notice that, where appropriate, R automatically reduces a matrix to a vector

or scalar when you subset it. You can override this using the optional drop
argument.

16

23
> A[2,2:ncol(A),drop=FALSE] # returns a matrix

[,1] [,2]
[1,] 5 8

You can stitch matrices together using the rbind() and cbind() functions.
These employ vector recycling:

> cbind(A, t(A))

[,1] [,2] [,3] [,4] [,5] [,6]

[1,] 1 4 7 1 2 3
[2,] 2 5 8 4 5 6
[3,] 3 6 10 7 8 10

> rbind(A, 1, 0)

[,1] [,2] [,3]

[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 10
[4,] 1 1 1
[5,] 0 0 0

Exercise 2.8. Construct the following matrices:

(a)  
1 3 5 7
2 4 6 8

(b)  
1 −1 1 · · · −1
1 −1 1 · · · −1
(dimensions 15 × 10).
 
 .. .. .. 
. . . 
1 −1 1 · · · −1

(c) The 5 × 15 matrix with three 1s in shifting positions:

 
1 1 1 0 0 ··· 0 0
0 0 0 1 1 · · · 0 0
(dimensions 5 × 15).
 
 .. .. .. .. 
. . . .
0 0 0 0 0 ··· 1 1

17

24
 [Hint: use column subsetting.]

(d)  
1 2 3 · · · 9 10
2 3 4
 · · · 10 11 
 .. 
3 4 5 .
 ;
 .. .. .. 
. . . 17
 
 9 10 17 18
10 11 · · · 17 18 19
[Look at the outer() function.]

(e)  
1 2 3 ··· 9
4
2 3 4 1
 
 .. .. 
3 4 . .
;
 
 ..
4 . 6
 ..
 

. 6 7
9 1 ··· 6 7 8
[The modular arithmetic operator %% may be useful here.]

(f)  
I5 1
0 −I6
where Ik is the k × k-identity matrix, and 1 and 0 are matrices with all
entries 1 and 0 respectively.

Exercise 2.9. Solve the following system of simultaneous equations using

matrix methods.

a + 2b + 3c + 4d + 5e = −5
2a + 3b + 4c + 5d + e = 2
3a + 4b + 5c + d + 2e = 5
4a + 5b + c + 2d + 3e = 10
5a + b + 2c + 3d + 4e = 11

Don’t just create your matrix by hand!

Exercise 2.10. In this section we’ve seen that the behaviour of the function
diag() depends upon its inputs. Can you think of some examples where
this might cause a problem?

18

25
2.6 Lists

Other than vectors and matrices, the main object for holding data in R is a
list1 . These are a bit like vectors, except that each entry can be any other
R object, even another list.

> x <- list(1:3, TRUE, "Hello", list(1:2, 5))

Here x has 4 elements: a numeric vector, a logical, a string and another list.
We can select an entry of x with double square brackets:

> x[[3]]

[1] "Hello"

To get a sub-list, use single brackets:

> x[c(1,3)]

[[1]]
[1] 1 2 3

[[2]]
[1] "Hello"

Notice the difference between x[[3]] and x[3].

We can also name some or all of the entries in our list, by supplying argu-
ment names to list():

> x <- list(y=1:3, TRUE, z="Hello")

> x

$y
[1] 1 2 3

[[2]]
[1] TRUE

$z
[1] "Hello"
1
Technically speaking, lists are also a kind of vector in R, but not every object in them
has to have the same type; ordinary logical, numeric or character vectors are known as
atomic vectors.

19

26
Notice that the [[1]] has been replaced by $y, which gives us a clue as to
how we can recover the entries by their name. We can still use the numeric
position if we prefer:

> x$y

[1] 1 2 3

> x[[1]]

[1] 1 2 3

The function names() can be used to obtain a character vector of all the
names of objects in a list.

> names(x)

[1] "y" "" "z"

You’ve seen most standard R objects now: almost all the more complicated
ones are just lists! We’ll see this in the next section.

20

27
3 Data

R comes with many datasets built-in, particularly in the MASS package. A

package is a collection (or library) of functions, datasets, and other objects;
most packages are not loaded automatically, so you have to do it yourself:

> library(MASS)

You can now access various datasets from this package. Try looking at the
dataset called hills.

> head(hills)

To find out what the data represent, use the help function ?hills.

3.1 Data Frames

The object hills is something called a data frame. A data frame is a

series of records represented by rows (in this case one for each race), each
containing values in several fields (in this case dist, climb, time).
You can check that hills is a data frame by inspecting its class(es)

> class(hills)

[1] "[Link]"

or more reliably by using an is() command

> is(hills, "[Link]")

[1] TRUE

We’ll talk more about classes later in the course.

Data frames share many of the characteristics of matrices. We can select
rows or columns in the same way:

> hills[3,]

dist climb time

Craig Dunain 6 900 33.65

21

28
> hills[hills$dist >= 12,]

dist climb time

Bens of Jura 16 7500 204.62
Lairig Ghru 28 2100 192.67
Seven Hills 14 2200 98.42
Two Breweries 18 5200 170.25
Moffat Chase 20 5000 159.83

However, they also behave like lists indexed by the columns:

> hills$time

[1] 16.08 48.35 33.65 45.60 62.27 73.22 204.62 36.37 29.75 39.75
[11] 192.67 43.05 65.00 44.13 26.93 72.25 98.42 78.65 17.42 32.57
[21] 15.95 27.90 47.63 17.93 18.68 26.22 34.43 28.57 50.50 20.95
[31] 85.58 32.38 170.25 28.10 159.83

The truth is that, like almost all complicated objects in R, data frames
are lists with some additional structure. Formally speaking, they are not
matrices, but they do behave similarly in certain circumstances.

Exercise 3.1. How do the results of the following commands differ from
what we would expect if hills were a matrix?

> hills[1,]
> hills[3]
> hills %*% c(1,2,4)
> mean(hills)

3.2 Manipulating Data using with()

We often want to use functions on the columns of a data frame, and it quickly
becomes inconvenient to repeatedly type (for example) hills$ before every
such event. For example, the command below will give a scatter plot of the
race times against climbs, amongst only those races less than 10 miles long.

> plot(hills$climb[hills$dist < 10], hills$time[hills$dist < 10])

22

29
The with() function allows us to refer to the names of objects in a data
frame (or, in fact, any list) without having to keep referring to the data
frame itself. For example, the command above becomes

> with(hills, plot(climb[dist < 10], time[dist < 10]))

If you just type climb or dist on their own, R won’t know what object
you’re referring to. Technically with() alters the scope of the expression
being evaluated (i.e. the code given in the second argument) so that it can
‘see’ the columns of the data frame as objects. We’ll learn a bit more about
scope when we talk about functions later on.

Exercise 3.2. Using with(), find the mean of the average speeds (in miles
per hour) for races which are between 5 and 10 miles long

3.3 Creating Data Frames

The command [Link]() is used to create a data frame, each argument

representing a column.

> books <- [Link](author=c("Ripley", "Cox", "Snijders", "Cox"),

+ year=c(1980, 1979, 1999, 2006),
+ publisher=c("Wiley", "Chapman", "Sage", "CUP"))
> books

author year publisher

1 Ripley 1980 Wiley
2 Cox 1979 Chapman
3 Snijders 1999 Sage
4 Cox 2006 CUP

Exercise 3.3. (a) Create a small data frame representing a database of

films. It should contain the fields title, director, year, country, and
at least three films.

(b) Create a second data frame of the same format as above, but containing
just one new film.

(c) Merge the two data frames using rbind().

(d) Try sorting the titles using sort(): what happens?

23

30
3.4 Factors

There are two main types of data which you will encounter this year: nu-
merical and categorical. We’ve seen how to create numerical vectors already.
Suppose we have the heights of 100 individuals, the first 50 male and the
rest female.

> [Link](1442) # fixes the random numbers

> height = round(rnorm(100, mean=rep(c(170,160),each=50), sd=10))
> sex = rep(c("M", "F"), each=50)
> head(sex)

[1] "M" "M" "M" "M" "M" "M"

We can tell R to treat sex as a categorical variable:

> Sex = [Link](sex)

> head(Sex)

[1] M M M M M M
Levels: F M

Note that it is displayed slightly differently. The new variable Sex is called a
factor; a factor is a categorical variable which takes various discrete levels,
in this case M and F for male and female.
R knows to do sensible things with factors:

> plot(Sex, height)

24

31
 200
190
180
170
160
150
140

F M

What happens if you try to plot sex against height instead? The distinction
between categorical and non-categorical data is especially important if we
have numbered groups.
The information in a factor is stored as a vector of integers:

> [Link](Sex)

[1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[36] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[71] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

Just as a data frame is really a list, a factor is really a vector of integers

(for levels) together with some extra information giving each level a names.
The additional information is contained within a list of attributes. You
can view this list directly.

> attributes(Sex)

$levels
[1] "F" "M"

25

32
$class
[1] "factor"

The attributes in this case are its class (you’ll see this in many objects)
and a vector of the level names. The class tells R that this object should be
treated as a factor so that, for example, it will be displayed to you in the
right way.
You may find that sometimes data are stored as a factor when you don’t
want them to be (see the exercise in the previous section). You can turn a
factor back in to a character vector easily enough:

> [Link](Sex)

[1] "M" "M" "M" "M" "M" "M" "M" "M" "M" "M" "M" "M" "M" "M" "M" "M" "M"
[18] "M" "M" "M" "M" "M" "M" "M" "M" "M" "M" "M" "M" "M" "M" "M" "M" "M"
[35] "M" "M" "M" "M" "M" "M" "M" "M" "M" "M" "M" "M" "M" "M" "M" "M" "F"
[52] "F" "F" "F" "F" "F" "F" "F" "F" "F" "F" "F" "F" "F" "F" "F" "F" "F"
[69] "F" "F" "F" "F" "F" "F" "F" "F" "F" "F" "F" "F" "F" "F" "F" "F" "F"
[86] "F" "F" "F" "F" "F" "F" "F" "F" "F" "F" "F" "F" "F" "F" "F"

Exercise 3.4. (a) Sample the numbers {1, 2, 3} uniformly with replace-
ment 50 times; use this to create a factor with levels Yes, No and Maybe.

(b) Create a subvector by removing the Maybe entries from the factor above.
What levels does the new factor have?

(c) Use the command droplevels() to remove the level Maybe.

Exercise 3.5. Take a look at the birthwt data from the MASS package.
How is race stored in these data? Is this sensible?
Define a factor based on race:

> Race = factor(birthwt$race)

Compare the effect of the commands summary(), plot() and mean() on

each of Race and birthwt$race. Which do you find more useful?

3.5 Row and Column Names

The labels above and to the left of the values in hills are not part of the
data itself, but can be accessed:

26

33
> names(hills)

[1] "dist" "climb" "time"

> [Link](hills)

[1] "Greenmantle" "Carnethy" "Craig Dunain"

[4] "Ben Rha" "Ben Lomond" "Goatfell"
[7] "Bens of Jura" "Cairnpapple" "Scolty"
[10] "Traprain" "Lairig Ghru" "Dollar"
[13] "Lomonds" "Cairn Table" "Eildon Two"
[16] "Cairngorm" "Seven Hills" "Knock Hill"
[19] "Black Hill" "Creag Beag" "Kildcon Hill"
[22] "Meall Ant-Suidhe" "Half Ben Nevis" "Cow Hill"
[25] "N Berwick Law" "Creag Dubh" "Burnswark"
[28] "Largo Law" "Criffel" "Acmony"
[31] "Ben Nevis" "Knockfarrel" "Two Breweries"
[34] "Cockleroi" "Moffat Chase"

As we saw above, in a data frame the column names can be used for indexing
(e.g. hills$time); the row names cannot be used in this way.
This additional information is stored as attributes, which are in a separate
list2 attached to the object hills:

> attributes(hills)

We could add an attribute to hills if we wanted:

> attributes(hills) <- c(attributes(hills), list(type="races"))

> attributes(hills)

Note that type(hills) doesn’t access hills$type, the most important

attributes such as names and class happen to have functions named after
them, which can be used to extract relevant information.

3.6 Reading in Data

Whilst R provides many interesting datasets, it is often necessary to load

data externally. The main commands for doing this are [Link]() and
[Link]().
2
Actually they’re not stored as a list (see ?attributes), but they behave very similarly.

27

34
Here is an example using the [Link] dataset, available on the course
website.

> dat <- [Link]("[Link]", header=TRUE)

> head(dat)

STATE CIG BLAD LUNG KID LEUK

1 AL 18.20 2.90 17.05 1.59 6.15
2 AZ 25.82 3.52 19.80 2.75 6.61
3 AR 18.24 2.99 15.98 2.02 6.94
4 CA 28.60 4.46 22.07 2.66 7.06
5 CT 31.10 5.11 22.83 3.35 7.20
6 DE 33.60 4.78 24.55 3.36 6.45

> class(dat)

[1] "[Link]"

What happens if header=TRUE is omitted?

When you specify the file name, be sure to use the double quotes (") around
it. You also need to give the correct path to the file. R will automatically
look for the file in its working directory. You can check what this is:

> getwd()

[1] "/data/redcrest/evans/Dropbox/Teaching/R Programming/2014"

Then if your file is in a subfolder called files, you need to write (for exam-
ple)

> dat <- [Link]("files/[Link]", header=TRUE)

In some systems you can use [Link]() to get the full path to a file.
In particular this works well on R GUI for Windows or OS X. For example:

> dat <- [Link]([Link](), header=TRUE)

will bring up a window for you to choose the file from.

*Exercise 3.6. Look at the documentation for [Link](). Use the
function to read in only lines 11 to 20 (Indiana up to Minnesota).

28

35
4 Functions

Everything which is done in R is done by functions. A function in a pro-

gramming langauge is much like its mathematical equivalent: it has some
inputs called arguments, and an output called the return value. In R a
function can only return a single object. If you type a function’s name at
the console, you can see its structure:

> setdiff

function (x, y)
{
x <- [Link](x)
y <- [Link](y)
unique(if (length(x) || length(y))
x[match(x, y, 0L) == 0L]
else x)
}
<bytecode: 0x2a49e98>
<environment: namespace:base>

There are two important parts to the function: the signature, which in
this case is function(x, y), and the body, which is the code between the
curly brackets. Broadly speaking, when a function is called, it takes the
information in the arguments, applies the code in the body to them, and
then spits out the final expression in the function. In this case that’s the
complex looking expression unique( ).

4.1 Arguments

You can look at the arguments for a function by typing args:

> args(setdiff)

function (x, y)
NULL

Arguments are a little complicated in R. You’ll notice that they have a

name: the arguments of setdiff() are called x and y. However, you don’t
usually have to specify an argument by name, because arguments also have
a position:

29

36
> a <- c(1,4,5,7)
> b <- c(1,2,5,9)
> setdiff(a, b) # everything in 'a' and not in 'b'

[1] 4 7

It assumes that the first argument supplied should be x, and the second y.
You can override this by specifying the name, and then the order doesn’t
matter:

> setdiff(y=b, x=a)

[1] 4 7

> setdiff(b, a)

[1] 2 9

If you specify some of the argument names but not all, then it will use the
ordering to deduce the others.

> setdiff(y=b, a) # 'a' must be the argument 'x'

[1] 4 7

Most functions don’t require all of their arguments to be specified.

> x <- rnorm(10)

> y <- x + rnorm(10)
> lm(y ~ x)

Call:
lm(formula = y ~ x)

Coefficients:
(Intercept) x
0.0212 0.7415

> args(lm)

30

37
function (formula, data, subset, weights, [Link], method = "qr",
model = TRUE, x = FALSE, y = FALSE, qr = TRUE, [Link] = TRUE,
contrasts = NULL, offset, ...)
NULL

The function lm() (which fits a linear model) only requires the single argu-
ment formula to run; the other arguments are optional: some of them have
default values which are shown in the signature (such as model = TRUE),
whereas others simply alter the behaviour of the function when specified
(such as data).
R also uses partial matching for arguments, so as long as you give enough
of the argument’s name to make in unambiguous which one you mean, it
will work:

> ## these two will do the same thing

> lm(form = y ~ x)
> lm(formula = y ~ x)

4.2 Writing Functions

To define your own function you just have to construct something in the
same format as above:

> square = function(x) {

+ x^2
+ }
> square(4)

[1] 16

Objects which are created inside a function do not exist outside it:

> mean2 <- function(x) {

+ n <- length(x)
+ sum(x)/n
+ }
>
> mean2(1:10)

[1] 5.5

31

38
> n

Error: object ’n’ not found

Clearly an object called n was used inside the function above, but it was
only inside the function’s namespace. Most functions in R do not have side
effects: they return a value, but do not change any of the objects which
you can reach at the console. In order to use a function, you usually have
to assign its output to something.

> x <- mean2(1:10)

> x

[1] 5.5

Exercise 4.1. The logit function is defined as

 
x
logit(x) = log , 0 < x < 1.
1−x

Write an R function in one argument to implement this. How does your

function behave for values of x such as 0, 1, or 2?

Exercise 4.2. Recall that the Taylor expansion of log(1 + x) is

∞
X xi
log(1 + x) = (−1)i+1
i
i=1

Write a function with arguments x and n, which calculates the Taylor ap-
proximation to log(1 + x) using n terms.
How many terms do you need to get within 10−6 of the correct solution
when x = 0.99?

Exercise 4.3. Given real vectors x, y of length n, the least squares slope
(α, β)T is given by
P
(xi − x̄)(yi − ȳ)
β = iP 2
i (xi − x̄)
α = ȳ − x̄β.

Write a function which takes two arguments, x and y, and returns a vector
of length 2 containing α and β. Verify that your function gives the correct
answer using R’s built-in function lm() [the syntax is lm(y~x)].

32

39
4.3 for() Loops

The most common way to execute a block of code multiple times is with a
for() loop. What’s going on in the code below?

> factorial2 = function(n) {

+ out = 1

+ for (i in 1:n) {
+ out = out*i
+ }

+ out
+ }
> factorial2(10)

[1] 3628800

You may have seen for() loops in other languages. The syntax in R is for
(i in x) for some vector (or list) x, where i will take each value in x. Most
commonly, x is a vector of the first n natural numbers.
i is a dummy variable, and can be called whatever you like, though it retains
its value outside the loop.

> for (sillyname in 1:4) print(sillyname)

[1] 1
[1] 2
[1] 3
[1] 4

> sillyname

[1] 4

Exercise 4.4. Write a function to perform matrix-vector multiplication. It

should take a matrix A and a vector b as arguments, and return the vector
Ab. Use two loops to do this, rather than %*% or any vectorization.

I generally recommend using seq_len() or seq_along() in for() loops,

because it always behaves the way you want (and runs quicker than seq()):

33

40
> n = 0
> 1:n # not a sequence of length n=0

[1] 1 0

> seq(n) # ditto

[1] 1 0

> seq_len(n) # better!

integer(0)

> for (i in seq_len(n)) print(i)

> for (i in seq(n)) print(i)

[1] 1
[1] 0

4.4 Conditional Code

It’s extremely common to need code to do different things depending upon

the number given to it. Let’s write a short function to find the absolute
value of a number.

> abs2 = function(x) {

+ if (x < 0) out = -x
+ else out = x

+ out
+ }
> abs2(-4)

[1] 4

The if() function will only execute the code which follows if the expression
in parentheses evaluates to TRUE 3 . When the expression is FALSE the code
3
Actually, any non-zero number will act the same way as TRUE, but it’s safer to only
use logicals.

34

41
following the else statement will be used instead. There is no need to
include an else, in which case the program will do nothing if the condition
in FALSE.
Take care not to allow the logical expression following the if() to be a
vector, or R will spit out a warning.

> abs2(c(1,-3)) # won't work properly.

Warning: the condition has length > 1 and only the first element
will be used

[1] 1 -3

If you only need to assign a single value (whether a number, logical, or

string) based on a condition, you can use ifelse():

> ifelse(TRUE, 94, "hello")

[1] 94

> ifelse(FALSE, 94, "hello")

[1] "hello"

4.5 while() loops

Here is a short function to check whether an integer is prime.

> isPrime = function(n) {

+ i = 2
+ if (n < 2) return(FALSE)

+ while (i < sqrt(n)) {

+ if (n %% i == 0) return(FALSE)
+ i = i+1
+ }
+ return(TRUE)
+ }
> isPrime(10)

35

42
[1] FALSE

> isPrime(37)

[1] TRUE

This illustrates several points. First, we don’t need to wait until reaching
the end of a function to return a value; we can use the return keyword
instead.
The other feature is the while() loop. This will keep running until the
expression in the parenthesis becomes false.

4.6 Writing for Speed

Loops are slow in R! It is usually much better to try to ‘vectorize’ a

function rather than calling it lots of times with a loop. Of course, for
squaring this is already built in to R.

> [Link](for (i in 1:1e6) i^2)

user system elapsed

0.177 0.012 0.196

> [Link](seq_len(1e6)^2)

user system elapsed

0.025 0.001 0.026

We can write a second function to do matrix-vector multiplication (see Ex-

ercise 4.4), but this time replacing the inner loop by a vectorized function
to take dot products.

> mult2 = function(A, b) {

+ n1 = nrow(A)
+ n2 = ncol(A)

+ out = numeric(n1)

+ for (i in 1:n1) {
+ out[i] = sum(A[i,]*b)

36

43
+ }
+ out
+ }

Now suppose we create a large matrix, and look at the difference in timing.

> A = matrix(rnorm(1e6), 1e3, 1e3)

> b = rnorm(1e3)
> [Link](mult(A,b))

user system elapsed

1.498 0.039 1.572

> [Link](mult2(A,b))

user system elapsed

0.013 0.005 0.018

> [Link](colSums(t(A)*b)) # can you see why this works?

user system elapsed

0.012 0.000 0.012

> [Link](A %*% b)

user system elapsed

0.004 0.000 0.004

The difference is dramatic. The moral of this is that it’s usually better to
use a built-in function, and almost always better to vectorize. The reason
%*% is so fast is that R calls underlying FORTRAN routines which have been
optimized over decades.

4.7 Recursion

Functions can recurse, which means they call themselves; here is a function
which calculates the entry Fn in the Fibonacci sequence with F0 = F1 = 1,
and Fk = Fk−1 + Fk−2 for k ≥ 2:

37

44
> fib = function(n) {
+ if (n < 2) return(1)
+ else return(fib(n-1) + fib(n-2))
+ }
> fib(10)

[1] 89

Recursion can be very slow though, so try to avoid it if possible. You can
also use Recall() instead of writing the function’s name in order to recurse.
Exercise 4.5. The number of moves required to complete a Towers of Hanoi
puzzle with k pieces is Hk = 2Hk−1 + 1 if k > 1, with H1 = 1. Write a
recursive function to evaluate Hk .
Exercise 4.6. Write a function to calculate a Fibonacci sequence using a
loop instead of a recursion. Compare the execution time to fib() (above)
for calculating the 30th term in the sequence using [Link]().

4.8 Scope

When a function is called, the code inside it is run in a separate environ-

ment to the code you run directly at the command line. This means it’s
possible for a variable inside a function to have the same name as something
at the command line without causing any problems:

> x <- 3
> f = function(y) {
+ x <- 5
+ x + y
+ }
> f(4)

[1] 9

> x # still the same value as before

[1] 3

However, if a function fails to find a variable withing its own environment,

then it will look to the parent environment for such a value: this is either
the function which called the current function, or the global environment
(i.e. the one you use at the command line).

38

45
> x <- 3
> g = function(y) {
+ x + y
+ }
> g(4)

[1] 7

Whilst this sort of behaviour can sometimes seem helpful, it is much better
to avoid writing confusing code like this. You are strongly recommended to
write functions which only require the information in their own arguments
to run.
This is the same principle used by the functions with() and subset():
they create an environment for the data frame (or list) you give as their
first argument; if any names supplied don’t match columns within the data
frame, R searches in the global environment:

> conv = 1.609

> with(hills, mean(dist/conv)) # where do conv and dist come from?

[1] 4.679

Exercise 4.7. What will happen if I create an object called dist before
running the commands above?

39

46
5 Graphics

Graphics and ‘data visualization’ are an integral part of statistics, and R

makes it easy to produce common plots quickly, as well as giving a power-
ful interface for more esoteric output. The basic command is the generic
function plot(). This will try to do the most sensible thing for the kind of
data you provide.
We have already seen that plot(x) and plot(x,y) will produce different
plots depending upon the class of the inputs. This is very partially summa-
rized in the following table:

y
x (missing) numeric factor
numeric series plot scatter plot spine plot
factor bar chart box plots spine plot

In fact there are many more plotting methods, most of which you will rarely
use.

5.1 One dimensional plots

For graphical summaries of one dimensional data we have already seen box-
plots and (in the practical) a time series for random walks. Among the most
useful is the histogram:

> hist(nlschools$lang, breaks=25, col=2)

40

47
 150 Histogram of nlschools$lang
Frequency

100
50
0

10 20 30 40 50 60

nlschools$lang

Note that the optional argument breaks chooses (approximately) how many
bins the histogram should have, and col alters the colour of the bars. Of
course, all plots should have properly labelled axes and a title, which can
be easily added.

> hist(nlschools$lang, breaks=25, col=2, xlab="Score",

+ main="Language test scores of Dutch 8th grade pupils")

Even the simple plot command for a single numeric vector comes with a
large range of options.

> x = cumsum(rnorm(250))
> plot(x, type="l", col=3)

41

48
 0
−5
x

−10
−15
−20

0 50 100 150 200 250

Index

Try this with type="b" or type="h" and see what happens. You can only
find out about a few of the graphics options with the documentation for
plot(). Try looking at ?par to find the real detail.

5.2 Adding to Plots

Consider the following simple scatter plot, augmented with the line y = x.

> x = rnorm(300)
> y = x + rnorm(300)
> plot(x,y, pch=20, col=4, cex=0.5)
> abline(a=0, b=1, lty=4, lwd=1.5)

42

49
 4

y=x
line of best fit
2
0
y

−2
−4

−4 −3 −2 −1 0 1 2

What does cex=0.5 do?

The function abline() plots lines of the form y = a + bx where a and b are
specified. As you can see, the width and appearance of the line is adjusted
with the options lty (line type), lwd (line width) and col (colour).
You can also use abline() in conjunction with output from a simple linear
model. Adding the following gives the red line:

> abline(lm(y ~ x), col=2)

5.3 Legends

The legend() command can be used to provide additional information

about your plots. The basic syntax is

> legend(x=-4, y=4, legend=c("y=x","line of best fit"),

+ lty=c(4,1), lwd=c(1.5,1), col=1:2)

where (x, y) is the top-left hand corner of the box, legend is a character
vector of annotations, and the other options are used to describe what to
display. Many other options are available via the help file.

43

50
5.4 Formulae and Boxplots

R contains a special object type called a formula which can be used to

represent statistical models compactly. Try typing some algebraic expression
at the console separated by the tilde operator ~ (this is on the # key on a
UK keyboard, and left of 1 on a US keyboard).

> x ~ a + b*c

x ~ a + b * c

The formula object can be used to express relationships between variables

under the convention that the left-hand side ‘is modelled by’ the right-hand
side.
This can be used, for example, when producing plots:

> data(genotype)
> head(genotype)

Litter Mother Wt
1 A A 61.5
2 A A 68.2
3 A A 64.0
4 A A 65.0
5 A A 59.7
6 A B 55.0

> boxplot(Wt ~ Litter, data=genotype)

44

51
 70
65
60
55
50
45
40
35

A B I J

The function interprets the formula as requiring that the left-hand side be
summarized in a way which is broken down by the right. Note that Wt and
Litter are contained within genotype, and are not recognized at the con-
sole4 , but the argument data=genotype ensures that the boxplot() func-
tion knows where to look for genotype$Wt.

5.5 Lattice Graphics

The plots above are all found in the base package of R, which is to say that
they are all preloaded functions. A very popular and powerful extension
to R’s graphics capabilities is made using the package lattice. The range
of plots which can be produced even using lattice’s default methods is
staggering, and we will show only a few small examples here.
The basic command is xyplot(), whose first argument is usually a formula.

> library(lattice)
> head(crabs)

sp sex index FL RW CL CW BD
4
That is, they are not in the global environment

45

52
1 B M 1 8.1 6.7 16.1 19.0 7.0
2 B M 2 8.8 7.7 18.1 20.8 7.4
3 B M 3 9.2 7.8 19.0 22.4 7.7
4 B M 4 9.6 7.9 20.1 23.1 8.2
5 B M 5 9.8 8.0 20.3 23.0 8.2
6 B M 6 10.8 9.0 23.0 26.5 9.8

> form <- log(FL) ~ log(RW) | sp*sex

> form

log(FL) ~ log(RW) | sp * sex

> xyplot(form, data=crabs)

2.0 2.2 2.4 2.6 2.8 3.0

M M
B O
3.0
2.8
2.6
2.4
2.2
2.0
log(FL)

F F
B O
3.0
2.8
2.6
2.4
2.2
2.0

2.0 2.2 2.4 2.6 2.8 3.0

log(RW)

The formula form in this case has three parts. The left-hand side log(FL)
is to be plotted against the right log(FL); since both these variables are
continuous, we will obtain a scatter plot. The conditioning bar ‘|’ indicates
that we wish the information to be broken down by the third term, sp*sex
(i.e. by species and by sex). Hence lattice produces four separate scatter
plots, each with the same axes.

46

53
The most common use of the lattice package is to produce these trellis
plots for representing multivariate data. A few more examples you might
find useful:

> histogram(~ height | [Link], data=singer)

60 65 70 75

Soprano 2 Soprano 1
40
30
20
10
0
Tenor 1 Alto 2 Alto 1
Percent of Total

40
30
20
10
0
Bass 2 Bass 1 Tenor 2
40
30
20
10
0
60 65 70 75 60 65 70 75

height

Notice that the command is histogram(), not hist(), and the plotting
options are different.

> library(MASS)
> densityplot(galaxies)

47

54
 0.00015

0.00010
Density

0.00005

0.00000

10000 20000 30000

galaxies

5.6 Function Plots

When applied to a vectorized function of one argument, plot() will produce

a graph in the specified range.

> plot(sin, -2*pi, 2*pi)

48

55
 1.0
0.5
0.0
sin

−0.5
−1.0

−6 −4 −2 0 2 4 6

For non-vectorized functions, this has to be done manually:

> func = function(x, n=10) {

+ idx = 1:n
+ return(sum((-1)^(idx+1)*x^idx/idx))
+ }
>
> plot(function(x) log(1+x), -0.5, 1.2)
>
> for (i in 2:10) {
+ x <- seq(from=-0.5, to=1.2, length=1000)
+ y <- sapply(x, func, n=i)
+ points(x, y, type="l")
+ }

49

56
 0.5
function(x) log(1 + x)

0.0
−0.5

−0.5 0.0 0.5 1.0

Using the lattice package you can also produce surface plots for functions
with two arguments using wireframe plots.

> func = function(x,y) (x^3-x)*sin(x+y)

> xs = ys = seq(-5, 5, [Link]=100)
> out = outer(xs, ys, func)
> wireframe(out)

50

57
 out

column
row

See also contourplot() and levelplot().

5.7 Customized Plots

To draw a plot from scratch, use the [Link]() command:

> [Link](1328) # to get the same values as me

> x = rnorm(100) # generate data
>
> [Link]()
> [Link](xlim=c(-3,3), ylim=c(-0.1,0.5))
> axis(side=1, pos=-0.1)
> hist(x, breaks=15, add=TRUE, freq=FALSE, col=2)
> plot(dnorm, -3, 3, add=TRUE)
> points(x, rep(-0.05,100), pch="|")
> title(main="Normal random variables")

[Link]() is used to control the range of the plot, the axis() function
draws on the axes, and title() is used to annotate with text. Try each of
the above commands in turn to see what they do.

51

58
5.8 Exporting Plots

You will likely need to use R plots in LaTeX documents for your practicals
and projects. If you are using a GUI such as the default R interface on
OS X or Windows, then select the window and go to File > Save As. I
recommend saving plots in PDF format, as this makes it easiest to integrate
with a LaTeX document. Other interfaces such as RStudio make it similarly
easy to create plots.
You can also save plots from the command line. The way to do this is to tell
R to send your plot commands to a file, instead of to the screen. This means
you won’t be able to see your plot whilst you produce it (but presumably
you’ll have already checked what it looks like!) Type pdf("[Link]")
to start, then run your chosen plot commands. Then finish with [Link]()
to go back to the default state. For example:

> pdf("[Link]")
> plot(hills)
> [Link]()

52

59
6 The apply Family of Functions

Much coding involves the repeated application of the same function to sev-
eral different pieces of data in a vector or list. For this reason, R has a series
of functions for performing such tasks, which results in much simpler and
easier to understand code.

6.1 apply()

If we want to perform a function on every row or column in a matrix (or on

an array, see next section), we can use the apply() function. The syntax
is apply(x, d, f), which will apply the function f to each of the dth di-
mension of x. If d=1 this corresponds to rows, and d=2 to the columns of a
matrix.

> A = cbind(1:10, (1:10)^2, (1:10)^3)

> apply(A, 2, sum)

[1] 55 385 3025

It will also work for ‘matrix-like’ objects, such as data frames (although see
also sapply() below).

> library(MASS)
> apply(hills, 2, mean)

dist climb time

7.529 1815.314 57.876

> apply(hills, 2, sd)

dist climb time

5.524 1619.151 50.041

Exercise 6.1. Using apply(), write a function which, given an (I × J)-

matrix X = (xij ) computes the magnitude of each row, that is
q
x2i1 + x2i2 + · · · + x2iJ , for each i = 1, . . . , I

and returns the results as a vector.

53

60
Exercise 6.2. Write a function which, given an (I × J)-matrix X, returns
a vector of length I with entries
x̄i+
yi =
si
where x̄i+ and si are respectively the sample mean and sample standard
deviation of entries in the ith row of X. [The mean divided by the standard
deviation is sometimes called the coefficient of variation.]
What happens if you use apply() with a function like range(), which re-
turns more than one value?
Exercise 6.3. Take a look at the data set EuStockMarkets (this is in the
datasets package, which should be already loaded). Find the mean absolute
change in returns from one day to the next for each stock (that is, the average
of |xi+1 − xi | over all days i). [Hint: recall the diff() function.]
Bonus*: Think of a more sensible measure of the volatility than this and
implement it [hint: one that doesn’t depend upon the scale].

Note that apply() does not run substantially faster than writing a loop to
do the same thing, it is simply easier to code up and to read.
For the particular task of sums or means of rows or columns in a matrix, R
contains special functions rowSums(), colSums(), rowMeans(), colMeans().
These are all much faster than the equivalent apply() commands.

> # 2000 x 2000 random matrix

> x = matrix(rnorm(4e6), 2000, 2000)
> [Link](apply(x,1,sum))

user system elapsed

0.121 0.004 0.127

> [Link](rowSums(x))

user system elapsed

0.01 0.00 0.01

Exercise 6.4. Write a function to renormalize the columns of a matrix so

that they sum to 1.
Exercise 6.5. Write a function to perform the same task as in Exercise 6.1,
but this time using rowSums(). Compare the speed of these two methods
using [Link]().

54

61
6.2 sapply() and lapply()

If we want to apply a function to every entry in a list or vector, we can

use lapply(). The syntax is just lapply(x, f) for a vector or list x and a
function f. It returns a list of the same length as x containing the results.
The following example uses lapply() twice, first on a vector, and then on
the resulting list.

> mu = c(-2,-1,0,1,2)
> out = lapply(mu, function(x) rnorm(100, mean=x))
> lapply(out, mean)

[[1]]
[1] -2.169

[[2]]
[1] -1.141

[[3]]
[1] 0.1753

[[4]]
[1] 1.105

[[5]]
[1] 2.008

Note that in the previous example, we didn’t want the final answer as a list,
so we might use the unlist() command to turn the results into a vector.
The sapply() function does this automatically when appropriate, but is
otherwise the same as lapply().

> sapply(out, mean)

[1] -2.1688 -1.1414 0.1753 1.1050 2.0080

If the function being applied returns something more complicated than a

single number, you should use lapply() instead.

6.3 replicate()

Sometimes we wish to repeat exactly the same operation multiple times,

without having a different input (this is typically used in conjunction with

55

62
the generation of random numbers).

> out = replicate(20, rnorm(100), simplify=FALSE)

Now out is a list of 20 independent data sets, each consisting of 100 standard
normal random variables.
Exercise 6.6. Suppose we wish to investigate the distribution of the maxi-
mum of 10 Poisson random variables with parameter λ = 5. Generate 1000
independent data sets consisting of such Poisson random variables (see the
command rpois()), find the maximum of each, and plot as a histogram.

6.4 tapply()

In the same vein as the previous examples, we may wish to evaluate a

function on data in a vector or list which are ‘grouped’ according to the
levels of some other factor. For example,

> library(MASS)
> head(genotype)

Litter Mother Wt
1 A A 61.5
2 A A 68.2
3 A A 64.0
4 A A 65.0
5 A A 59.7
6 A B 55.0

> with(genotype, tapply(Wt, Mother, mean))

A B I J
55.40 58.70 53.36 48.68

returns a vector of means. If the function provided gives more than a single
value, then tapply() will adapt accordingly:

> with(genotype, tapply(Wt, Mother, summary))

$A
Min. 1st Qu. Median Mean 3rd Qu. Max.
36.3 49.0 58.2 55.4 62.1 68.2

56

63
$B
Min. 1st Qu. Median Mean 3rd Qu. Max.
42.0 55.2 59.8 58.7 63.5 69.8

$I
Min. 1st Qu. Median Mean 3rd Qu. Max.
39.7 48.9 54.2 53.4 57.5 61.8

$J
Min. 1st Qu. Median Mean 3rd Qu. Max.
39.6 42.9 50.0 48.7 53.5 61.0

It is also possible to provide more than one grouping in the form of a list or
data frame, in which case the data are broken down by both:

> tapply(genotype$Wt, genotype[,1:2], mean)

Mother
Litter A B I J
A 63.68 52.40 54.12 48.96
B 52.33 60.64 53.92 45.90
I 47.10 64.37 51.60 49.43
J 54.35 56.10 54.53 49.06

Exercise 6.7. Find the heaviest rats born to each mother in the genotype()
data.

6.5 mapply()

Sometimes it may be useful to apply a function of several arguments repeat-

edly, where more than one argument can change.

> mapply(seq, from=c(1,4,-3), to=c(2,9,0), by=0.5)

[[1]]
[1] 1.0 1.5 2.0

[[2]]
[1] 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0

[[3]]
[1] -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0

64
57

64
`

Simple Linear Regression Analysis with the help of R

**What is the relationship between cows’ food intake and milk yield?

65
`

R- code for Parameter Estimation and plot

x=c(4,6,10,12)
y=c(3,5.5,6.5,9)
plot(x,y)

x_sq=x^2
y_sq=y^2
xy=x*y
t1=sum(x)
t2=sum(y)
t3=sum(x_sq)
t4=sum(y_sq)
t5=sum(xy)
t1
[1] 32
t2
[1] 24
t3
[1] 296
t4
[1] 162.5
t5
[1] 218
t6=t5-(t1*t2)/4
t7=t3-(t1^2)/4
beta1=t6/t7
beta1
[1] 0.65
t8=mean(x)
t8
[1] 8
t9=mean(y)
t9
[1] 6
beta0=t9-beta1*t8
beta0
[1] 0.8

hence our required regression model is given as

y_hat = 0.80 + 0.65*x

66
`

R code regression Analysis:

x=c(4,6,10,12)
y=c(3,5.5,6.5,9)
model=lm(y~x)
model

Call:
lm(formula = y ~ x)

Coefficients:
(Intercept) x
0.80 0.65

summary(model)

Call:
lm(formula = y ~ x)

Residuals:
1 2 3 4
-0.4 0.8 -0.8 0.4

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.8000 1.2166 0.658 0.5784
x 0.6500 0.1414 4.596 0.0442 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.8944 on 2 degrees of freedom
Multiple R-squared: 0.9135, Adjusted R-squared: 0.8703
F-statistic: 21.12 on 1 and 2 DF, p-value: 0.04422

a=[Link](x=14)
predict(model,a)
[1] 9.9

Predicted value of y at x=14 is 9.9

anova(model)
Analysis of Variance Table

Response: y
Df Sum Sq Mean Sq F value Pr(>F)
x 1 16.9 16.9 21.125 0.04422 *
Residuals 2 1.6 0.8
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Results:
1. Slope(𝛽1 ) milk yield (y) is expected to increase by 0.65lb
for each 1lb increases in food intake(x).
2. Y-Intercept(𝛽0 ) average milk yield(y) is expected to be
0.8lb when food intake (x) is 0.

67