Books and online help

• Introductory Statistics with R (Springer, Dalgaard).

• A first course in statistical programming with R (CUP, Braun and
Murdoch).
Scientific Programming with R • Computational Genome Analysis: An Introduction (Springer, Deonier,
Tavaré and Waterman).
• S programming (Springer, Venables and Ripley).
Stephen Eglen
• R programming for Bioinformatics (CRC Press, Gentleman).
Laurent Gatto
• Scientific programming and simulation using R (CRC, Jones, Maillardet
and Robinson).
September 6, 2015 • The art of R programming (No Starch Press, Matloff).
• Writing Scientific Software (WSS) (CUP, Oliveira and Stewart).
• [Link], [Link], [Link]
• R -help mailing list.

• Eglen (2009) [Link]

1 / 191 2 / 191

Aims of course Part 2: Scientific computing issues

This course aims to teach R as a general-purpose programming language.
Issues specific to Computational Biology (e.g. Bioconductor packages) are In part 2 of the course1 , we will explore various other topics, building on
covered in other course modules. core knowledge of R .
In part 1, topics to be mastered in this course include: • Numerical integration
• Interactive use of R . • Phase plane analysis
• Basic data types: vector, matrix, list, [Link], factor, • Handling large files/data bases
character. • String processing (e.g. for genomic data)
• Writing scripts. • Advanced graphing / presenting results
• Graphical facilities. • Reproducible research
• Writing your own functions. • Future directions (R and generally)
• File input/output.
• (Object-oriented programming)
• Control-flow statements, looping.
• (Package development)
• Vectorization.
• (R profiling)
• Numerics issues.
• Debugging. 1
tentative
3 / 191 4 / 191
Outline
Introduction
Scientific programming/software
Vectors
Calling functions • Different from software engineering (but should try to adhere to SE
Scripts best practice, of course).
Matrices • Moving target.
Boolean logic • Domain scientists write the code (some argue this is a weakness).
Lists
• Has of course to be accurate, user-friendly (no GUI vs CLI ranting
Factors
here), usable and useful, flexible, efficient and open, owned by the
Character arrays
community, facilitate reproducible research.
Objects in your environment
• Contribute to users education (i.e not be a black box), in terms data
Basic plotting
requirements, the data processing and result interpret.
Reading/writing data to file system
Importance of documentation.
Writing functions
Conditionals and looping See Gentleman et al. (2004) Genome Biology (Bioconductor paper) for an
Vectorization example of successful scientific software.
Random number generation
Debugging
5 / 191 6 / 191

What is R ?
Trustworthy software
The complexity of the data processes and of the computations applied to
them mean that those who receive the results of modern data analysis have
limited opportunity to verify the results by direct observation. Users of the
• Statistical computing environment, and programming language.
analysis have no option but to trust the analysis, and by extension the
software that produced it. Both the data analyst and the software provider • Very popular in many areas of statistics, computational biology.
therefore have a strong responsibility to produce a result that is trustworthy, • “Programming with data” (Chambers)
and, if possible, one that can be shown to be trustworthy. • Approach: command-line for one-liners; interactive usage; write
scripts/functions for larger work (edit/run cycle); develop package for
This places an obligation on all creators of software to program in such a consolidation and distribution.
way that the computations can be understood and trusted.

John M. Chambers, Software for Data Analysis (Springer)

7 / 191 8 / 191
History Strengths of R

• GPL’d, available on many platforms.

• S language came from Bell Labs (Becker, Chambers and Wilks). • Excellent development team with yearly release cycle.
Commercial version S-plus (1988). • Source always available to examine/edit.
• R emerged as a combination of S and Scheme: Ross Ihaka and Robert • Fast for vectorized calculations.
Gentleman (NZ). • Foreign-language interface (C/Fortran) when speed crucial, or for
• 1993: first announcement. interfacing with existing code.
• 1995: 0.60 release, now under GPL. • Good collection of numerical/statistical routines.
• 2014-10-07: release 3.1.1. Stable, multi-platform. Major release every • Comprehensive R Archive Network (CRAN) ∼ 5913 packages
March/April. [2014-10-07] (cf 1000 in April 2007).
• R-core now 20 people, key academics in field, including John Chambers. • On-line doc, with examples.
• High-quality graphics (pdf, postscript, quartz, x11, bitmaps). Often
used just for plotting . . .

9 / 191 10 / 191

Graphics example Weaknesses of R

• Loops are slow. Learn how to vectorize solutions.

• No fast compiler yet, and unlikely to happen due to nature of language.
Byte compiler available in compiler package.
• No (decent) endorsed GUI built-in to R . Tk is available within base R ,
and packages for other graphical tooklits (e.g. Gtk2, Qt) are also
available.
“Programming Graphical User Interfaces with R ”, M. F. Lawrence and
J. Verzani

Jean YH Yang; gpQuality

[Link]
11 / 191 12 / 191
Brief comparison to matlab Using R

• Flexible language, similar to matlab, but definitely not “everything is a • Start-up: type ‘R’ at command line.
matrix”. Frames, lists, vectors . . . • Type commands interactively, and get results.
• From matlab to R : • Type commands into a file; source( ' myfile .R'); edit file . . .
[Link] • Mac/Win has a GUI for interactive use, with internal editors.
• Comprehensive matlab and R guide: http: • All platforms have a command-line interface
//[Link]/faculty/hiebeler/comp/[Link] • Many external editors have support for R , including
• Use x[i] not x(i) for indexing vectors. • Emacs through ESS ([Link]
• Making vectors: x <- c(10, 9, 5, 1) • Eclipse IDE ([Link]
• Rstudio ([Link]
• Assignment: best to use <- rather than =.
• ...

13 / 191 14 / 191

My very first R session Interacting with R

x <- rnorm(50, mean=4)

x • Can use up/down arrow keys to go through command history. Within a
mean(x) command, use left/right arrow keys to edit.
range(x) • History can be saved over sessions (?history).
hist(x) • Multiple commands can be put onto one line, using ; as separator
## check help -- how to change title? between lines, e.g. x<-10; y<-3; a <- 5.
?hist • TAB can do object/file completion.
hist(x, main="my first plot")
q()

15 / 191 16 / 191
Objects and Functions Objects and Functions

> x <- 200

> half.x <- x/2
R manipulates objects. Each object has a name and a type (vector, > threshold <- 95.0
matrix, list, . . . ) > age <- c(15, 19, 30)
Name of an object: letters (upper/lower case are distinct), digits, period. > age[2] ## [] for accessing element.
Start with a letter.
Objects set by way of assignement. Use the <- assignment operator rather [1] 19
than = wherever possible. (Does i = i+1 make sense?)
> length(age) ## () for calling function.

[1] 3

17 / 191 18 / 191

What’s up with the assignment and underscore? (Advanced) Outline

Introduction
Vectors
Historically, underscore was used in S for assignment (because an old system Calling functions
keyboard had a key equivalent to the ASCII underscore that generated a Scripts
back arrow). Hence underscore was not used within variables.
Matrices
More recently, = is now available as an assignment operator (similar to
Boolean logic
languages like C), but is frowned upon as it can be confusing.
Lists
What does i = i+1 imply mathematically?
Factors
Better to stick to i <- i + 1 and use equals just within calls to functions,
Character arrays
e.g. runif(max=3).
Note also that assignments return values: Objects in your environment
Basic plotting
> y <- 1 + ( x <- 9 ) Reading/writing data to file system
> a <- b <- 0 Writing functions
Conditionals and looping
[Link] Vectorization
Random number generation
Debugging
19 / 191 20 / 191
Vectors Vectors
Some operations work element by element, others on the whole vector,
Vectors are a fundamental object for R . Scalars are treated as vector of
compare:
length 1.
> y <- c(20, 49, 16, 60, 100)
> y <- c(10, 20, 40) > min(y)
> y[2]
[1] 16
[1] 20
> range(y)
> length(y)
[1] 16 100
[1] 3
> sqrt(y)
> x <- 5
> length(x) [1] 4.472136 7.000000 4.000000 7.745967 10.000000

[1] 1 > log(y)

21 / 191
[1] 2.995732 3.891820 2.772589 4.094345 4.605170 22 / 191

Generating vectors Accessing and setting elements

> x <- seq(from=100, by=1, length=20)
> x[3] ## just element 3.
Many short hand methods for regular sequences; c() for irregular.
[1] 102
> x <- seq(from=1, to=9, by=2)
> y <- seq(from=2, by=7, length=3) > x[c(12,14)] ## element 12 and 14
> z <- 4:8
> a <- [Link](5) ## fast for integers [1] 111 113
> b <- c(3, 9, 2)
> x[1:5]
> d <- c(a, 10, b)
> e <- rep( c(1,2), 3)
[1] 100 101 102 103 104
> f <- integer(7)
> bad <- 1:4
> x[-bad] ## exclude elements

[1] 104 105 106 107 108 109 110 111 112 113 114 115 116 117
23 / 191
[15] 118 119 24 / 191
Accessing and setting elements Accessing and setting elements

Can also provide a logical vector of same length as vector (logical values
explained later). Elements can be set in several ways

> x <- c(5, 2, 9, 4) > x <- rep(0,10)

> v <- c(TRUE, FALSE, FALSE, TRUE) > x[1:3] <- 2
> x[v] > x[5:6] <- c(-5, NA)
> x[7:10] <- c(1,9) ## recycling.
[1] 5 4

25 / 191 26 / 191

Recycling rule (Advanced) Recycling rule (Advanced)

Recycling is convenient, but dangerous; when vectors are of different lengths,

the shorter one is often recycled to make a vector of the same length.

> a <- c(1,5) + 2

> x <- c(1,2); y <- c(5,3,9,2) > x <- 1:10
> x + y > y <- x * 2
> z <- x^2
[1] 6 5 10 4 > y + z
> x + 1:2
> x + c(y,1) ## odd recycling, warning. > x + 1:3

Warning in x + c(y, 1): longer object length is not a

multiple of shorter object length

[1] 6 5 10 4 2

27 / 191 28 / 191
Recycling rule (Advanced) Naming indexes of a vector
> x <- 1:10
> y <- x * 2 > joe <- c(24, 1.70)
> z <- x^2 > joe
> y + z
[1] 24.0 1.7
[1] 3 8 15 24 35 48 63 80 99 120
> names(joe)
> x + 1:2
NULL
[1] 2 4 4 6 6 8 8 10 10 12
> names(joe) <- c("age", "height") ## replacement function
> x + 1:3 > joe

Warning in x + 1:3: longer object length is not a multiple age height

of shorter object length 24.0 1.7
[1] 2 4 6 5 7 9 8 10 12 11
29 / 191 30 / 191

Naming indexes of a vector Common functions for vectors

• length()
> joe["height"] == joe[2]
• rev()
height • sum(), cumsum(), prod(), cumprod()
TRUE • mean(), sd(), var(), median()
• min(), max(), range(), summary()
Refering to index by name rather than by position can make code more
readable, and flexible. Cannot do things like x[1:4] easily though, since • exp(), log(), sin(), cos(), tan() (radians, not degrees)
you need to name all four elements you want. • round(), ceiling(), floor(), signif()
Although extremly useful, names have a cost when processing large objects. • sort(), order(), rank()
• which(), [Link]()
Note: in second use of names() above, we are actually using the
• any(), all()
replacement function names<-, see later.

31 / 191 32 / 191
Functions as function args Outline
Introduction
Vectors
Calling functions
Functions can be called within function calls; the following are equivalent: Scripts
Matrices
> x <- c(3, 2, 9, 4) Boolean logic
> Lists
> y <- exp(x); z1 <- which(y > 20) ## case 1 Factors
> z2 <- which ( exp(x) > 20) ## case 2 Character arrays
> Objects in your environment
> [Link](z1, z2) Basic plotting
Reading/writing data to file system
[1] TRUE
Writing functions
Conditionals and looping
Vectorization
Random number generation
Debugging
33 / 191 34 / 191

Default values for function arguments Default values for function arguments

> x <- c(2.091, 4.126, 7.925)

> round() ## required arg is missing
A function will error if not all required arguments are provided. Some
functions have both required and optional arguments. If the optional Error in eval(expr, envir, enclos): 0 arguments passed to
arguments are not provided, they are either ignored, or they take a default ’round’which requires 1 or 2 arguments
value.
> round(x)
Usage:
round(x, digits = 0) [1] 2 4 8

> round(x, digits = 2)

[1] 2.09 4.13 7.92

Let’s see how this works in more detail.

35 / 191 36 / 191
Argument matching Argument matching
Typical calls are as follows:
> seq(1, 3, 0.5) ## positional matching
R has a flexible method for specifying arguments to function. We can either
provide an actual value for a formal argument, or give arguments as [1] 1.0 1.5 2.0 2.5 3.0
key=value (or formal=actual).
> seq(1, 5,[Link]=3) ## can skip args (e.g. by)
As an example, let’s look at help for seq:
[1] 1 3 5
seq(from = 1, to = 1, by = ((to - from)/([Link] - 1)),
[Link] = NULL, [Link] = NULL, ...) > seq(to=5) ## order not important.

[1] 1 2 3 4 5
(NB: in seq(from=x), from is the formal argument of the function, and
here x is the actual value.) > seq(f=5,t=1) ## abbrev tags.

The ... notation allows for other arguments to be passed, which are not [1] 5 4 3 2 1

used by this function. > seq(len=5, 1,2) ## tags removed before positional matching

[1] 1.00 1.25 1.50 1.75 2.00

37 / 191 38 / 191

. . . in function calls (Advanced) Replacement functions (Advanced)

Why do some functions, like sqrt, require only one argument, yet others
take many arguments?
Functions like c, cbind, have ... in the arguments: > x <- 1:5
Usage: > x

c(..., recursive=FALSE) [1] 1 2 3 4 5

Arguments: > length(x)

[1] 5
...: objects to be concatenated.
The ... indicate any number of objects may be passed, not just (say) one > length(x) <- 2
or two. > x
The result of c() is to combine them all into one long vector, taking into
account if the keyword “recursive” is provided (when args are first flattened). [1] 1 2
The ... can also indicate that other arguments can be provided which are
not processed directly by this function, but may be useful for other functions
(e.g. popular when plotting).
39 / 191 40 / 191
Replacement functions (Advanced) Getting help: key commands

• help(hist) to see help file (or ?hist).

Normally length(x) would return a value, rather than you assigning a
• args(hist) to see arguments of a function.
value to the function! These are replacement functions, see help page:
• example(boxplot) run examples in help page.
Usage: • options(help type="html") will then use web-browser for help.
• [Link]("histogram")
length(x)
• demo() to list all demos, e.g. demo(graphics)
length(x) <- value
NB: In the R terminal ?command works as shorthand for help("command")
except for a small number of commands, e.g. if, while. Use the longhand
for these.

41 / 191 42 / 191

Help pages Numbers and special values

• numeric (floating-point, double): 12, 4.92, 1.5e3 – [Link]()
(integers converted to f.p.)
• integers 1L – [Link]()
• What you can expect to find: • complex: 3+2i – [Link]()

• Description – one line summary > typeof(1)

• Usage – formal arguments
• Arguments – interpretation of arguments [1] "double"
• Details – what the function does
> typeof(1L)
• Value – return value.
• References – documentation [1] "integer"
• See also – helps you find related pages
• Examples – guaranteed to run: example(hist) > [Link](1)

[1] FALSE

> [Link](1L)

[1] TRUE

43 / 191 44 / 191
Numbers and special values Numbers and special values

> typeof(NA)

Special values: [1] "logical"

• NA: not available. (Often used to represent missing data point) –
> typeof(NaN)
[Link]()
• NaN: not a number. e.g. 0/0 – [Link]() [1] "double"
• Inf, -Inf: ±∞ – [Link]()
> typeof(Inf)
You will also meet:
• NULL: often, list of zero length – [Link]() [1] "double"

> typeof(NULL)

[1] "NULL"

45 / 191 46 / 191

Operator precedence ?Syntax Operator precedence ?Syntax

> 3 * 4 + 2 != 3 * (4 + 2)

[1] TRUE
> 3 * 4 + 2 != 3 * (4 + 2)
> 2^3+1 != 2^(3+1) > 2^3+1 != 2^(3+1)
> 1:5-1
[1] TRUE

> 1:5-1

[1] 0 1 2 3 4

47 / 191 48 / 191
Operator precedence ?Syntax Operators
Subset taken from ?Syntax, see that page for full list. Highest precedence
at top.
'[ [[' indexing Most operators will be familiar, but some may not:
'$ @' component / slot extraction
'^' exponentiation (right to left) x <- 10
'- +' unary minus and plus x == 4 ## test for equality
':' sequence operator x != 10 ## not equal?
'%any%' special operators 7 %/% 2 ## division, ignoring remainder. (3)
'* /' multiply, divide 7 %% 2 ## remainder (1)
'+ -' (binary) add, subtract x <- 9 ## assignment
'< > <= >= == !=' ordering and comparison x <<- 9 ## assign x to 9 in the global env. (BAD)
'!' negation ## Raising to a power can be done in two ways.
'& &&' and [Link]( 10.1 ** 2.5, 10.1^2.5 )
'| ||' or
'<- <<-' assignment (right to left)
'?' help (unary and binary)
Bottom line: use parentheses to order preference.
49 / 191 50 / 191

When things go wrong Types of parentheses

Syntax errors are those where you’ve just made a typing mistake.
Logical errors are harder to find!

Common problems: • f(3,4) – call the function f, with arg1=3, arg2=4.

• missing close bracket leads to continuation line. • a + (b*c) – use to enforce order over which statements are executed.
> x <− ( 1 + ( 2 ∗ 3 ) • { expr1; expr2; ...; exprn } – group a set of expressions into
+ one compound expression. Value returned is value of last expression;
Hit Ctrl C (below) or keep typing! used in looping/conditionals.
• too many parens: • x[4] – get the 4th element of the vector x.
2 + (2 ∗ 3)) • l[[3]] – get the 3rd element of some list l, and return it. (compare
• wrong/mismatched brackets (see next slide). with l[3] which returns a list with just the 3rd element inside – will
• Likewise, do not mix double quotes and single quotes. see list objects later).
(Unless you need to quote within quotes.)
• ...
• wrong variable name (not syntax error)
• When things seem to take too long, try C-c [Ctrl and C, together]
51 / 191 52 / 191
From interactive to source files Outline
Introduction
Vectors
• Typing in commands interactively is good for one-liners, but soon you
Calling functions
will want to switch to putting your sequence of commands into a script
file, and then ask R to run (source) those commands. Scripts
Matrices
• This leaves to a rapid edit–run–edit cycle.
Boolean logic
• e.g. type these commands into a file trig.R: Lists
x <- seq(from=0, to=2*pi, length=100) Factors
y <- sin(x) Character arrays
z <- cos(2*x) Objects in your environment
z ## will not appear when source'd Basic plotting
print(y[1:10]) ## should use print() Reading/writing data to file system
plot(x, y, type='l') Writing functions
lines(x, z, type='l', col='red') Conditionals and looping
Vectorization
• Eval within R using source(’trig.R’).
Random number generation
Debugging
53 / 191 54 / 191

Scripts Why are scripts a good thing?

• Use source(’trig.R’, echo=TRUE) to see commands and output.
Or use print(x) to print an object within a script.
• Keep your code open in the editor in one window, and keep R running
in another window. • You don’t have to remember what commands you ran, they are saved
• Are you in the right directory? Check that you can see your script file in the file.
in the same directory as where R is currently. Check dir(), and setwd, • This corresponds to the “source is real” philosophy of using S and R .
see later. • You can easily give your work to others, by passing them the file.
• On unix, the initial directory is the directory from where you started R .
• You can eventually run your scripts in BATCH, i.e. non-interactively.
On windows, the initial directory might be “My Documents”. You may Good for long jobs which you can leave overnight.
need to change directory (setwd) first.
• Use a good editor that helps you spot mistakes (e.g. paren matching,
syntax highlighting). Examples: Emacs/ESS, gedit, Rstudio.
• Use .R or .r as the filename suffix. Avoid any temptation to put
spaces (although R does not mind) in your filenames!

55 / 191 56 / 191
Running scripts in batch (Advanced) Rscript
• Rscript works as an interpreter. It is quicker to start than R (fewer
packages are loaded?).
• At the command line, type R CMD BATCH trig.R. R will start up,
• Can use interactively or for standalone scripts.
process your commands and then quit.
• Output is stored in the file [Link] $ Rscript -e 'round(runif(10))'
• If there were no errors, the last line of the output file shows the time
[1] 0 1 0 0 1 0 1 1 0 0
taken to run the script.
• Any output is not shown on the screen but sent to a PDF called #!/usr/bin/env Rscript
[Link]. args <- commandArgs(TRUE)
• This is a GREAT way of testing your scripts, since R starts with an
stopifnot(length(args)==3)
empty workspace, you will see if you have all the steps needed. args = [Link](args)
• Aim to always leave your scripts in a working state at the end of a n = args[1]; mean = args[2]; sd = args[3]
session, so that a few days later you don’t have to remember why it rnorm(n, mean, sd)
wasn’t working!
$ ./simple_rnorm.R 5 10 2
[1] 8.651807 13.372094 10.063347 10.703155 6.351886
57 / 191 58 / 191

Commenting your work Line wrapping

• Do not be shy when putting comments into your code.
• Meaningful variable names help, but do document. At a bare minimum,
each file should state at the top what the purpose of the file. Important • Line-wrapping. Do not write beyond around column 72, for readability.
variables and functions should be clearly documented. You can break long expressions at suitable points.
• You may think it obvious how your code works, but try looking at it a
• End of line shold not look like end of an expression. Compare:
week or a month later and then see if you clearly understand it. If in
doubt, document it. x <- sqrt( c(100, 200, 300, 400, 500) ) + 10
• Describe what your code is doing, not how it is doing it (WSS, p79). x <- sqrt( c(100, 200, 300, 400, 500) )
Compare the following two: + 10
x <- sqrt( c(100, 200, 300, 400, 500) ) +
> s <- s + 1 ## prepare to process next subject
10
> j <- j + 1 ## increment j by 1.
x <- sqrt( c(100, 200, 300, 400, 500) ) +
• Comments can be put before commands, if you temporarily do not 10
want to run that command; remove the comments when you want to
run the command again, or delete the line.
> ## x <- c(x, c(1,2,3))
59 / 191 60 / 191
Line wrapping Line wrapping

> ## 1: ok - all fits onto one line, just. > ## 3: solved, by moving the operator (+) up.
> (x <- sqrt( c(100, 200, 300, 400, 500) ) + 10) > x <- sqrt( c(100, 200, 300, 400, 500) ) +
+ 10
[1] 20.00000 24.14214 27.32051 30.00000 32.36068
> x
> ## 2: not okay -- first line is seen as complete.
[1] 20.00000 24.14214 27.32051 30.00000 32.36068
> x <- sqrt( c(100, 200, 300, 400, 500) )
> + 10 > ## 4: as 3, but indentation makes it clearer.
> x <- sqrt( c(100, 200, 300, 400, 500) ) +
[1] 10
+ 10
> x
> x
[1] 20.00000 24.14214 27.32051 30.00000 32.36068
[1] 10.00000 14.14214 17.32051 20.00000 22.36068

61 / 191 62 / 191

Outline Matrices
Introduction A matrix is just a vector with some additional markup to reformat it. Matrix
Vectors stored in column-major order (like fortran, unlike C).
Calling functions > x <- 1:6
Scripts > [Link](x)

Matrices [1] FALSE

Boolean logic
Lists > dim(x) <- c(2,3)
> [Link](x)
Factors
Character arrays [1] TRUE
Objects in your environment
> x
Basic plotting
Reading/writing data to file system [,1] [,2] [,3]
Writing functions [1,] 1 3 5
[2,] 2 4 6
Conditionals and looping
Vectorization > dim(x)
Random number generation
[1] 2 3
Debugging
63 / 191 64 / 191
Matrices Matrices

> x[2,2] ## extracting a value.

[1] 4
> x[, 2, drop=TRUE] ## default
> x[1,] ## extracting row
[1] 3 4
[1] 1 3 5
> x[, 2, drop=FALSE] ## gotcha!
> x[1:2, 2:3]
[,1]
[,1] [,2] [1,] 3
[1,] 3 5 [2,] 4
[2,] 4 6

> x[,2] ## not column vector!

[1] 3 4

65 / 191 66 / 191

Typical matrix construction methods Typical matrix construction methods

• matrix()
• cbind() > matrix( floor(runif(6, max = 50)), nrow = 3 ) ## ncol = 2
• rbind()
[,1] [,2]
[1,] 43 37
> m <- matrix( floor(runif(6, max=50)), nrow=3) [2,] 6 8
> x <- rbind( c(1,4,9), c(2,6,8), c(3,2,1)) [3,] 47 48
> y <- cbind( c(1,2,3), 5, c(4,5,6))

67 / 191 68 / 191
Typical matrix construction methods Typical matrix construction methods

> rbind( c(1,4,9), c(2,6,8), c(3,2,1) ) > cbind( c(1,2,3), 5, c(4,5,6) ) ## recycling again

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

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

69 / 191 70 / 191

Typical matrix construction methods Common matrix operations

Note that matrix indices can also be named: • diagonal: diag(x) – watch if x matrix or scalar!
• matrix multiplication: %*% vs * (element-wise)
> dimnames(m) <- list(student = c("ann", "bob", "joe"),
+ exam = c("math", "french"))
> x <- matrix(1:4, 2,2)
> m
> i <- diag(2) ## 2 x 2 identity matrix
exam > x %*% i ## should be x
student math french
[,1] [,2]
ann 19 21
[1,] 1 3
bob 35 3
[2,] 2 4
joe 12 12
> x * i ## not x!
> m["bob", ] ## get bob's scores
[,1] [,2]
math french
[1,] 1 0
35 3
[2,] 0 4
71 / 191 72 / 191
Common matrix operations Arrays
Arrays as extension of matrices to multiple dimensions.
> array(1:12, c(2,2,3))
• transpose: t(x)
, , 1
• dim, row, ncol
[,1] [,2]
• inverse: solve(x), [1,] 1 3
[2,] 2 4

> (x %*% solve(x)) == diag(nrow(x)) , , 2

[,1] [,2]
[,1] [,2] [1,] 5 7
[2,] 6 8
[1,] TRUE TRUE
[2,] TRUE TRUE , , 3

[,1] [,2]
[1,] 9 11
[2,] 10 12

73 / 191 74 / 191

Outline Boolean values – ?logical

Introduction Logical values TRUE/FALSE (avoid abbrev to T/F).
Vectors
Calling functions TRUE/FALSE equivalent to 1/0
Scripts
Matrices > [Link](TRUE)
Boolean logic
[1] 1
Lists
Factors > [Link](FALSE)
Character arrays
Objects in your environment [1] 0
Basic plotting
> [Link](1)
Reading/writing data to file system
Writing functions [1] TRUE
Conditionals and looping
Vectorization > [Link](0)
Random number generation
Debugging [1] FALSE
75 / 191 76 / 191
Boolean values – ?logical Boolean values – ?logical

> d[which(d>5.0)]
> d <- c(3.2, 1.0, 4.0, 9.2, 2.3, 8.1, 6.3) [1] 9.2 8.1 6.3
> d > 5.0
> [Link] <- (d > 3.0) & (d< 5.0)
[1] FALSE FALSE FALSE TRUE FALSE TRUE TRUE > d[[Link]]
> d[d> 5.0] [1] 3.2 4.0
[1] 9.2 8.1 6.3 > d[![Link]]
> which(d>5.0) [1] 1.0 9.2 2.3 8.1 6.3
[1] 4 6 7 > ifelse(d > 3.0, 1.0, 0.0) ## same as [Link](d > 3)

[1] 1 0 1 1 0 1 1

77 / 191 78 / 191

Boolean values – ?logical Boolean logic: issues

Key operators for handling boolean values: a & b (same for a | b) is an elementwise operation, with a result the
same length as the longer of a, b (recycling is used if one vector is shorter).
> !TRUE ## negation: swap T -- F.
a && b examines only the first element of a and b, returning one logical
[1] FALSE value. Lazy evaluation is used: we calculate only what’s needed to
determine result.
> TRUE & FALSE ## and: both must be true.
> TRUE || [Link]()
[1] FALSE
[1] TRUE
> FALSE | TRUE ## or: one must be true.
> TRUE && stop("no")
[1] TRUE
Error in eval(expr, envir, enclos): no
> xor(TRUE, TRUE) ## xor: only one is true.
> FALSE && stop("no")
[1] FALSE
[1] FALSE
79 / 191 80 / 191
Lazy (Advanced) Boolean logic: issues
R uses Lazy evaluation, which delays the evaluation of an expression (here
the argument) until its value is actually required [2 ]:
Comparing numbers: When testing numbers for equality, can use x == y
> f <- function(x) { 10 } when x, y are integers, otherwise use [Link](x,y). See later on
> [Link](f([Link](3))) numerics.

user system elapsed Avoid using F

0 0 0
> F <- 3
> f <- function(x) { force(x); 10 } > F == FALSE
> [Link](f([Link](3)))
[1] FALSE
user system elapsed
0.000 0.000 3.001

2
example from Hadley Wickham’s devtools
81 / 191 82 / 191

Outline What is a list?

Introduction
Vectors A list is used to collect a group of objects of different sizes and types. Very
Calling functions flexible. Often returned as the result of a complex function (e.g. model fit)
Scripts to return all relevant information in one object.
Matrices
Boolean logic > l <- list(who='joe', height=1.70, dob=c(1960, 12, 1))
Lists > l
Factors
$who
Character arrays
[1] "joe"
Objects in your environment
Basic plotting
$height
Reading/writing data to file system
[1] 1.7
Writing functions
Conditionals and looping $dob
Vectorization [1] 1960 12 1
Random number generation
Debugging
83 / 191 84 / 191
What is a list? What is a list?

> length(l)

[1] 3 List elements can either be accessed by name (e.g. l$height or

l[[’height’]] – if named list) or by position (l[[2]]).
> names(l) ## show components
When using numbers to index list, compare l[2] (a list with one element)
[1] "who" "height" "dob" with l[[2]]. You can therefore do l[2:3] but not l[[2:3]].
> l$height ## access an element.

[1] 1.7

85 / 191 86 / 191

What is a list? Modifying lists (Advanced)

> unlist(l) ## opposite of list

We can append new items to list either by making a new list from the old
who height dob1 dob2 dob3 one (e.g. 1) , or directly by assigning new element (e.g. 2):
"joe" "1.7" "1960" "12" "1"
> l1 <- list(who="fred")
> str(l) ## structure of l
> l1 <- c(l1, height=1.8) ## e.g. 1
> l1[["dob"]] <- c(1965, 10, 17) ## e.g. 2
List of 3
$ who : chr "joe"
$ height: num 1.7
$ dob : num [1:3] 1960 12 1

87 / 191 88 / 191
Modifying lists (Advanced) Modifying lists (Advanced)

Deleting list items:

Finally, for completeness, here is a way to predefine a list of given length
> l1["height"] <- NULL and gradually fill it in:
> str(l1)
> empty <- vector("list", 3) ## Prealloc to given length.
List of 2 > names(empty) <- c("who", "height", "dob")
$ who: chr "fred" > empty[["height"]] <- 1.8
$ dob: num [1:3] 1965 10 17

89 / 191 90 / 191

Initialisation Data frames

Initialising R object (list, vectors, . . . ) is much faster than creating and
extending these objects on the fly
A data frame is like a matrix, but each column can be of a different type.
> n <- 1e4
> l <- list() Data frame is stored as a list, with each element a vector of same length.
> [Link](for (i in 1:n) l[[i]] <- rnorm(1e3) ) Useful for reading in tabular data from a file (see [Link]).
user system elapsed > nms <- c("joe", "fred", "harry")
2.728 0.045 2.776 > a <- c(24, 19, 30)
> ht <- c(1.7, 1.8, 1.75)
> l <- vector("list", n)
> s <- c(TRUE, FALSE, TRUE)
> [Link](for (i in 1:n) l[[i]] <- rnorm(1e3) )

user system elapsed

1.211 0.005 1.215

91 / 191 92 / 191
Data frames Data frames

> cbind(a, ht) ## matrix, preserves > cbind(nms, ht) ## matrix, ht col becomes strings.

a ht nms ht
[1,] 24 1.70 [1,] "joe" "1.7"
[2,] 19 1.80 [2,] "fred" "1.8"
[3,] 30 1.75 [3,] "harry" "1.75"

> class(cbind(a, ht)) > class(cbind(nms, ht))

[1] "matrix" [1] "matrix"

93 / 191 94 / 191

Data frames Data frames

> d$age ## same as d[, "age"]

> d <- [Link](name = nms, [1] 24 19 30

+ age = a,
+ height = ht, > names(d)
+ student = s)
> class(d) [1] "name" "age" "height" "student"

[1] "[Link]" > d[2,] ## access 2nd row.

name age height student

2 fred 19 1.8 FALSE

95 / 191 96 / 191
Data frames Outline
Compare how a data frame (d) is printed, compared to printing Introduction
[Link](d). Vectors
> d Calling functions
Scripts
name age height student
1 joe 24 1.70 TRUE Matrices
2 fred 19 1.80 FALSE Boolean logic
3 harry 30 1.75 TRUE
Lists
> [Link](d) Factors
Character arrays
$name
[1] joe fred harry Objects in your environment
Levels: fred harry joe Basic plotting
$age Reading/writing data to file system
[1] 24 19 30 Writing functions
$height Conditionals and looping
[1] 1.70 1.80 1.75 Vectorization
$student Random number generation
[1] TRUE FALSE TRUE Debugging
97 / 191 98 / 191

Factors (Advanced) Factors (Advanced)

> scores1 <- c('good', 'poor', 'bad', 'poor',

+ 'bad', 'bad', 'good')
> scores <- factor(scores1)
> scores
(Mostly seen when reading in data frames)
Factors internally code categorical variables with a number. e.g. 1=Sunday, [1] good poor bad poor bad bad good
2=Monday, . . . 7=Saturday. For large vectors, this is more efficient storage, Levels: bad good poor
especially when character strings repeat. Can also make code more readable. > levels(scores)

[1] "bad" "good" "poor"

Also useful in many statistical functions, using the formula interface.
> [Link](scores) ## integer representation

[1] 2 3 1 3 1 1 2

> [Link](scores) ## show strings

[1] "good" "poor" "bad" "poor" "bad" "bad" "good"

99 / 191 100 / 191
Factors (Advanced) Outline
Introduction
Vectors
Calling functions
> which(scores1 == 'bad')
Scripts
[1] 3 5 6 Matrices
Boolean logic
> ## Can do further comparisons with an ordered factor. Lists
> ## Levels are now ordered, as shown by "<" in levels. Factors
> Character arrays
> s2 <- factor(scores1, Objects in your environment
+ levels = c('poor', 'bad', 'good'), Basic plotting
+ ordered = TRUE) Reading/writing data to file system
> s2[1] > s2[2] Writing functions
Conditionals and looping
[1] TRUE Vectorization
Random number generation
Debugging
101 / 191 102 / 191

Strings / character arrays Strings / character arrays

Within a script, easy way to generate output:

> cat("Now computing the steady-state\n")

Character arrays are vectors of strings. Now computing the steady-state

Use single (’) or double (”) quotes to mark strings, but don’t mix:
> x <- 134
> x <- 'good' > cat("sqrt of", x, "is", sqrt(x), "\n")
> z <- "no' # need to match"
> z <- "it's working" sqrt of 134 is 11.57584

> cat("sqrt of", x, "is", sqrt(x), "\n", sep='__')

sqrt of__134__is__11.57584__

See also message, warning and stop for communicating diagnostic

messages – cat and print are generally used when displaying an object.
103 / 191 104 / 191
Strings / character arrays paste0
blackslash characters allow you to generate control characters, importantly:
newline: \n, tab: \t.

> cat("5\t9\n_")

5 9
_ Don’t forget paste0() where sep=”.
(Humour): [Link]
paste() returns a string, e.g. for assignment.
paste0-is-statistical-computings-most-influential-contribution-of-
> x <- 1:5; [Link] <- '/home/stephen/res'
> file <- paste([Link], '/expt_res', x, '.dat', sep='')
> file

[1] "/home/stephen/res/expt_res1.dat"
[2] "/home/stephen/res/expt_res2.dat"
[3] "/home/stephen/res/expt_res3.dat"
[4] "/home/stephen/res/expt_res4.dat"
[5] "/home/stephen/res/expt_res5.dat"
105 / 191 106 / 191

Strings Strings

• Just as R stores vectors of numbers, it also stores vectors of strings.

> grep('e', s)
• Pattern matching facilities are available, based on Unix terms (grep, [1] 1 2 5
regular expressions). These are worth learning:
> grep('^e', s) ## regexps...
> s <- c('apple', 'bee', 'cars', 'danish', 'egg')
[1] 5
> nchar(s)
> sub('e', '_', s)
[1] 5 3 4 6 3
[1] "appl_" "b_e" "cars" "danish" "_gg"
> substr(s, 2,3)
> gsub('e', '_', s) ## global sub, watch "bee"
[1] "pp" "ee" "ar" "an" "gg"
[1] "appl_" "b__" "cars" "danish" "_gg"

107 / 191 108 / 191

Strings Outline
Introduction
Vectors
Calling functions
Scripts
> toupper(s) Matrices
Boolean logic
[1] "APPLE" "BEE" "CARS" "DANISH" "EGG" Lists
Factors
> sprintf('name %s len %d', s, nchar(s)) ## C users!
Character arrays
[1] "name apple len 5" "name bee len 3" Objects in your environment
[3] "name cars len 4" "name danish len 6" Basic plotting
[5] "name egg len 3" Reading/writing data to file system
Writing functions
Conditionals and looping
Vectorization
Random number generation
Debugging
109 / 191 110 / 191

Environments (Advanced) > e <- [Link]()

> e

An environment is a frame, or collection of named objects (variables), and <environment: 0x7ff9e22741f0>

a pointer to an enclosing environment.
> ls(e)
The working environment of your interactive R session is the
<R GlobalEnv>. character(0)

> x <- 2 > [Link](e)

> ls()
<environment: R_GlobalEnv>
[1] "x"
> e$x <- 1
> ## current environment > ls(e)
> environment()
[1] "x"
<environment: R_GlobalEnv>
> e$x != x

[1] TRUE
111 / 191 112 / 191
Inspecting variables and the environment Converting an object from one type to another

Some ways of finding out what kind of object you have and converting.
x <- 9
y<- c(2,4,5) mode(y)
m <- matrix(2:5, 12,10) typeof(y)
objects() ## what vars do I have? class(y)
ls() ## shorthand for objects. [Link](y)
str(m) ## Display compact representation [Link](y) ## [Link]() check if type XYZ
head(m) [Link](y)
tail(m) [Link](m) ## convert object to type XYZ
rm(list = ls()) ## clear up the working environment l <- list(x=c(1,4), y=c(6,9,12))
ls() unlist(l)

113 / 191 114 / 191

What is an object? OO programming in R

• An object is typically either a variable or a function.
• You can use the same name for a function and a variable, and R uses
context to decide which you mean: Data abstraction; object manipluation becomes independent of the
implementation details.
> sum <- 3 + 4 + 5
> total <- sum(1:4) There are several OO programming frameworks in R . The main ones,
> total supported in base R are S3, S4 and S4 ReferenceClasses.

[1] 10 Bioconductor provides many ad hoc classes to store, manipulate and process
microarray, RNA Seq, proteomics, flow cytometry, . . . data.
> sum

[1] 12 It there is interest, we could have a session about OO programming at the

end of the course.
> sum(sum) ## can get confusing!

[1] 12
115 / 191 116 / 191
Outline Basic plotting
Introduction
Vectors
Calling functions
Scripts
Matrices
Boolean logic • Basic x,y plots (scatter plots)
Lists • Multiple plots in one figure
Factors • Saving your plots
Character arrays
Objects in your environment This section will just introduce the mechanics of making basic plots, rather
Basic plotting than worry about interpreting them.
Reading/writing data to file system
Writing functions
Conditionals and looping
Vectorization
Random number generation
Debugging
117 / 191 118 / 191

Basic plotting Basic plotting

> ## Expand on previous plot ...
> x <- seq(from=0, to=2*pi, len=1000) > plot(x,y, main='cos(2x)', type='l', lty=1, bty='n')
> y <- cos(2*x) > y2 <- sin(2*x)
> lines(x,y2, type='l', lty=2)
> plot(x,y) ## just provide data; sensible labelling > same <- which( abs(y - y2) < 0.01)
> points(x[same], y[same], pch=19, col='red', cex=3)
> legend('bottomleft', c("cos(2x)", "sin(2x)"), bty='n', lty=c(1,2))
1.0

●
●
●
●●
●
●
●●
●
●
●●
●
● ●
●
●
●
●
●
●●
●
●
●●
●
●
●●
●
●
●●
●
●
●●
●
●
●● ●●
●
●
●
●
●
●
●●
●
●
●●
● cos(2x)
●●
●
● ●
●● ●
●
●● ●●
●
●
●●
● ●●
● ●
●
● ●
●
●
●
● ●
●
● ●
●
● ●
●
●
●
● ●
●
● ●● ●
●
●● ●● ●
● ●
●
● ● ● ●

1.0
●
● ●
● ●
●
● ●
●
●
● ●
● ●● ●
●
●
● ●
● ●
● ●
●
● ●
● ●● ●
●
●
● ● ● ●
●
0.5

●
● ●● ●
● ●
●
●
●●
●
●
●●
● ●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●●
●
●
●
●
●
●
● ● ●
● ●

0.5
●● ● ●
● ●
● ●● ● ●
●
●
● ● ●
● ●
●
● ● ●
● ●
● ●
● ● ●●
●● ● ●
● ●
0.0

●
● ● ●● ●
● ●
y

●● ● ●
● ●
●
● ●● ●● ●
●
●● ● ●
● ●
0.0
● ● ● ●
● ●
y

●
● ● ● ●●
●
● ●
● ●
● ●
● ● ●
● ●
●● ●
● ● ●
●
●
● ● ●● ●
●
● ● ● ●
−0.5

●
● ●
● ●
● ●
●
● ●
● ●
● ●
●
●
● ●
● ●
● ●
●
●
−0.5

●
● ●
● ●
● ●
●
● ●
● ●
● ●
●
●
● ●
● ●
● ●
●
●
● ●
● ●
● ●
●
●
● ● ●
● ●
●
●
●
●
●●
●
●
●●
●
●
●● ●
●
●
●
●
●
●
●
●
●
●
●
● ●
●●
●
●
●●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
cos(2x)
● ●
−1.0

●
●
●●
● ●
●
●● ●
●
●
● ●●
●
●
●
●
●●
●
● ●
●
●
●● ●
●
●
●● ●
●●
●
●
−1.0

●●
●
●
●●
●
●
●●
●
●
● ●
●
●●
●
●
●●
●
●
●●
● sin(2x)

0 1 2 3 4 5 6 0 1 2 3 4 5 6

x 119 / 191 x 120 / 191

Options controlling the plot Multiple data sources on one plot
When you wish to have multiple data sources on one plot (e.g. two
time-series plots), the approach is to draw the first using plot and then draw
subsequent features using lines or points.
Axes are not rescaled, so draw the bigger plot first.
par() outputs the (long) list of options that control plotting behaviour. > x <- 1:30
> y <- sqrt(x); z <- log(x)
Read ?par for all the details! > plot(x,y); lines(x,z, col='red')
Common options to explore: > plot(x,z, type = "l", col = "red"); points(x,y) ## some data missing
• mfrow, mfcol: multiple plots in figure ●
●
●
●

●
●

• mar, oma: margins around plot and figure. ●

●

3.5
●
●
● ●
●
● ●
●

• ask: whether to hit RETURN between pages of figures.

3.0
●

5
● ●
● ●
●
●
● ●

2.5
●
● ●
●
●

4
● ●
●

2.0
● ●
●
● ●

z
●

1.5
●

3
● ●
●
●

1.0
●
●
●

0.5
●
●

0.0
●

1
0 5 10 15 20 25 30 0 5 10 15 20 25 30

x x

121 / 191 122 / 191

Multiple plots in one figure Saving your plots

mfrow and mfcol are useful parameters within par(), but margins often
R can save plots in many formats, including PDF, postscript, PNG, JPEG.
need to be changed to maximise space.
sin (x) sin (2x) Best to use vector formats (PDF, postscript) for graphs and bitmap formats
1.0

1.0

(png, jpeg) for images.

0.5

0.5
sin(2 * x)

R has output devices, only one of which is active, [Link]().

sin(x)
0.0

0.0

pdf(file = 'mfrow_eg.pdf',
width = 4, height = 6)
−1.0

−1.0

par(mfrow = c(3,2))
par(mar = c(3.5, 3.5, 1.5, 0.5),
mgp = c(2.5, 1, 0))
0 1 2 3
x
4 5 6 0 1 2 3
x
4 5 6 > [Link]()
x <- seq(from = 0, to = 2*pi, sin (3x) cos (x)
> pdf(file='[Link]', width=7, height=7) ## inch
1.0

1.0

len = 100)
> [Link]()
0.5

0.5

plot(x, sin(x), main = "sin (x)",

sin(3 * x)

cos(x)

type = 'l')
0.0

0.0

plot(x, sin(2*x), main = "sin (2x)", > hist( rnorm(9999) )

type = 'l')
> [Link]() ## close device
−1.0
−1.0

plot(x, sin(3*x), main = "sin (3x)",

type = 'l') 0 1 2 3 4 5 6 0 1 2 3 4 5 6
plot(x, cos(x), main = "cos (x)", x
cos (2x)
x
cos (3x)
> png(file='[Link]', w=600, h=600) ## pixels
type = 'l')
1.0

1.0

plot(x, cos(2*x), main = "cos (2x)", > hist( rnorm(9999) )

0.5

0.5

type = 'l')
> [Link]()
cos(2 * x)

cos(3 * x)

plot(x, cos(3*x), main = "cos (3x)",

0.0

0.0

type = 'l')
[Link]()
−1.0

−1.0

0 1 2 3 4 5 6 0 1 2 3 4 5 6
Zoom in on text of PNG to see limitations of this format.
x x

123 / 191 124 / 191

Next steps with plotting (Advanced)

Estimated Density 0 1
60 65 70 75 80

Soprano 2 Soprano 1 6
0.25
0.20
●
0.15 ●
●
0.10
0.05

R has a vast range of functions for plotting particular data types. You may
●●●●●
●●●●● ● ●●●●●●
●●●● 0.00
Alto 2 Alto 1
0.25

read about different packages for plotting: 0.20

0.15
0.10
4
●

factor(cyl)

with Normal Fit

Kernel Density
● ●
●
0.05 ● ● 4

• base graphics (or “traditional”)

wt
● ●
0.00 ●●●●● ●●
● ●●●●●●●●●●● ● ●
● ● ●● ● ● 6
Tenor 2 Tenor 1
●
● ●
● 8
●
0.25

• lattice/grid (lattice is built upon grid)

0.20 ●
●
0.15 ●
●
0.10
●
0.05

• ggplot2 – [Link]
●
● ●●●● ● ● ●●●●●●●●●●● ● 0.00 ●
●
Bass 2 Bass 1 2 ●
●
0.25
0.20 ●
●
0.15
0.10
0.05
0.00 ●●●●●●● ●● ● ●●●●●●●●
●

60 65 70 75 80
10 15 20 25 20 25 30
Height (inches) mpg

125 / 191 126 / 191

Outline
Introduction
Vectors
Calling functions
Here are some starting points to explore: Scripts
• demo(graphics) to see diversity of plots. Matrices
• low-level functions: symbols(), rect(), segments(), abline(). Boolean logic
Lists
• curve(x*log(x))
Factors
• R graphics gallery
Character arrays
[Link] Objects in your environment
[Link]
Basic plotting
Reading/writing data to file system
Writing functions
Conditionals and looping
Vectorization
Random number generation
Debugging
127 / 191 128 / 191
Reading/writing data to file system Interacting with the file system

• where am I currently? getwd()

• What’s my current directory? dir, getwd, setwd • change me to a new directory: setwd("/tmp")
• scan, readLines (GUIs have chooser for interactively changing directory.)
• [Link], [Link], [Link] • What files are in my [current] directory?
• RData files – save and load
> dir()
• Further I/O functions > dir("/tmp")
> dir(pattern="\\.R$") ## regexps, see later.

129 / 191 130 / 191

Scan, write, readLines Scan, write, readLines

For basic reading/writing of data, use scan/write. Filenames are specified
relative to current directory. Can even give URL as a file. Files often have a
header which can be skipped over.
> x <- scan('files/[Link]', skip=1)
> [Link] <- round( runif(100, min=5, max=10), 2)
> summary(x) > tf <- tempfile()
> tf
Min. 1st Qu. Median Mean 3rd Qu. Max.
14.00 19.25 27.25 31.58 39.25 65.00 [1] "/var/folders/rw/2g0_whns55x580hlgf1vjs8m0000gn/T//RtmpwUmOlm/file128
> s1 <- readLines('files/[Link]') ## treats as strings
> summary(s1)
> write([Link], tf)
> s <- scan(tf)
Length Class Mode > [Link](s, [Link])
5 character character
[1] TRUE
> h <- scan('[Link] ## [Link]
> summary(h)

Min. 1st Qu. Median Mean 3rd Qu. Max.

136.2 158.3 164.9 163.0 169.2 178.7

131 / 191 132 / 191

[Link] et al. Rdata files
Text files are useful for portably storing data, so that they can be read
If data are tabular, [Link] or [Link] is often useful. (Useful for
across applications. R has its own format for efficiently storing objects. Files
importing spreadsheets; just save as a comma separated value file, CSV.)
much smaller than text files. However, this format is not universally known.

x <- [Link]('./files/[Link]', sep = '\t', > n <- 99999; x <- rnorm(n)

header = TRUE) > [Link] <- tempfile()
names(x) > [Link] <- tempfile()
head(x) > write(x, n, file = [Link])
x[2,] > save(x, n, file = [Link])
x$Goals >
[Link](x) > ## Compare sizes of files with the object.
tf <- tempfile() > [Link](x)
[Link](x, tf, [Link]=FALSE) > [Link]([Link])$size
## sort by goals scored. > [Link]([Link])$size ## compressed
x[order(x$Goals, decreasing = TRUE), ] >
> rm(x,n)
See ?[Link]. > load([Link]) ## reload data.
133 / 191 134 / 191

Saving your workspace with .RData files Further I/O functions (Advanced)

When you quit R , you are asked:

> q()
Save workspace image? [y/n/c]:
If you answer y, all objects (variables and functions) in your global
R has many facilities for I/O. See for example the following help topics.
environment are saved for future use, using [Link]. From ?save:
• ?connections — interface to files, pipes, sockets, compressed files . . .
'[Link]()' is just a short-cut for "save my current
• ?sink — divert R output to a connectin
workspace", i.e., 'save(list = ls(all=TRUE), file = ".RData")'.
It is also what happens with 'q("yes")'. • ?dget / ?dput — read/write ASCII representation of an R object.

Warning: If an .RData file is present in the current directory when starting

R , it is silently loaded. I think it can be dangerous, as you may not realise Also: XML, DBMS, SQLite, netCDF, hdf5, . . .
what values have been silently loaded. Better to instead be explicit:
> [Link](file='[Link]') ## keep everything
> save(x, y, z, file='[Link]') ## or just key objects
> ## ...
> load('[Link]') ## reload objects

or start R using R --no-restore

135 / 191 136 / 191
Outline Writing functions: overview
Introduction
Vectors
Calling functions
Scripts
Matrices
Boolean logic • Why bother?
Lists • How to write (local args, return value; cannot change value)
Factors • Example: computing std. deviation
Character arrays • Local variables within functions
Objects in your environment
• Recursion.
Basic plotting
Reading/writing data to file system
Writing functions
Conditionals and looping
Vectorization
Random number generation
Debugging
137 / 191 138 / 191

Functions Functions

• How to define a new function:

• Functions promote code reuse. > [Link] <- function(arg1, arg2) {
• Black-box approach; given inputs, what output should I expect? This + ## Doc string here.
requires good documentation of what your function does. Can it be + x <- arg1 * 2
described without having to look at the code? + y <- sqrt(arg2) + 5
• Finding the right level of definition for a function is hard, and how to + z <- x * y
modularise comes with experience. Typically rewrite many times before + ## last value is the return value of the function.
getting final solution + ## Use a list to return several items.
+ z ## same as return(x)
+ }

139 / 191 140 / 191

Example of writing a new function Terminology of variables within functions
Compute the standard deviation of a vector of numbers:
sP Pn
n 2 • In [Link], x is the name of a formal argument. In the following, y is
i=1 (xi − x̄) xi
[Link] = where x̄ = i=1 called the actual argument (doesn’t have to be named x – can be
n−1 n
named however you wish).

> [Link] <- function(x) { > n <- 5

+ ## Return std dev of X. > y <- c(9, 2, 7, 10)
+ n <- length(x) > [Link](y)
+ xbar <- sum(x)/n [1] 3.559026
+ diff <- x - xbar
+ [Link] <- sum( diff^2) > print(n) ## should still be 5, not 4.
+ var <- [Link] / (n-1) [1] 5
+ ## last value calculated is return value.
• Local variables within function are not available outside of function.
+ sqrt(var)
+ }

141 / 191 142 / 191

Handling unbound variables

Variables created by assignment within a function are known as local
• Any change to formal args within a function does not change value of variables (e.g. y below). If a variable is not local, or a formal argument, it is
actual argument outside the function: an unbound variable. It may be found in the enclosing environment
(typically the global workspace), or an error is generated – this is bad
> [Link] <- function(x) {
practice!
+ x <- x^2 ## change internally
+ sum(x) > fn1 <- function(x) {
+ } + y <- x^2
> y <- c(4, 5, 6) + res <- sum( (y - thresh)^2 )
> [Link](y) + res
+ }
[1] 77 > dat <- 1:5
> fn1(dat) ## case 1
> y
[1] 4 5 6 Error in fn1(dat): object ’thresh’ not found

> thresh <- 10

> fn1(dat) ## case 2

[1] 379
143 / 191 144 / 191
Handling unbound variables (2) The ... arguments
In this case, better to define thresh as an argument of the function, and
provide a default value:

> fn1 <- function(x, thresh=10) { From R language definition

+ y <- x^2
The ... argument is special and can contain any number of
+ res <- sum( (y - thresh)^2 )
arguments. It is generally used if the number of arguments is
+ res
unknown or in cases where the arguments will be passed on to
+ }
another function.
> fn1(dat) ## case 3
See ?cat for an example of number of arguments is unknown.
[1] 379

Advanced: use codetools::checkUsage() to find unbound vars;

codetools::findGlobals() for globals.

145 / 191 146 / 191

The ... arguments Writing a replacement function (Advanced)

> p <- function(x, y, ...) { Convention for a replacement function is that the name should end with <-.
+ if (diff(range(x)) > diff(range(y))) {
The last argument of the replacement function must be called value and is
+ plot(x, y, ...)
+ } else { the RHS of the assignment.
+ plot(y, x, ...)
> "threshold<-" <- function(x, value) {
+ } + ## X is the object to update
+ } + ## VALUE is the value on the RHS.
> p(rnorm(10, 0, 1), rnorm(10, 0, 5), main = "Expect swap", + y <- ifelse(x > value, 1, 0)
+ return(y)
+ col = "red") + }
> p(rnorm(10, 0, 5), rnorm(10, 0, 1), main = "No swap", pch = 19) > x <- c(0.3, 0.1, 0.6, 0.7, 0.9, 0.2)
> threshold(x) <- 0.4
Expect swap No swap
> x
2.5

● ●
0.5

● ● ●
[1] 0 0 1 1 1 0
1.5

●
●
0.0
x

● ●
0.5

In general, replacement functions are used to set variable attributes or

−0.5

● ●
●
● ●
−0.5

object slots.
●
●
−1.0

●
● ●

−10 −5 0 5 −5 0 5 10

y x

147 / 191 148 / 191

Tips for writing functions Outline
Introduction
Vectors
• Can you think of a way to break down the problem so that a team can
Calling functions
work on the problem, with each person assigned to a independent
Scripts
piece? “Divide + conquer”.
Matrices
• Each function should be easy to test, then you can “freeze” it. Write
Boolean logic
test cases, which can be automatically checked.
Lists
> [Link]([Link](100,200), 300) Factors
Character arrays
• Rule of thumb: each function should be no more than a page or two of
Objects in your environment
code. Basic plotting
• For large projects, avoid mixing computation and plotting in the same Reading/writing data to file system
function – separate the two jobs; this makes it easier to run in batch. Writing functions
> res <- [Link](par1, par2, par3) Conditionals and looping
> [Link](res) Vectorization
Random number generation
Debugging
149 / 191 150 / 191

Control-flow constructs if / if ... else ...

> x <- 8;
>
> if (x > 10) {
• if + ## condition was true
• switch + cat("x is bigger than 10\n")
+ } else {
• for
+ cat("x is 10 or less\n")
• while
+ }
• Vectorization
• simple applications – numerics x is 10 or less

Notes:
“else ...” can be omitted if you do not need it.
if returns a value, which can be assigned, e.g. y <- if (x <10) 40
else 20. A better solution in this case however is the vectorized form
y <- ifelse(x<10, 40, 20)
151 / 191 152 / 191
Braces in conditional constructs switch (Advanced)

Curly braces not needed if there is only one expression in the if clause: Nested if ... else commands can get a bit messy. Like other languages, R has
a switch construct. From ?switch:
if ( x > 10 ) { if ( x > 10 )
y <- 1 y <- 1 > centre <- function(x, type) {
} ## + switch(type,
+ mean = mean(x),
But braces are needed in multiline if/else statement: + median = median(x),
+ trimmed = mean(x, trim = .1))
+ }
if ( x > 10 ) { if ( x > 10 ) > x <- rcauchy(10)
y <- 1 y <- 1 > centre(x, "mean")
} else { else
y <- 0 ## OK y <- 0 ## NOT OK [1] -2.801097
} ##
> centre(x, "median")
From ?Control: Note that it is a common mistake to forget to put
braces ('{ .. }') around your statements, e.g., after 'if(..)' or [1] 0.08003974
'for(....)'. In particular, you should not have a newline between
> centre(x, "trimmed")
'}' and 'else' to avoid a syntax error in entering a 'if ... else'
construct at the keyboard or via 'source'. For that reason, one
[1] -0.05755348
(somewhat extreme) attitude of defensive programming is to always
use braces, e.g., for 'if' clauses. 153 / 191 154 / 191

Recursive functions Looping constructs

Here is an example of using conditionals with a divide and conquer
approach; quicksort in a few lines (albeit not very efficient).
> qsort <- function(data) {
+ ## Sort data into ascending order. Looping constructs allow you to repeat calculations as many times as you
+ n <- length(data)
+ if (n <= 1) {
wish. This is why computers are so useful – it is just as easy (usually) to
+ data repeat something 1000 times as 10 times.
+ } else {
+ pivot <- data[floor(n/2)]
+ less <- data[which(data < pivot)]
+ equal <- data[which(data == pivot)] e.g. if you want to simulate flipping a (biased) coin 100 times, and counting
+ greater <- data[which(data > pivot)]
+ c( qsort(less), equal, qsort(greater)) the number of heads, no problem. If you want to repeat this process 1000
+ } times, no problem. (See later.)
+ }
> all(replicate(99, {
+ data <- runif(2000, max=10)
+ [Link](qsort(data), sort(data)) }))

[1] TRUE

155 / 191 156 / 191

for loops while loops
for (var in seq) command while (condition) {
command
seq is a vector; var is set in turn to each value in the vector, and then command
command executed. Multiple commands can be given within braces. }
e.g. So the commands are executed until the condition is no longer true.
Typically then one of the commands will change the condition.
> x <- 6 e.g. print all the Fibonacci numbers (f[i] = f[i-1] + f[i-2]) less than
> for (i in 1:3) { 100.
+ res <- x * i
+ cat(x, "*", i, "=", res, "\n") > n1 <- 0; n2 <- 1
+ } > while (n2 < 100) {
+ print(n2)
6 * 1 = 6 + old <- n2
6 * 2 = 12 + n2 <- n2 + n1
6 * 3 = 18 + n1 <- old
+ }
157 / 191 158 / 191

Breaking out of loops Breaking out of loops

repeat expr will repeatedly execute expr until you break out of the loop. next allows you to skip to next iteration of a loop. Both next and break
can be used within other loops (while, for).
> i <- 3
> repeat { > for (i in 1:10) {
+ if (i == 5) { + if ((i %% 2) == 0)
+ break + next
+ } else { + print(i)
+ cat("i is", i, "\n") + }
+ i<- i+1
+ } [1] 1
+ } [1] 3
[1] 5
i is 3 [1] 7
i is 4 [1] 9

159 / 191 160 / 191

A word on indentation Outline
Introduction
Indentation helps you see the flow of the logic, rather than flattened version. Vectors
(Use tab key to indent). Reformatting tools are available (e.g. within Calling functions
Emacs). Scripts
Matrices
## version 1. ## version 2. Boolean logic
i <- 3 i <- 3 Lists
repeat { repeat { Factors
if (i==10) { if (i==10) { Character arrays
break break Objects in your environment
} else { } else { Basic plotting
cat("i is", i, "\n") cat("i is", i, "\n") Reading/writing data to file system
i<- i+1 i<- i+1 Writing functions
} } Conditionals and looping
} } Vectorization
Indentation helps to show structure, and match braces. Random number generation
Debugging
161 / 191 162 / 191

Vectorization Vectorization example

When possible, operate on vectors, rather than using for loops. Q: Flip a biased coin [p=0.6 of heads] 100 times; how many heads do you
Rewrite code, but beware sometimes not possible (Fibonacci). get? Repeat this for 1000 trials.
Compute difference between times of events, e. Given n events, there will
> n <- 100 ## number of coin flips in trial
be n-1 inter-event times. interval[i] <- e[i+1] - e[i]
> diff1 <- function(e) {
> p <- 0.6 ## prob of getting heads
+
+
n <- length(e)
interval <- rep(0, n-1) ## good to pre-alloc!
> ntrials <- 1000
+ for (i in 1:(n-1)) {
+ interval[i] <- e[i+1] - e[i]
+ }
+ interval > trial1 <- function(n, [Link]) {
+ }
> diff2 <- function(e) { + count <- 0
+ n <- length(e) + for (i in 1:n) {
+ e[-1] - e[-n] + if (runif(1) < [Link])
+ }
> e <- c(2, 5, 10.2, 12, 19) + count <- count +1
> diff2(e) + }
+ count
[1] 3.0 5.2 1.8 7.0 + }
>
> [Link](diff1(e), diff2(e)) > res <- rep(0, ntrials)
> for (j in 1:ntrials) {
[1] TRUE + res[j] <- trial1(n, p)
+ }
Advantages: shorter, more readable and faster (no loops). > hist(res)
163 / 191 164 / 191
Vectorization example Vectorization example
Q: Flip a biased coin [p=0.6 of heads] 100 times; how many heads do you Q: Flip a biased coin [p=0.6 of heads] 100 times; how many heads do you
get? Repeat this for 1000 trials. get? Repeat this for 1000 trials.
> n <- 100 ## number of coin flips in trial > n <- 100 ## number of coin flips in trial
> p <- 0.6 ## prob of getting heads > p <- 0.6 ## prob of getting heads
> ntrials <- 1000 > ntrials <- 1000

> trial1 <- function(n, [Link]) { > trial2 <- function(n, [Link]) { > trial1 <- function(n, [Link]) { > trial2 <- function(n, [Link]) {
+ count <- 0 + [Link] <- runif(n) + count <- 0 + [Link] <- runif(n)
+ for (i in 1:n) { + sum( [Link] < [Link]) + for (i in 1:n) { + sum( [Link] < [Link])
+ if (runif(1) < [Link]) + } + if (runif(1) < [Link]) + }
+ count <- count +1 > + count <- count +1 >
+ } > res <- replicate(ntrials, + } > res <- replicate(ntrials,
+ count + trial2(n, p)) + count + trial2(n, p))
+ } > hist(res) + } > hist(res)
> >
> res <- rep(0, ntrials) > res <- rep(0, ntrials)
> for (j in 1:ntrials) { > for (j in 1:ntrials) {
+ res[j] <- trial1(n, p) + res[j] <- trial1(n, p)
+ } + }
> hist(res) > hist(res)
165 / 191 In this case, hist( rbinom(1000, 100, 0.6)) would also work! 166 / 191

apply family apply family

e.g. how to compute sum of each row of a matrix? sum(A) will normally Other functions: lapply (apply to list), sapply (simplify), replicate, . . .
return the sum of all elements of A.
> lapply(ls(), [Link])
apply(X, MARGIN, FUN, ...) > sapply(ls(), [Link])
MARGIN = 1 for row, 2 for cols. > tapply(1:20, factor(rep(letters[1:5], each = 4)), sum)
FUN = function to apply > ## see also ?by
... = extra args to function. > mapply(sum, 1:5, 1:5, 1:5)
> hist( replicate(200, mean(rnorm(100)) ))

> A <- matrix(1:6, 2,3) How to *apply – [Link]

> [Link] <- apply(A, 1, mean) ## or rowMeans
> [Link] <- apply(A, 2, sum, [Link]=TRUE) ## colSums Exercise: Why is this better than writing a for loop?
(Eglen, 2009); parallel package: mclapply(), parLapply(), . . . .

167 / 191 168 / 191

Anonymous functions (Advanced) Fibonacci sequence

Sometimes you don’t want to pollute name space by defining a new function,
so just use an “anonymous function”, i.e. a function without a name. 0, 1, 1, 2, 3, 5, 8, 13, 21, . . .
Particularly useful e.g. in an apply call.
f [n] = f [n − 1] + f [n − 2]
> [Link] <- matrix(1:10, ncol=5)
How to vectorize?
> apply([Link], 2, function(x) { sum(x^2) + 10 })
Exercise: write a function, fibonnaci(n) that returns the nth element of
the sequence. Assume that fibonnaci(1) == 0, fibonacci(2) == 1.
Since functions are just objects, anonymous functions are just objects Exercise: use fibonacci() to estimate the golden ratio.
without names, similar to ’anonymous numbers’ like a+b in an expression
a+b+c.

169 / 191 170 / 191

Efficiency A final word on efficiency

“premature optimization is the root of all evil” (Knuth, 1974).

Examples adopted from [Link]/res/code_segments

f1 is bad; should pre-allocate vector, rather than rely on R to allocate
memory repeatedly3 .
• Rule 1 of optimization: don’t bother (Kernighan).
> f1 <- function() { > f2 <- function() { • For loops are not always a bad thing. See last example in Help Desk
+ n <- 1e4; decay <- 0.9995 + n <- 1e4; decay <- 0.99995 article (May 2008).
+ out <- rep(0, n) ##pre-alloc
+ out <- 1.0 + out[1] <- 1.0 [Link]
+ for (i in 2:n) + for (i in 2:n)
+ out[i] <- out[i-1] * decay + out[i] <- out[i-1] * decay
+ out + out
+ } + }
> [Link](o1 <- f1()) > [Link](o2 <- f2())

user system elapsed user system elapsed

0.396 0.128 0.525 0.036 0.001 0.037

3
This is not optimisation, but rather good coding practice.
171 / 191 172 / 191
Numerics issues Numerics issues
Although integer arithmic is reliable, floating-point arithmetic is to be
treated with care! (All R ’s calculations are in what C programmers call
“double precision”.) Compare:
> 1.0 + 2.0 == 3.0

[1] TRUE Solution - testing near equality

> 0.1 + 0.2 == 0.3 > [Link]( 0.1 + 0.2, 0.3)

[1] FALSE [1] TRUE

From FAQ (7.31?)

> a <- sqrt(2)
> a * a == 2

[1] FALSE

> a * a - 2

[1] 4.440892e-16

173 / 191 174 / 191

How big is infinity? How small is epsilon?

How big can  be such that 1 +  = 1? (Taken from Goldberg (1991) ACM
Use while loop to estimate it: article, p220).
> x <- 1 > eps <- 1
> while ( [Link](x*2) ) { > while (eps + 1 > 1) {
+ x <- x*2 + eps <- eps * 0.5
+ } + }
> > eps ##1.110223e-16
> x ## 8.988466e+307
[1] 1.110223e-16
[1] 8.988466e+307
> 1 + eps ## 1
> x*2 ## Inf
[1] 1
[1] Inf
> (1 + eps == 1) ## TRUE
> (x*2)/2 ## Inf
[1] TRUE
[1] Inf
> 1 + (2*eps) ##1
> .Machine$[Link]
[1] 1
[1] 1.797693e+308
> (1 + (2*eps) == 1) ## FALSE

[1] FALSE
175 / 191 176 / 191
Outline Random number generation
Introduction
Vectors
Computers usaully generate “pseudo-random numbers”. They are generated
Calling functions
based on some iterative formula:
Scripts
Matrices xnew = f (xold ) mod N
Boolean logic
Lists where modulo operation provides the “remainder” division.
Factors To generate the first random number, you need a seed.
Character arrays Setting the seed allows you to reliably generate the same sequence of
Objects in your environment numbers, which can be useful when debugging programs.
R has many routines for generating random samples from various
Basic plotting
distributions, but for now we will just use runif(), (and maybe rnorm()).
Reading/writing data to file system
Exercise: write a random number generator. See: “Randu: a bad random
Writing functions
number generator”. [Link]
Conditionals and looping
Exercise: Apply the central limit theorem to generate samples from a
Vectorization
normal distribution by adding together samples from a uniform distribution.
Random number generation
Debugging
177 / 191 178 / 191

Outline Debugging (Advanced)

Introduction
Vectors
Calling functions
Scripts See An introduction to the Interactive Debugging Tools in R , Roger
Matrices D Peng for detailed usage.
Boolean logic [Link]
Lists • warnings vs errors; converting warnings to errors; stopifnot().
Factors • what to do when I get an error: traceback()
Character arrays
• simple print statements are often useful.
Objects in your environment
• Use of browser() at key points in code.
Basic plotting
Reading/writing data to file system • debug(fn), undebug(fn)
Writing functions • Using recover() rather than browser()
Conditionals and looping
Vectorization
Random number generation
Debugging
179 / 191 180 / 191
Warnings and errors Traceback
• A warning is softer than an error; if a warning is generated your
program will still continue, whereas an error will stop the program.
> log(c(2, 1, 0, -1, 2)) When your program generates an error, use traceback() to find out where
it went wrong:
Warning in log(c(2, 1, 0, -1, 2)): NaNs produced
[1] 0.6931472 0.0000000 -Inf NaN 0.6931472 > start <- function() { go( sqrt(10)) }
> go <- function(x) { inner(x, '-13')}
> xor( c(TRUE, FALSE))
> inner <- function(a, b) {
Error in xor(c(TRUE, FALSE)): argument "y" is missing, + c <- sqrt(b)
with no default + a * log(c)
• If you try to isolate warnings, you can change warnings to errors: + }
options(warn=2). See ?options for further details. > start() ## error
• Add warnings and errors to your code using warning(), stop(). > traceback() ## postmortem debugging
• Can add “assertions” into your code to check that certain values hold.
> stopifnot(x>0)
• Other useful safety checks: all(x>0), any(x>0)
181 / 191 182 / 191

Single-stepping through your code Safety-checks: browser

Here’s a possible usage of browser() that I have in my code:
f i n d . h i g h <− f u n c t i o n ( x , t ) {
Use browser() to single-step through your code. Place it within your ## R e t u r n s a m p l e s i n x b i g g e r t h a n t .
## ( B e t t e r t o u s e x [ x>t ] i n r e a l − l i f e ! )
function at the point you want to examine (e.g.) local variables. max . l e n g t h <− 100 ## s h o u l d be u p p e r l i m i t . . .
r e s u l t s <− r e p ( 0 , max . l e n g t h )
c o u n t e r <− 0
Can use debug([Link]) to step through entire function. for ( i in x) {
undebug() will remove that debug call. if ( i > t) {
Within the browser, you can enter expressions as normal, or you can give a c o u n t e r <− c o u n t e r + 1
i f ( c o u n t e r > max . l e n g t h ) {
few debug commands: browser ()
} else {
• n: single-step r e s u l t s [ c o u n t e r ] <− i
}
• c: exit browser and continue }
• Q: exit browser and abort, return to top-level. }
r e s u l t s [ 1 : counter ]
• where: show stack trace. }
x <− rnorm ( 1 0 0 )
Debug on stddev.R f i n d . high (x , 0.7)

x <− rnorm ( 1 0 0 0 )
(1− pnorm ( 0 . 7 ) ) ∗ l e n g t h ( x ) ## e x p e c t e d .
f i n d . high (x , 0.7)
183 / 191 184 / 191
recover Packages

recover() is like browser(), except you can choose which level to inspect, • R has a packaging system for external code.
rather than the level at which browser was called. • A package is loaded from a library using library("[Link]").
Following allows recover() to be launched when you hit an error:
• Beware: don’t call a package a library! A library is a group of folders
o p t i o n s ( e r r o r=r e c o v e r ) where packages are stored . . .
Here we simply tell R that when an error is generated, we call the function
“recover”. The default is NULL, in which case stop is called. > library() ## view available packages
> library(help="cluster") ## what's in this?
From ?options: > library("cluster") ## load package
Note that these need to specified as e.g. > example(pam) ## can use pam and friends.
options(error=utils::recover) in startup files such as .Rprofile. > detach("package:cluster") ## remove pkg.

185 / 191 186 / 191

CRAN: Comprehensive R Archive Network Managing libraries

CRAN: Site(s) for downloading R , and also its many contributed
*packages*.
If you do not have write access to the default library, R will create a local
Mac/Win have a GUI for installing packages, or it can be done on the
one in your home directory. Multiple libraries are supported:
command line:
> .libPaths()
[Link](c("splancs", "sp"))
R CMD INSTALL [Link] ## from shell [1] "/Users/stephen/NOBACKUP/RLIB"
[2] "/Library/Frameworks/[Link]/Versions/3.2/Resources/librar
If asked to selected a CRAN mirror, in UK use:
[Link] You can also set R libraries as a global shell environment in your .bashrc
file (e.g. PWF linux)
For Bioconductor
export R_LIBS=$HOME/NOBACKUP/RLIB
> source("[Link]
(Be careful! Check that you are not overwriting an existing R LIBS setting.)
> library("BiocInstaller")
> biocLite("Biobase")

187 / 191 188 / 191

Bioconductor Other topics of interest (Advanced)

A success story of R . Started 2001 with aims to:

• S3, S4 and S4 Reference classes for OOP.
• provide access to stat/graphical methods for analysis of genomic data.
• Building your own packages. Useful for packaging up your code, data
• link seamlessly to on-line databases (PubMed/GenBank). sets and documentation. You may wish to do this for large projects
• allow rapid development of extensible software. that you wish to share with others. Read Writing R Extensions manual
• provide training in methods (short courses). and see [Link] to get started.
• promote software with high quality docs and reproducible research • Access to databases. Computational Biology datasets are often quite
(vignettes) . . . large, and you might wish to access data via databases. R package DBI
• Gentleman et al. (2004) Genome Biology 5:R80. provides common interface to SQLite, MySQL, Oracle. See Gentleman
[Link] (2008), Chapter 8.

189 / 191 190 / 191

Reproducible research: Sweave, knitr and vignettes

• Use one file to store code and document. Best shown by way of
example... [Link]

• Vignettes often used in Bioconductor to document packages.

• Interactively explore vignettes.

> library("tkWidgets")
> vExplorer()
• knitr is the natural successor to Sweave. It rocks!

191 / 191