Permission is granted to make and distribute verbatim copies of this manual provided the copyright notice and this permission notice are preserved on all copies. Permission is granted to copy and distribute modified versions of this manual under the conditions for verbatim copying, provided that the entire resulting derived work is distributed under the terms of a permission notice identical to this one. Permission is granted to copy and distribute translations of this manual into another language, under the above conditions for modified versions, except that this permission notice may be stated in a translation approved by the R Development Core Team. ISBN 3-900051-12-7

Suggestions to the reader
Most R novices will start with the introductory session in Appendix A. We have made a number of small changes to reflect differences between the R and S programs. page 67 on the graphics facilities can be read at almost any time and need not wait until all the preceding sections have been digested.org. Smith when at the University of Adelaide. Please address email correspondence to R-core@R-project.Preface
1
Preface
This introduction to R is derived from an original set of notes describing the S and SPlus environments written in 1990–2 by Bill Venables and David M. Comments and corrections are always welcome. Many users will come to R mainly for its graphical facilities. Chapter 12 [Graphics]. We would like to extend warm thanks to Bill Venables (and David Smith) for granting permission to distribute this modified version of the notes in this way. This should give some familiarity with the style of R sessions and more importantly some instant feedback on what actually happens.
. and for being a supporter of R from way back. and expanded some of the material. In this case.

and documentation for S/S-Plus can typically be used with R. integrated collection of intermediate tools for data analysis. Chambers and Trevor J. The new features of the 1991 release of S are covered in Statistical Models in S edited by John M. However. as is frequently the case with other data analysis software. • a suite of operators for calculations on arrays. It has developed rapidly. The formal methods and classes of the methods package are based on those described in Programming with Data by John M. John M. and • a well developed.
1. most programs written in R are essentially ephemeral. in particular matrices. • a large. Chambers. and also forms the basis of the S-Plus systems. user defined recursive functions and input and output facilities.1 The R environment
R is an integrated suite of software facilities for data manipulation. loops. and has been extended by a large collection of packages. Chambers and Allan R. coherent. rather than an incremental accretion of very specific and inflexible tools. Becker.Chapter 1: Introduction and preliminaries
2
1 Introduction and preliminaries
1. For R. • graphical facilities for data analysis and display either directly at the computer or on hardcopy. Among other things it has • an effective data handling and storage facility. There are now a number of books which describe how to use R for data analysis and statistics. We prefer to think of it of an environment within which many classical and modern statistical techniques have been implemented. page 102. written for a single piece of data analysis.
1. calculation and graphical display. keeping the differences between the S implementations in mind. A few of these are built into the base R environment. the basic reference is The New S Language: A Programming Environment for Data Analysis and Graphics by Richard A. but many are supplied as packages.2 Related software and documentation
R can be regarded as an implementation of the S language which was developed at Bell Laboratories by Rick Becker. There are about 25
. The evolution of the S language is characterized by four books by John Chambers and coauthors. Hastie. See Appendix F [References]. See Section “What documentation exists for R?” in The R statistical system FAQ. Wilks. yet many people use R as a statistics system. for precise references. John Chambers and Allan Wilks. R is very much a vehicle for newly developing methods of interactive data analysis.3 R and statistics
Our introduction to the R environment did not mention statistics.) The term “environment” is intended to characterize it as a fully planned and coherent system. simple and effective programming language (called ‘S’) which includes conditionals. (Indeed most of the system supplied functions are themselves written in the S language.

Thus whereas SAS and SPSS will give copious output from a regression or discriminant analysis. More details on packages are given later (see Chapter 13 [Packages].
1. Setting up a workstation to take full advantage of the customizable features of R is a straightforward if somewhat tedious procedure. it is easy to change to a different R prompt if you wish. To quit the R program the command is > q() At this point you will be asked whether you want to save the data from your R session. but users may need to be prepared to do a little work to find it. and on others you will receive a text prompt to which you can respond yes. In S a statistical analysis is normally done as a series of steps. In particular we will occasionally refer to the use of R on an X window system although the vast bulk of what is said applies generally to any implementation of the R environment. If you are running R under Windows or Mac OS you will need to make some small adjustments. say ‘work’. Users in difficulty should seek local expert help.5 Using R interactively
When you use the R program it issues a prompt when it expects input commands. page 82). 4. This guide is aimed at users who have this facility. as we shall see. We will assume that the UNIX shell prompt is ‘$’.4 R and the window system
The most convenient way to use R is at a graphics workstation running a windowing system. Most classical statistics and much of the latest methodology is available for use with R. R will give minimal output and store the results in a fit object for subsequent interrogation by further R functions. Start the R program with the command $ R 3. Most users will find it necessary to interact directly with the operating system on their computer from time to time. we mainly discuss interaction with the operating system on UNIX machines. The default prompt is ‘>’. and so it may appear that nothing is happening. with intermediate results being stored in objects. In using R under UNIX the suggested procedure for the first occasion is as follows: 1.R-project.
1.Chapter 1: Introduction and preliminaries
3
packages supplied with R (called “standard” and “recommended” packages) and many more are available through the CRAN family of Internet sites (via http://CRAN.org) and elsewhere. and will not be considered further here. This will be the working directory whenever you use R for this particular problem. Create a separate sub-directory. which on UNIX might be the same as the shell prompt. There is an important difference in philosophy between S (and hence R) and the other main statistical systems. $ mkdir work $ cd work 2. However. On some systems this will bring up a dialog box. no or cancel (a single letter abbreviation will
. At this point R commands may be issued (see later). to hold data files on which you will use R for this problem. In this guide.

Then launch R by double clicking on the icon.start() is particularly useful as it is contains a high-level concept list which searches though available functions. terminating with the q() command at the end of the session. 1.start() which will launch a Web browser that allows the help pages to be browsed with hyperlinks.
1. Data which is saved will be available in future R sessions. making it a “character string”: This is also necessary for a few words with syntactic meaning including if.6 An introductory session
Readers wishing to get a feel for R at a computer before proceeding are strongly advised to work through the introductory session given in Appendix A [A sample session]. It can be a great way to get your bearings quickly and to understand the breadth of what R has to offer. quit without saving. Our convention is to use double quote marks for preference. for and function. For example. or return to the R session. Make ‘work’ the working directory and start the program as before: $ cd work $ R 2. To get more information on any specific named function.Chapter 1: Introduction and preliminaries
4
do) to save the data before quitting. and set that in the ‘Start In’ field in your R shortcut.7 Getting help with functions and features
R has an inbuilt help facility similar to the man facility of UNIX. Further R sessions are simple. To use R under Windows the procedure to follow is basically the same. > ??solve Try ?help.search command (alternatively ??) allows searching for help in various ways.
1. page 84. The help. The ‘Search Engine and Keywords’ link in the page loaded by help. as in the string "It’s important". > help("[[") Either form of quote mark may be used to escape the other. Create a folder as the working directory. the argument must be enclosed in double or single quotes. for example solve. Use the R program. subsequent help requests are sent to the HTML-based help system. The examples on a help topic can normally be run by
. the command is > help(solve) An alternative is > ?solve For a feature specified by special characters. On most R installations help is available in HTML format by running > help. On UNIX.search for details and more examples.

Comments can be put almost2 anywhere. The set of symbols which can be used in R names depends on the operating system and country within which R is being run (technically on the locale in use). It is case sensitive as are most UNIX based packages. it is evaluated. An assignment also evaluates an expression and passes the value to a variable but the result is not automatically printed. and amongst those which do some will silently discard the excess and some will use it as the start of the next line. If a command is not complete at the end of a line. not inside strings. so A and a are different symbols and would refer to different variables.Chapter 1: Introduction and preliminaries
5
> example(topic ) Windows versions of R have other optional help systems: use > ?help for further details. Once a command is located in this way. page 95.
1 2 3
For portable R code (including that to be used in R packages) only A–Za–z0–9 should be used.
. Elementary commands consist of either expressions or assignments. If an expression is given as a command. or by a newline.’ the second character must not be a digit. and if it starts with ‘. starting with a hashmark (‘#’).’).’ and ‘_’.8 R commands. by default + on second and subsequent lines and continue to read input until the command is syntactically complete. case sensitivity.9 Recall and correction of previous commands
Under many versions of UNIX and on Windows. Elementary commands can be grouped together into one compound expression by braces (‘{’ and ‘}’). Commands are separated either by a semi-colon (‘. and characters can be removed with the DEL key or added with the other keys. Command lines entered at the console are limited3 to about 4095 bytes (not characters). and the value is lost. R will give a different prompt.
1. The vertical arrow keys on the keyboard can be used to scroll forward and backward through a command history. The recall and editing capabilities under UNIX are highly customizable.
Technically R is an expression language with a very simple syntax. This prompt may be changed by the user. We will generally omit the continuation prompt and indicate continuation by simple indenting. everything to the end of the line is a comment. R provides a mechanism for recalling and re-executing previous commands. More details are provided later: see Appendix C [The command-line editor]. the cursor can be moved within the command using the horizontal arrow keys. Normally all alphanumeric symbols are allowed1 (and in some countries this includes accented letters) plus ‘. with the restriction that a name must start with ‘.’ or a letter.
1. etc. You can find out how to do this by reading the manual entry for the readline library. printed (unless specifically made invisible). nor within the argument list of a function definition some of the consoles will not allow you to enter more.

> sink("record. Emacs Speaks Statistics) for working interactively with R. The collection of objects currently stored is called the workspace. foo. If you indicate that you want to do this. At the end of each R session you are given the opportunity to save all the currently available objects. ‘record. the objects are written to a file called ‘. During an R session. objects are created and stored by name (we discuss this process in the next session). It is recommended that you should use separate working directories for analyses conducted with R. say ‘commands. junk. z. The function sink. See Section “R and Emacs” in The R statistical system FAQ. When R is started at later time from the same directory it reloads the workspace from this file. ink. but it can be quite hard to decide what they might be when the several analyses have been conducted in the same directory. functions.
.Chapter 1: Introduction and preliminaries
6
Alternatively. At the same time the associated commands history is reloaded.RData’5 in the current directory.R’ in the working directory ‘work’. The leading “dot” in this file name makes it invisible in normal file listings in UNIX. To remove objects the function rm is available: > rm(x.11 Data permanency and removing objects
The entities that R creates and manipulates are known as objects.
1.
1. character strings. These may be variables. y. Names like this are often meaningful in the context of a single analysis.Rhistory’.
4 5
of unlimited length. they may be executed at any time in an R session with the command > source("commands. The R command > objects() (alternatively. bar) All objects created during an R sessions can be stored permanently in a file for use in future R sessions.lis") will divert all subsequent output from the console to an external file. temp. ls()) can be used to display the names of (most of) the objects which are currently stored within R. The command > sink() restores it to the console once again. and the command lines used in the session are saved to a file called ‘. the Emacs text editor provides more general support mechanisms (via ESS. arrays of numbers.lis’. or more general structures built from such components.10 Executing commands from or diverting output to a file
If commands4 are stored in an external file.R") For Windows Source is also available on the File menu. It is quite common for objects with names x and y to be created during an analysis.

c(10. 0. the value of the expression is a vector with the same length as the longest vector which occurs in the expression. Shorter vectors in the expression are recycled as often as need be (perhaps fractionally) until they match the length of the longest vector. the action of c() is rather different. 6.4. <-.2 Vector arithmetic
Vectors can be used in arithmetic expressions.4 and 21.1.4.1 [Concatenating lists].Last.1 Vectors and assignment
R operates on named data structures.6. can be thought of as a syntactic short-cut to this.6. So now if we were to use the command > 1/x the reciprocals of the five values would be printed at the terminal (and the value of x. the value is printed and lost 2 . 6. So with the above assignments the command
1 2
With other than vector types of argument. 3. In particular a constant is simply repeated.6.4.value before any other statements are executed.4.
. 5. 3. 21. Vectors occurring in the same expression need not all be of the same length.4. which consists of the two characters ‘<’ (“less than”) and ‘-’ (“minus”) occurring strictly side-by-side and it ‘points’ to the object receiving the value of the expression. c(10. In most contexts the ‘=’ operator can be used as a alternative. use the R command > x <.7. x) would create a vector y with 11 entries consisting of two copies of x with a zero in the middle place. say.1 A number occurring by itself in an expression is taken as a vector of length one. unchanged). 21. numbers and vectors
2. See Section 6. 5.4. The further assignment > y <. 3.4. Actually. The simplest such structure is the numeric vector. in which case the operations are performed element by element. namely 10.1. such as list mode arguments.2. page 29. 5.
2. 6.c(x.1. of course. An equivalent way of making the same assignment as above is with: > assign("x". Assignment can also be made using the function assign(). it is still available as . 3. 21. which is a single entity consisting of an ordered collection of numbers.7) -> x If an expression is used as a complete command. If they are not. 5.1. numbers and vectors
7
2 Simple manipulations.7)) The usual operator. consisting of five numbers. Assignments can also be made in the other direction.7) This is an assignment statement using the function c() which in this context can take an arbitrary number of vector arguments and whose value is a vector got by concatenating its arguments end to end. using the obvious change in the assignment operator.Chapter 2: Simple manipulations.6. To set up a vector named x. 6. Notice that the assignment operator (‘<-’). So the same assignment could be made using > c(10.

and to many other R functions. max(x)). In addition all of the common arithmetic functions are available. and if these are the only two arguments given the result is the same as the colon operator. sort(x) returns a vector of the same size as x with the elements arranged in increasing order. For example 1:30 is the vector c(1. so. which is the same as sum(x)/length(x). all have their usual meaning. 2*x repeated 2. 4. 29. *. To work with complex numbers.3 Generating regular sequences
R has a number of facilities for generating commonly used sequences of numbers. The elementary arithmetic operators are the usual +. and prod(x) their product. Internally calculations are done as double precision real numbers. The parallel maximum and minimum functions pmax and pmin return a vector (of length equal to their longest argument) that contains in each element the largest (smallest) element in that position in any of the input vectors. .. If the argument to var() is an n-by-p matrix the value is a p-by-p sample covariance matrix got by regarding the rows as independent p-variate sample vectors.2 times.10) is the same vector as 2:10. however there are other more flexible sorting facilities available (see order() or sort. range is a function whose value is a vector of length two. but sqrt(-17+0i) will do the computations as complex numbers. The first two parameters may be
. / and ^ for raising to a power... sqrt. It has five arguments. Parameters to seq(). namely c(min(x). Two statistical functions are mean(x) which calculates the sample mean. or double precision complex numbers if the input data are complex. in which case the order in which they appear is irrelevant.list() which produce a permutation to do the sorting). supply an explicit complex part. specify the beginning and end of the sequence. For most purposes the user will not be concerned if the “numbers” in a numeric vector are integers. for example 2*1:15 is the vector c(2. The function seq() is a more general facility for generating sequences. if given. and so on. element by element. 30). 28. sum(x) gives the total of the elements in x. numbers and vectors
8
> v <. exp. max and min select the largest and smallest elements of a vector respectively. The first two arguments. can also be given in named form. 2.. and 1 repeated 11 times.. tan. The construction 30:1 may be used to generate a sequence backwards. . reals or even complex. Put n <. y repeated just once.. and var(x) which gives sum((x-mean(x))^2)/(length(x)-1) or sample variance.2*x + y + 1 generates a new vector v of length 11 constructed by adding together. Note that max and min select the largest and smallest values in their arguments. only some of which may be specified in any one call.
2. -. cos. That is seq(2. length(x) is the number of elements in x. sin.Chapter 2: Simple manipulations. even if they are given several vectors.10 and compare the sequences 1:n-1 and 1:(n-1). The colon operator has high priority within an expression. 30). log. Thus sqrt(-17) will give NaN and a warning.

.. which specify a step size and a length for the sequence respectively. then c1 & c2 is their intersection (“and”). The simplest form is > s5 <. In addition if c1 and c2 are logical expressions. The logical operators are <. 5. length(vector ).. you should always use TRUE and FALSE. the default by=1 is assumed. -4. In general any operation
. However there are situations where logical vectors and their coerced numeric counterparts are not equivalent. thus seq(1.
2. For example > temp <. The first two are often abbreviated as T and F. numbers and vectors
9
named from=value and to=value . FALSE. Logical vectors are generated by conditions. When an element or value is “not available” or a “missing value” in the statistical sense. Note however that T and F are just variables which are set to TRUE and FALSE by default. times=5) which will put five copies of x end-to-end in s5.6. . If neither of these is given. Logical vectors may be used in ordinary arithmetic. >=.. each=5) which repeats each element of x five times before moving on to the next. in which case they are coerced into numeric vectors. Similarly > s4 <. or the empty sequence if the vector is empty (as it can be).2) -> s3 generates in s3 the vector c(-5.rep(x. see below). >. and NA (for “not available”. from=-5.0). c1 | c2 is their union (“or”). by=. 4.8.seq(length=51. and !c1 is the negation of c1. The next two parameters to seq() may be named by=value and length=value . == for exact equality and != for inequality. to=30) and seq(to=30.x > 13 sets temp as a vector of the same length as x with values FALSE corresponding to elements of x where the condition is not met and TRUE where it is. but are not reserved words and hence can be overwritten by the user.
2. A related function is rep() which can be used for replicating an object in various complicated ways. FALSE becoming 0 and TRUE becoming 1. Hence. -4. .6. 2. which if used must be the only parameter.. The elements of a logical vector can have the values TRUE. a place within a vector may be reserved for it by assigning it the special value NA. and creates a sequence 1.Chapter 2: Simple manipulations.8. The fifth parameter may be named along=vector . 5. respectively. seq(from=1.4 Logical vectors
As well as numerical vectors.0.5 Missing values
In some cases the components of a vector may not be completely known.. <=. R allows manipulation of logical quantities.rep(x. for example see the next subsection. For example > seq(-5. from=1) are all the same as 1:30.30). 4. Another useful version is > s6 <. by=.2) generates the same vector in s4.

na(x) since NA is not really a value but a marker for a quantity that is not available. using \ as the escape character. NaN. Where needed they are denoted by a sequence of characters delimited by the double quote character.
2. the so-called Not a Number. The arguments are by default separated in the result by a single blank character. is.is. Examples are > 0/0 or > Inf . "X3". e. "Y10")
. The motivation for this rule is simply that if the specification of an operation is incomplete.nan(xx) is only TRUE for NaNs.na(x) gives a logical vector of the same size as x with value TRUE if and only if the corresponding element in x is NA. in the same way they would be if they were printed. "New iteration results". 1:10. Note that there is a second kind of “missing” values which are produced by numerical computation.6 Character vectors
Character quantities and character vectors are used frequently in R. The function is. In summary. sep="") makes labs into the character vector c("X1". "Y4". "Y6". which changes it to string . Other useful escape sequences are \n. ind <. Character vectors may be concatenated into a vector by the c() function. "X9". Thus x == NA is a vector of the same length as x all of whose values are NA as the logical expression itself is incomplete and hence undecidable."Y"). backspace—see ?Quotes for a full list. but are printed using double quotes (or sometimes without quotes).NA). so \\ is entered and printed as \\. numbers and vectors
10
on an NA becomes an NA. To differentiate these.na(xx) is TRUE both for NA and NaN values. examples of their use will emerge frequently. the result cannot be known and hence is not available. for example as plot labels.g. that is. "X5". Character strings are entered using either matching double (") or single (’) quotes. Missing values are sometimes printed as <NA> when character vectors are printed without quotes. "X7". \t. "Y8". and inside double quotes " is entered as \". For example > labs <. newline.Inf which both give NaN since the result cannot be defined sensibly. but this can be changed by the named parameter. possibly empty.Chapter 2: Simple manipulations. is.c(1:3..paste(c("X". values. tab and \b. sep=string . "x-values".na(z) Notice that the logical expression x == NA is quite different from is. > z <. They use C-style escape sequences. Any numbers given among the arguments are coerced into character strings in the evident way. "Y2". The paste() function takes an arbitrary number of arguments and concatenates them one by one into character strings.

4.7 Index vectors. times=4)] (an admittedly unlikely thing to do) produces a character vector of length 16 consisting of "x". . A vector of character strings."orange")]
3
paste(. "peach") > lunch <. 2. "x" repeated four times. > fruit <. In this case a sub-vector of the names vector may be used in the same way as the positive integral labels in item 2 further above. 1. Note that if x has missing values.c(5. In this case the index vector must be of the same length as the vector from which elements are to be selected..x[-(1:5)] gives y all but the first five elements of x. In this case the values in the index vector must lie in the set {1. For example > y <.1).3
2. in that order. in the same order. collapse=ss ) joins the arguments into a single character string putting ss in between.. A logical vector.na(x)) & x>0] -> z creates an object z and places in it the values of the vector x+1 for which the corresponding value in x was both non-missing and positive. The corresponding elements of the vector are selected and concatenated. y will be shorter than x. Such an index vector specifies the values to be excluded rather than included. numbers and vectors
11
Note particularly that recycling of short lists takes place here too. A vector of positive integral quantities. Thus > y <..2. . Also > c("x". 10. A vector of negative integral quantities. This possibility only applies where an object has a names attribute to identify its components. The index vector can be of any length and the result is of the same length as the index vector. length(x)}. selecting and modifying subsets of a data set
Subsets of the elements of a vector may be selected by appending to the name of the vector an index vector in square brackets. Such index vectors can be any of four distinct types. see the help for sub and substring.
. thus c("X". "apple". 1.c("orange". .na(x)] creates (or re-creates) an object y which will contain the non-missing values of x. 20) > names(fruit) <.Chapter 2: Simple manipulations. in the result. For example x[6] is the sixth component of x and > x[1:10] selects the first 10 elements of x (assuming length(x) is not less than 10). "banana". 3. 2. "y". Also > (x+1)[(!is. More generally any expression that evaluates to a vector may have subsets of its elements similarly selected by appending an index vector in square brackets immediately after the expression.2.x[!is. "Y") is repeated 5 times to match the sequence 1:10."y")[rep(c(1.fruit[c("apple". "y". . Values corresponding to TRUE in the index vector are selected and those corresponding to FALSE are omitted. There are more tools for character manipulation.

See Chapter 4 [Factors].Chapter 2: Simple manipulations. See Section 6. See Section 6. This option is particularly useful in connection with data frames. in which case the assignment operation is performed only on those elements of the vector. • data frames are matrix-like structures. • matrices or more generally arrays are multi-dimensional generalizations of vectors. but there are several others which we will meet more formally in later sections.8 Other types of objects
Vectors are the most important type of object in R. numbers and vectors
12
The advantage is that alphanumeric names are often easier to remember than numeric indices. • functions are themselves objects in R which can be stored in the project’s workspace.na(x)] <.1 [Lists].abs(y)
2. they are vectors that can be indexed by two or more indices and will be printed in special ways. page 16. Many experiments are best described by data frames: the treatments are categorical but the response is numeric. as we shall see later. See Chapter 10 [Writing your own functions]. Lists provide a convenient way to return the results of a statistical computation. This provides a simple and convenient way to extend R. • factors provide compact ways to handle categorical data. For example > x[is. See Chapter 5 [Arrays and matrices]. The vector assigned must match the length of the index vector. and in the case of a logical index vector it must again be the same length as the vector it is indexing.0 replaces any missing values in x by zeros and > y[y < 0] <. Think of data frames as ‘data matrices’ with one row per observational unit but with (possibly) both numerical and categorical variables. page 19.
. An indexed expression can also appear on the receiving end of an assignment. In fact.3 [Data frames].-y[y < 0] has the same effect as > y <. page 28. page 44. • lists are a general form of vector in which the various elements need not be of the same type. in which the columns can be of different types. and are often themselves vectors or lists. The expression must be of the form vector[index_vector ] as having an arbitrary expression in place of the vector name does not make much sense here. page 29.

see Section 3.character(z) after which digits is the character vector c("0". Another property of every object is its length.g.. character or raw.0:9 we could put > digits <. (The only apparent exception to this rule is the special “value” listed as NA for quantities not available. By the mode of an object we mean the basic type of its fundamental constituents. A further coercion. logical.. (and a few where it might not be). their modes and attributes
13
3 Objects. which are of mode list. their modes and attributes
3. Expressions as objects form an advanced part of R which will not be discussed in this guide. R also operates on objects called lists. complex.as. Because of this. Vectors must have their values all of the same mode.
.1 Intrinsic attributes: mode and length
The entities R operates on are technically known as objects.Chapter 3: Objects. reconstructs the numerical vector again:
1 2
numeric mode is actually an amalgam of two distinct modes. "9"). "1". For example with > z <. or mode. page 14. then in an expression mode(z) is the character string "complex" and length(z) is 100. R caters for changes of mode almost anywhere it could be considered sensible to do so. except indirectly when we discuss formulae used with modeling in R. if z is a complex vector of length 100. or change of mode.3 [Getting and setting attributes]. For example the empty character string vector is listed as character(0) and the empty numeric vector as numeric(0). These are ordered sequences of objects which individually can be of any mode. lists are known as “recursive” rather than atomic structures since their components can themselves be lists in their own right. vectors of logical values and vectors of character strings. Further properties of an object are usually provided by attributes(object ). namely numeric 1 . Note however that length(object ) does not always contain intrinsic useful information. numeric.. e. which we discuss in some detail later. Thus any given vector must be unambiguously either logical. mode and length are also called “intrinsic attributes” of an object. Note that a vector can be empty and still have a mode. when object is a function. For example. This is a special case of a “property” of an object. .. The other recursive structures are those of mode function and expression. Examples are vectors of numeric (real) or complex values. "2". namely integer and double precision. as explained in the manual. complex. Functions are the objects that form part of the R system along with similar user written functions. character and raw. These are known as “atomic” structures since their components are all of the same type. but in fact there are several types of NA). The functions mode(object ) and length(object ) can be used to find out the mode and length of any defined structure2 .

(the first two components of which are at this point both NA). This applies to any structure at all. their modes and attributes
14
> d <. "dim") <. however. The function attr(object. because of roundoff errors in the character representation. Similarly character() is a empty character vector.3 and vectors can be extended (by missing values) in the same way. This automatic adjustment of lengths of an object is used often. (see Section 7. Once an object of any size has been created.
3.3 There is a large collection of functions of the form as.Chapter 3: Objects. except in rather special circumstances when some new attribute is being created for some particular purpose. coercion from numeric to character and back again will not be exactly reversible. is very important.integer(digits) Now d and z are the same.) We can then retain just the first three values by > length(alpha) <. (The old indices are not retained. Some care should be exercised when assigning or deleting attributes since they are an integral part of the object system used in R. then > alpha <. page 33. Thus > e[3] <. The concept.as. When it is used on the left hand side of an assignment it can be used either to associate a new attribute with object or to change an existing one. For example > e <.
. or for investing an object with some other attribute it may not already possess. These functions are rarely used.
3
In general.numeric() makes e an empty vector structure of mode numeric. and so on.17 now makes e a vector of length 3. provided the mode of the additional component(s) agrees with the mode of the object in the first place. of course. Hence if alpha is an object of length 10.alpha[2 * 1:5] makes it an object of length 5 consisting of just the former components with even index.c(10.3 Getting and setting attributes
The function attributes(object ) returns a list of all the non-intrinsic attributes currently defined for that object.something () for either coercion from one mode to another.
3. name ) can be used to select a specific attribute.10) allows R to treat z as if it were a 10-by-10 matrix. The reader should consult the different help files to become familiar with them. For example > attr(z. new components may be added to it simply by giving it an index value outside its previous range.2 Changing the length of an object
An “empty” object may still have a mode. for example in the scan() function for input.2 [The scan() function]. for example to associate a creation date or an operator with an R object.) Conversely to truncate the size of an object requires only an assignment to do so.

and other so-called generic functions such as summary() will react to it as an argument in a way sensitive to its class. but only briefly. page 51. but one is when you are learning to come to terms with the idea of class and generic functions. "array".Chapter 3: Objects. but "matrix". the plot() function will display it graphically in a certain way. which is rather like a matrix.
4
A different style using ‘formal’ or ‘S4’ classes is provided in package methods. their modes and attributes
15
3. For example if winter has the class "data. Only in rather special situations do you need to use this facility. For simple vectors this is just the mode. "logical". it will be printed in a certain way. use the function unclass(). for example "numeric". reported by the function class. To remove temporarily the effects of class.
. "factor" and "data. "character" or "list". whereas > unclass(winter) will print it as an ordinary list.4 The class of an object
All objects in R have a class. A special attribute known as the class of the object is used to allow for an object-oriented style4 of programming in R.frame" then > winter will print it in data frame form.frame".frame" are other possible values. Generic functions and classes will be discussed further in Section 10.9 [Object orientation]. For example if an object has class "data.

1 [Contrasts]. 54. the Northern Territory.
.500 57. Tasmania. "qld". "vic". South Australia. "vic". 59. "vic". statef. "vic". "wa". While the “real” application of factors is with model formulae (see Section 11. 60.Chapter 4: Ordered and unordered factors
16
4 Ordered and unordered factors
A factor is a vector object used to specify a discrete classification (grouping) of the components of other vectors of the same length. we have a sample of 30 tax accountants from all the states and territories of Australia1 and their individual state of origin is specified by a character vector of state mnemonics as > state <.600 55. New South Wales. A factor is similarly created using the factor() function: > statef <. 59. "sa". 42. "tas". 48. 58.500 53. 52. "nsw". "qld". 49. 58. "vic".500 56. page 56). namely the Australian Capital Territory. “sorted” means sorted in alphabetical order. "nsw".c("tas". mean) giving a means vector with the components labelled by the levels act nsw nt qld sa tas vic wa 44. "nsw". suppose we have the incomes of the same tax accountants in another vector (in suitably large units of money) > incomes <. we here look at a specific example. Victoria and Western Australia. 40. 62. 64. for example.tapply(incomes. 61. "wa".000 60. "wa". 65. > levels(statef) [1] "act" "nsw" "nt" "qld" "sa" "tas" "vic" "wa"
4. 41.1 A specific example
Suppose. "sa".000 52. 61.c(60. "sa". "qld". "nsw". 46. "qld".2 The function tapply() and ragged arrays
To continue the previous example. 61. 43) To calculate the sample mean income for each state we can now use the special function tapply(): > incmeans <. 49. 56. 46. "nt". "qld".1.
4. 69. 70. 61. "act". "nt". 49.250
1
Readers should note that there are eight states and territories in Australia.factor(state) The print() function handles factors slightly differently from other objects: > statef [1] tas sa qld nsw nsw nt wa wa qld vic nsw vic qld qld sa [16] tas sa nt wa vic qld nsw nsw wa sa act nsw vic vic act Levels: act nsw nt qld sa tas vic wa To find out the levels of a factor the function levels() can be used. "wa".333 55. "nsw". 48. "nsw". "act")
Notice that in the case of a character vector. "sa". Queensland. 51. R provides both ordered and unordered factors.

statef. page 44. The combination of a vector and a labelling factor is an example of what is sometimes called a ragged array. to each group of components of the first argument. defined by the levels of the second component.) After this assignment. or in the order they were specified to factor if they were specified explicitly. (You could also investigate R’s facilities for t-tests.g. The result is a structure of the same length as the levels attribute of the factor containing the results. here incomes.3102 4. we might wish to split the tax accountants by both state and sex. since arguments are coerced to factors when necessary (using as. The values in the vector are collected into groups corresponding to the distinct entries in the factor.3 Ordered factors
The levels of factors are stored in alphabetical order. The ordered() function creates such ordered factors but is otherwise identical to factor. specified by the assignment: > stderr <.7386 0.) The function tapply() can also be used to handle more complicated indexing of a vector by multiple categories.factor()).
4. here mean(). the standard errors are calculated by > incster <. e.function(x) sqrt(var(x)/length(x)) (Writing functions will be considered later in Chapter 10 [Writing your own functions]. state)’. When the subclass sizes are all the same the indexing may be done implicitly and much more efficiently. stderr) and the values calculated are then > incster act nsw nt qld sa tas vic wa 1. For most purposes the only difference between ordered
2
Note that tapply() also works in this case when its second argument is not a factor. and this is true for quite a few other functions. For example. To do this you could use tapply() once more with the length() function to find the sample sizes. such a function is a very simple one liner. and in this case was unnecessary as R also has a builtin function sd().5 4.. and the qt() function to find the percentage points of the appropriate t-distributions. Sometimes the levels will have a natural ordering that we want to record and want our statistical analysis to make use of. However in this simple instance (just one factor) what happens can be thought of as follows. The value is a vector of function results.5 5. here statef2 . To do this we need to write an R function to calculate the standard error for any given vector.5 4.
. Since there is an builtin function var() to calculate the sample variance. Suppose further we needed to calculate the standard errors of the state income means.1061 2. since the subclass sizes are possibly irregular. as if they were separate vector structures. as we see in the next section. ‘tapply(incomes.6575 As an exercise you may care to find the usual 95% confidence limits for the state mean incomes.244 2. The reader should consult the help document for more details. labelled by the levels attribute of the factor. The function is then applied to each of these groups individually.Chapter 4: Ordered and unordered factors
17
The function tapply() is used to apply a function.tapply(incomes.

but the contrasts generated for them in fitting linear models are different.
.Chapter 4: Ordered and unordered factors
18
and unordered factors is that the former are printed showing the ordering of the levels.

4. if an array name is given with just one subscript or index vector. subsections of an array may be specified by giving a sequence of index vectors in place of subscripts. in this case the dimension vector is ignored.” with the first subscript moving fastest and the last subscript slowest. z is a vector of 1500 elements. then the full range of that subscript is taken. a[2. Also.2) and data vector containing the values c(a[2. This is not the case.2].2]. a[2. for example numeric. Arrays can be one-dimensional: such arrays are usually treated in the same way as vectors (including when printing).2].1]. a[2.
5. Suppose. as we shall see in Section 5. however. page 21.1].3.. separated by commas. e.2 Array indexing. a[2. a[2. say a.100) gives it the dim attribute that allows it to be treated as a 3 by 5 by 100 array..4.2. a[2. and in particular the special case of matrices. that is “column major order.1 Arrays
An array can be considered as a multiply subscripted collection of data entries.4.2) then there are 3×4×2 = 24 entries in a and the data vector holds them in the order a[1. however if any index position is given an empty index vector.3. .. which is the same as omitting the subscripts entirely and using a alone. A dimension vector is a vector of non-negative integers.4.4 [The array() function]. a[3.1].5.
.Chapter 5: Arrays and matrices
19
5 Arrays and matrices
5.2]) in that order. If its length is k then the array is k-dimensional.c(3...1. For example if the dimension vector for an array.1. For any array.4. if the single index is not a vector but itself an array. a[2.1. A vector can be used by R as an array only if it has a dimension vector as its dim attribute.1].2]. a[2. say Z. R allows simple facilities for creating and handling arrays. is c(3. Continuing the previous example. a[2.g. a[2.1. for example.] is a 4 × 2 array with dimension vector c(4. The values in the data vector give the values in the array in the same order as they would occur in FORTRAN. More generally. a[.1].] stands for the entire array. a matrix is a 2-dimensional array. but the exceptions can cause confusion.2. the dimension vector may be referenced explicitly as dim(Z) (on either side of an assignment). then the corresponding values of the data vector only are used. The dimensions are indexed from one up to the values given in the dimension vector. as we next discuss. The assignment > dim(z) <. Other functions such as matrix() and array() are available for simpler and more natural looking assignments.2].1]. Subsections of an array
Individual elements of an array may be referenced by giving the name of the array followed by the subscripts in square brackets.

array(c(1:3.Chapter 5: Arrays and matrices
20
5.3 Index matrices
As well as an index vector in any subscript position.1] [.] 1 5 0 13 17 [2. dim=c(4.4] [.] 2 0 10 14 18 [3. We could proceed as follows: > > > > Xb Xv ib iv <<<<matrix(0. The entries in the index matrix are the row and column indices for the doubly indexed array.0 # Replace those elements by zeros. n. n. > x [.2] [. varieties)
.] 0 7 11 15 19 [4.1] [. as in the following example. v) cbind(1:n.] 1 5 9 13 17 [2. suppose we wish to generate an (unreduced) design matrix for a block design defined by factors blocks (b levels) and varieties (v levels).3].2] [. and • Replace these entries in the array X by zeroes. Suppose for example we have a 4 by 5 array X and we wish to do the following: • Extract elements X[1. In this case we need a 3 by 2 subscript array.2] [1.4] [.] 2 6 10 14 18 [3. > x [. > x <. [. blocks) cbind(1:n. As a less trivial example.] 3 7 11 15 19 [4.1] as a vector structure.] 2 2 [3. or to extract an irregular collection as a vector.] 1 3 [2. an index matrix may be given consisting of two columns and as many rows as desired.3:1). In the case of a doubly indexed array.3] [. dim=c(3.2)) > i # i is a 3 by 2 index array.1] [.] 4 8 12 16 20 > Negative indices are not allowed in index matrices.5] [1.array(1:20.2] and X[3. Further suppose there are n plots in the experiment.] 4 8 12 16 20 > i <.3] [.5] [1.5)) # Generate a 4 by 5 array. a matrix may be used with a single index matrix in order either to assign a vector of quantities to an irregular collection of elements in the array. b) matrix(0. A matrix example makes the process clear. NA and zero values are allowed: rows in the index matrix containing a zero are ignored. X[2. and rows containing an NA produce an NA in the result.] 3 1 > x[i] # Extract those elements [1] 9 6 3 > x[i] <.

Chapter 5: Arrays and matrices

21

> Xb[ib] <- 1 > Xv[iv] <- 1 > X <- cbind(Xb, Xv) To construct the incidence matrix, N say, we could use > N <- crossprod(Xb, Xv) However a simpler direct way of producing this matrix is to use table(): > N <- table(blocks, varieties) Index matrices must be numerical: any other form of matrix (e.g. a logical or character matrix) supplied as a matrix is treated as an indexing vector.

5.4 The array() function
As well as giving a vector structure a dim attribute, arrays can be constructed from vectors by the array function, which has the form > Z <- array(data_vector, dim_vector ) For example, if the vector h contains 24 or fewer, numbers then the command > Z <- array(h, dim=c(3,4,2)) would use h to set up 3 by 4 by 2 array in Z. If the size of h is exactly 24 the result is the same as > Z <- h ; dim(Z) <- c(3,4,2) However if h is shorter than 24, its values are recycled from the beginning again to make it up to size 24 (see Section 5.4.1 [The recycling rule], page 21) but dim(h) <- c(3,4,2) would signal an error about mismatching length. As an extreme but common example > Z <- array(0, c(3,4,2)) makes Z an array of all zeros. At this point dim(Z) stands for the dimension vector c(3,4,2), and Z[1:24] stands for the data vector as it was in h, and Z[] with an empty subscript or Z with no subscript stands for the entire array as an array. Arrays may be used in arithmetic expressions and the result is an array formed by element-by-element operations on the data vector. The dim attributes of operands generally need to be the same, and this becomes the dimension vector of the result. So if A, B and C are all similar arrays, then > D <- 2*A*B + C + 1 makes D a similar array with its data vector being the result of the given element-by-element operations. However the precise rule concerning mixed array and vector calculations has to be considered a little more carefully.

5.4.1 Mixed vector and array arithmetic. The recycling rule
The precise rule affecting element by element mixed calculations with vectors and arrays is somewhat quirky and hard to find in the references. From experience we have found the following to be a reliable guide. • The expression is scanned from left to right.

Chapter 5: Arrays and matrices

22

• Any short vector operands are extended by recycling their values until they match the size of any other operands. • As long as short vectors and arrays only are encountered, the arrays must all have the same dim attribute or an error results. • Any vector operand longer than a matrix or array operand generates an error. • If array structures are present and no error or coercion to vector has been precipitated, the result is an array structure with the common dim attribute of its array operands.

5.5 The outer product of two arrays
An important operation on arrays is the outer product. If a and b are two numeric arrays, their outer product is an array whose dimension vector is obtained by concatenating their two dimension vectors (order is important), and whose data vector is got by forming all possible products of elements of the data vector of a with those of b. The outer product is formed by the special operator %o%: > ab <- a %o% b An alternative is > ab <- outer(a, b, "*") The multiplication function can be replaced by an arbitrary function of two variables. For example if we wished to evaluate the function f (x; y) = cos(y)/(1 + x2 ) over a regular grid of values with x- and y-coordinates defined by the R vectors x and y respectively, we could proceed as follows: > f <- function(x, y) cos(y)/(1 + x^2) > z <- outer(x, y, f) In particular the outer product of two ordinary vectors is a doubly subscripted array (that is a matrix, of rank at most 1). Notice that the outer product operator is of course non-commutative. Defining your own R functions will be considered further in Chapter 10 [Writing your own functions], page 44.

An example: Determinants of 2 by 2 single-digit matrices
As an artificial but cute example, consider the determinants of 2 by 2 matrices [a, b; c, d] where each entry is a non-negative integer in the range 0, 1, . . . , 9, that is a digit. The problem is to find the determinants, ad − bc, of all possible matrices of this form and represent the frequency with which each value occurs as a high density plot. This amounts to finding the probability distribution of the determinant if each digit is chosen independently and uniformly at random. A neat way of doing this uses the outer() function twice: > d <- outer(0:9, 0:9) > fr <- table(outer(d, d, "-")) > plot(as.numeric(names(fr)), fr, type="h", xlab="Determinant", ylab="Frequency") Notice the coercion of the names attribute of the frequency table to numeric in order to recover the range of the determinant values. The “obvious” way of doing this problem with for loops, to be discussed in Chapter 9 [Loops and conditional execution], page 42, is so inefficient as to be impractical.

Chapter 5: Arrays and matrices

23

It is also perhaps surprising that about 1 in 20 such matrices is singular.

5.6 Generalized transpose of an array
The function aperm(a, perm) may be used to permute an array, a. The argument perm must be a permutation of the integers {1, . . . , k}, where k is the number of subscripts in a. The result of the function is an array of the same size as a but with old dimension given by perm[j] becoming the new j-th dimension. The easiest way to think of this operation is as a generalization of transposition for matrices. Indeed if A is a matrix, (that is, a doubly subscripted array) then B given by > B <- aperm(A, c(2,1)) is just the transpose of A. For this special case a simpler function t() is available, so we could have used B <- t(A).

5.7 Matrix facilities
As noted above, a matrix is just an array with two subscripts. However it is such an important special case it needs a separate discussion. R contains many operators and functions that are available only for matrices. For example t(X) is the matrix transpose function, as noted above. The functions nrow(A) and ncol(A) give the number of rows and columns in the matrix A respectively.

5.7.1 Matrix multiplication
The operator %*% is used for matrix multiplication. An n by 1 or 1 by n matrix may of course be used as an n-vector if in the context such is appropriate. Conversely, vectors which occur in matrix multiplication expressions are automatically promoted either to row or column vectors, whichever is multiplicatively coherent, if possible, (although this is not always unambiguously possible, as we see later). If, for example, A and B are square matrices of the same size, then > A * B is the matrix of element by element products and > A %*% B is the matrix product. If x is a vector, then > x %*% A %*% x is a quadratic form.1 The function crossprod() forms “crossproducts”, meaning that crossprod(X, y) is the same as t(X) %*% y but the operation is more efficient. If the second argument to crossprod() is omitted it is taken to be the same as the first. The meaning of diag() depends on its argument. diag(v), where v is a vector, gives a diagonal matrix with elements of the vector as the diagonal entries. On the other hand
1

Note that x %*% x is ambiguous, as it could mean either xT x or xxT , where x is the column form. In such cases the smaller matrix seems implicitly to be the interpretation adopted, so the scalar xT x is in this case the result. The matrix xxT may be calculated either by cbind(x) %*% x or x %*% rbind(x) since the result of rbind() or cbind() is always a matrix. However, the best way to compute xT x or xxT is crossprod(x) or x %o% x respectively.

> solve(A.
2
Even better would be to form a matrix square root B with A = BB T and find the squared length of the solution of By = x. formally x = A−1 b where A−1 denotes the inverse of A. Numerically.
5. This consists of a matrix of orthonormal columns U with the same column space as M.values = TRUE)$values
5.eigen(Sm) will assign this list to ev. The result of this function is a list of two components named values and vectors.eigen(Sm.7.b) solves the system.b). The assignment > ev <. with their names. If the expression > eigen(Sm) is used by itself as a command the two components are printed.solve(A) %*% b instead of solve(A. D is actually returned as a vector of the diagonal elements. with evident meanings.eigen(Sm)$values evals now holds the vector of eigenvalues and the second component is discarded. The result of svd(M) is actually a list of three components named d. returning x (up to some accuracy loss). somewhat confusingly. For large matrices it is better to avoid computing the eigenvectors if they are not needed by using the expression > evals <. In R. should be computed by something like2 x %*% solve(A.4 Singular value decomposition and determinants
The function svd(M) takes an arbitrary matrix argument. it is both inefficient and potentially unstable to compute x <. This is the same convention as that used for diag() in Matlab. When after > b <. The quadratic form xT A−1 x which is used in multivariate computations.Chapter 5: Arrays and matrices
24
diag(M). Had we only needed the eigenvalues we could have used the assignment: > evals <.3 Eigenvalues and eigenvectors
The function eigen(Sm) calculates the eigenvalues and eigenvectors of a symmetric matrix Sm. and calculates the singular value decomposition of M. perhaps using the Cholesky or eigendecomposition of A.7. where M is a matrix.A %*% x only A and b are given. Then ev$val is the vector of eigenvalues of Sm and ev$vec is the matrix of corresponding eigenvectors.x). Also. if k is a single numeric value then diag(k) is the k by k identity matrix!
5. Note that in linear algebra. the vector x is the solution of that linear equation system. a second matrix of orthonormal columns V whose column space is the row space of M and a diagonal matrix of positive entries D such that M = U %*% D %*% t(V). M. gives the vector of main diagonal entries of M. u and v. only.
.7.2 Linear equations and inversion
Solving linear equations is the inverse of matrix multiplication. which can be computed by solve(A) but rarely is needed. rather than computing the inverse of A.

and another.5 Least squares fitting and the QR decomposition
The function lsfit() returns a list giving results of a least squares fitting procedure. that is.2 [Linear models].
5.diag() for. you might like to consider writing a function. page 54. determinant. Another closely related function is qr() and its allies. See the help facility for more details. or column-wise.
5. and rbind() vertically. say tr(). then. [Hint: You will not need to use an explicit loop. Further note that you almost always will prefer using lm(. .. y) res <.qr.. matrices can be built up from other vectors and matrices by the functions cbind() and rbind(). y) gives the results of a least squares fit where y is the vector of observations and X is the design matrix. An assignment such as > ans <.) (see Section 11. cbind() and rbind()
As we have already seen informally.lsfit(X.8 Forming partitioned matrices. As a further trivial but potentially useful example. to calculate the trace of a square matrix. to give the sign and modulus (optionally on log scale). including the sign. page 57) to lsfit() for regression modelling. regression diagnostics. y) fit <.coef(Xplus.7. y)
These compute the orthogonal projection of y onto the range of X in fit. low-level way to perform least squares calculations. Redundancies will be discovered and removed as they are found. and also for the follow-up function ls.function(M) prod(svd(M)$d) after which we could use absdet() as just another R function. arg_2. Look again at the diag() function. it is not hard to see that > absdetM <.] R has a builtin function det to calculate a determinant.qr(X) b <. b is essentially the result of the Matlab ‘backslash’ operator.Chapter 5: Arrays and matrices
25
If M is in fact square. Although still useful in some contexts. the projection onto the orthogonal complement in res and the coefficient vector for the projection in b. It is not assumed that X has full column rank.fitted(Xplus.qr. If this calculation were needed often with a variety of matrices it could be defined as an R function > absdet <.prod(svd(M)$d) calculates the absolute value of the determinant of M. or row-wise. arg_3. among other things. it would now generally be replaced by the statistical models features.)
. Consider the following assignments > > > > Xplus <.cbind(arg_1. Note that a grand mean term is automatically included and need not be included explicitly as a column of X. Roughly cbind() forms matrices by binding together matrices horizontally.qr. This alternative is the older. as will be discussed in Chapter 11 [Statistical models in R]. In the assignment > X <.resid(Xplus.

statef. forming the columns. but more convenient than. for example. In this case any vector argument.cbind(1. for example with the cut() function:
. the basic c() function does not. are of course taken as row vectors. The function table() allows frequency tables to be calculated from equal length factors. .
5. The assignment > statefr <.table(statef) gives in statefr a table of frequencies of each state in the sample. This is occasionally useful in its own right.
5. possibly cyclically extended. This simple case is equivalent to. The frequencies are ordered and labelled by the levels attribute of the factor. Suppose. If some of the arguments to cbind() are vectors they may be shorter than the column size of any matrices present. arg 2.vector() > vec <.9 The concatenation function.c(X) There are slight differences between the two. but ultimately the choice between them is largely a matter of style (with the former being preferable). that statef is a factor giving the state code for each entry in a data vector. X1.vector(X) However a similar result can be achieved by using c() with just one argument. and so on. together with an initial column of 1s we can use > X <.tapply(statef. simply for this side-effect: > vec <. . the result is a k-way array of frequencies. or matrices with the same column size. X2) The result of rbind() or cbind() always has matrix status.as.Chapter 5: Arrays and matrices
26
the arguments to cbind() must be either vectors of any length. but rather clears numeric objects of all dim and dimnames attributes. length) Further suppose that incomef is a factor giving a suitably defined “income class” for each entry in the data vector. If there are k factor arguments. . in which case they are cyclically extended to match the matrix column size (or the length of the longest vector if no matrices are given). To combine these by columns into a matrix X. with arrays
It should be noted that whereas cbind() and rbind() are concatenation functions that respect dim attributes. c(). The function rbind() does the corresponding operation for rows. > statefr <. The official way to coerce an array back to a simple vector object is to use as. Suppose X1 and X2 have the same number of rows.10 Frequency tables from factors
Recall that a factor defines a partition into groups. Hence cbind(x) and rbind(x) are possibly the simplest ways explicitly to allow the vector x to be treated as a column or row matrix respectively. that is the same number of rows. Similarly a pair of factors defines a two way cross classification. The result is a matrix with the concatenated arguments arg 1.

the names are transferred to the sublist. Lst[[3]] and Lst[[4]].1 Lists
An R list is an object consisting of an ordered collection of objects known as its components. Here is a simple example of how to make a list: > Lst <. This is especially useful. then the function length(Lst) gives the number of (top level) components it has. Lst$child. Lst[[4]] is a vector subscripted array then Lst[[4]][1] is its first entry. one can also use the names of the list components in double square brackets.e.]’ is a general subscripting operator. So in the simple example given above: Lst$name is the same as Lst[[1]] and is the string "Fred". no."name".ages[1] is the same as Lst[[4]][1] and is the number 4.ages=c(4.list(name="Fred". If. for example.
.9)) Components are always numbered and may always be referred to as such. This is a very useful convention as it makes it easier to get the right component if you forget the number. further. The names of components may be abbreviated down to the minimum number of letters needed to identify them uniquely. i. a complex vector. a character array. Lst[[x]] It is very important to distinguish Lst[[1]] from Lst[1]. Lst$wife is the same as Lst[[2]] and is the string "Mary". Lst[["name"]] is the same as Lst$name.7.. The vector of names is in fact simply an attribute of the list like any other and may be handled as such. If Lst is a list. more conveniently. a matrix. child. a logical value.Chapter 6: Lists and data frames
28
6 Lists and data frames
6. Additionally. when the name of the component to be extracted is stored in another variable as in > x <. and if it is a named list the name is not included. or. a function.. and so on.. Other structures besides lists may. The latter is a sublist of the list Lst consisting of the first entry only. Thus the former is the first object in the list Lst.. similarly be given a names attribute also. Lst[[2]]. of course. ‘[[. Components of lists may also be named. If it is a named list. wife="Mary".children=3. There is no particular need for the components to be of the same mode or type. Thus if Lst is the name of a list with four components.. these may be individually referred to as Lst[[1]]. whereas ‘[. by giving an expression of the form > name $component_name for the same thing. and in this case the component may be referred to either by giving the component name as a character string in place of the number in double square brackets. Thus Lst$coefficients may be minimally specified as Lst$coe and Lst$covariance as Lst$cov. a list could consist of a numeric vector.]]’ is the operator used to select a single element. and.

frame(home=statef. and matrix structures must all have the same row size.2 Constructing and modifying lists
New lists may be formed from existing objects by the function list().data.B. . lists. (which can be freely chosen). There are restrictions on lists that may be made into data frames.list(name_1 =object_1.Chapter 6: Lists and data frames
29
6. • Matrices. or variables.. the result is an object of mode list also.. > list. or other data frames.
6. shot=incomef) A list whose components conform to the restrictions of a data frame may be coerced into a data frame using the function as. For example > Lst[5] <. and data frames provide as many variables to the new data frame as they have columns. whose levels are the unique values appearing in the vector. numeric matrices. such as dim attributes. An assignment of the form > Lst <. If these names are omitted. A data frame may for many purposes be regarded as a matrix with columns possibly of differing modes and attributes. object m for the components and giving them names as specified by the argument names.. • Numeric vectors.frame". whose components are those of the argument lists joined together in sequence.A. list.list(matrix=Mat)
6. .1 Making data frames
Objects satisfying the restrictions placed on the columns (components) of a data frame may be used to form one using the function data. elements. • Vector structures appearing as variables of the data frame must all have the same length. logicals and factors are included as is. lists. Lists. The components used to form the list are copied when forming the new list and the originals are not affected.3.ABC <. can be extended by specifying additional components. .data.frame: > accountants <.3 Data frames
A data frame is a list with class "data. respectively.frame()
. In this case all other attributes. . are discarded. factors.2. character. It may be displayed in matrix form.
6. namely • The components must be vectors (numeric. like any subscripted object.c(list. and its rows and columns extracted using matrix indexing conventions. name_m =object_m ) sets up a list Lst of m components using object 1. or logical).1 Concatenating lists
When the concatenation function c() is given list arguments. loot=incomes. the components are numbered only.C) Recall that with vector objects as arguments the concatenation function similarly joined together all arguments into a single vector structure. list. and character vectors are coerced to be factors. .

without the need to quote the list name explicitly each time. The attach() function takes a ‘database’ such as a list or data frame as its argument. this statement detaches from the search path the entity currently at position 2.3. but rather masks it with another variable u in the working directory at position 1 on the search path. for list components is not always very convenient. To detach a data frame.
6. At this point an assignment such as > u <. and what is attached is a copy of the original object.2 attach() and detach()
The $ notation. except under the list notation as lentils$u and so on. v or w in position 1.v+w However the new value of component u is not visible until the data frame is detached and attached again. lentils$v. and provided there are no variables u. v and w would be no longer visible. such as accountants$statef. but the original list or data frame is unchanged.3 Working with data frames
A useful convention that allows you to work with many different problems comfortably together in the same working directory is • gather together all variables for any well defined and separate problem in a data frame under a suitably informative name. The attach > attach(lentils) places the data frame in the search path at position 2. Thus in the present context the variables u. the simplest way is to resort once again to the $ notation: > lentils$u <. A useful facility would be somehow to make the components of a list or data frame temporarily visible as variables under their component name.Chapter 6: Lists and data frames
30
The simplest way to construct a data frame from scratch is to use the read. • when working with a problem attach the appropriate data frame at position 2. To make a permanent change to the data frame itself.v+w does not replace the component u of the data frame. page 32. but it is much safer to always use a name. for example by detach(lentils) or detach("lentils") Note: In R lists and data frames can only be attached at position 2 or above.3. and then detach(). use the function > detach() More precisely.
6. You can alter the attached values via assign. • before leaving a problem. Thus suppose lentils is a data frame with three variables lentils$u. and use the working directory at level 1 for operational quantities and temporary variables.
. Entities at positions greater than 2 on the search path can be detached by giving their number to detach. u. lentils$w. This is discussed further in Chapter 7 [Reading data from files]. v and w are available as variables from the data frame in their own right. add any variables you wish to keep for future reference to the data frame using the $ form of assignment.table() function to read an entire data frame from an external file.

5 Managing the search path
The function search shows the current search path and so is a very useful way to keep track of which data frames and lists (and packages) have been attached and detached. Initially it gives > search() [1] ".GlobalEnv" "Autoloads" "package:base" 1 where .
6. In particular any object of mode "list" may be attached in the same way: > attach(any.
. > detach("lentils") > search() [1] ".old.Chapter 6: Lists and data frames
31
• finally remove all unwanted variables from the working directory and keep it as clean of left-over temporary variables as possible.list) Anything that has been attached can be detached by detach. y and z.
6.4 Attaching arbitrary lists
attach() is a generic function that allows not only directories and data frames to be attached to the search path.3. for example. all of which have variables named x. In this way it is quite simple to work with many problems in the same directory. we detach the data frame and confirm it has been removed from the search path. by position number or.GlobalEnv" "lentils" "Autoloads" "package:base" > ls(2) [1] "u" "v" "w" and as we see ls (or objects) can be used to examine the contents of any position on the search path. preferably.3. After lentils is attached we have > search() [1] ". Finally.GlobalEnv" "Autoloads" "package:base"
1
See the on-line help for autoload for the meaning of the second term. by name.GlobalEnv is the workspace. but other classes of object as well.

use the package argument. edit brings up a separate spreadsheet-like environment for editing. However. for example data(package="rpart") data(Puromycin.Chapter 7: Reading data from files
34
data() As from R version 2.frame()) to enter new data via the spreadsheet interface. In most cases this will load an R object of the same name. its datasets are automatically included in the search. which is equivalent to xold <. This is useful for making small changes once a data set has been read. for example data(infert) and this can still be used with the standard packages (as in this example).1 Loading data from other R packages
To access data from a particular package.0.edit(xold).3.
7. The command > xnew <. many packages still use the earlier convention in which data was also used to load datasets into R. and on completion the changed object is assigned to xnew.0 all the datasets supplied with R are available directly by name.edit(data. package="datasets") If a package has been attached by library. If you want to alter the original dataset xold.edit(xold) will allow you to edit your data set xold.4 Editing data
When invoked on a data frame or matrix. User-contributed packages can be a rich source of datasets.
.
7. in a few cases it loads several objects. so see the on-line help for the object to see what to expect. However. Use > xnew <. the simplest way is to use fix(xold).

The bandwidth bw was chosen by trial-and-error as the default gives too much smoothing (it usually does for “interesting” densities). bw=0. prob=TRUE) > lines(density(eruptions.6. seq(1.163 4.6000 2. 1st Qu.454 > fivenum(eruptions) [1] 1. make a plot of density > hist(eruptions.1585 4. 1. 5. 5.100
The decimal point is 1 digit(s) to the left of the | 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 | | | | | | | | | | | | | | | | | | 070355555588 000022233333335577777777888822335777888 00002223378800035778 0002335578023578 00228 23 080 7 2337 250077 0000823577 2333335582225577 0000003357788888002233555577778 03335555778800233333555577778 02222335557780000000023333357778888 0000233357700000023578 00000022335800333 0370
A stem-and-leaf plot is like a histogram.4585 5.1000 > stem(eruptions)
Max.Chapter 8: Probability distributions
36
8. The simplest is to examine the numbers.600 2.0000 4.1)) > rug(eruptions) # show the actual data points More elegant density plots can be made by density. (Better automated
. Median Mean 3rd Qu. > attach(faithful) > summary(eruptions) Min. Two slightly different summaries are given by summary and fivenum and a display of the numbers by stem (a “stem and leaf” plot).2).2.000 3. > hist(eruptions) ## make the bins smaller. and R has a function hist to plot histograms.488 4.2 Examining the distribution of a set of data
Given a (univariate) set of data we can examine its distribution in a large number of ways. and we added a line produced by density in this example. 0.

"pnorm".01052 and the Kolmogorov-Smirnov test > ks. df = 5) qqnorm(x).)
. x. p-value = 0.5
4. mean = mean(long). sd = sqrt(var(long))) One-sample Kolmogorov-Smirnov test data: long D = 0. df = 5). p-value = 0. we might want a more formal test of agreement with normality (or not).test(long) Shapiro-Wilk normality test data: long W = 0.4284 alternative hypothesis: two.Chapter 8: Probability distributions
38
par(pty="s") # arrange for a square figure region qqnorm(long).9793.0661. We can make a Q-Q plot against the generating distribution by qqplot(qt(ppoints(250).0
−2
−1
0
1
2
Theoretical Quantiles
x <.5
5.sided (Note that the distribution theory is not valid here as we have estimated the parameters of the normal distribution from the same sample. Let us compare this with some simulated data from a t distribution
Normal Q−Q Plot
Sample Quantiles
3.0
3.test(long. R provides the Shapiro-Wilk test > shapiro.rt(250.0
4. qqline(x) which will usually (if it is a random sample) show longer tails than expected for a normal. qqline(long) which shows a reasonable fit but a shorter right tail than one would expect from a normal distribution. xlab = "Q-Q plot for t dsn") qqline(x) Finally.

03 80.00 80.04 80.03 79.test(A.Chapter 8: Probability distributions
39
8. B) which indicates that the first group tends to give higher results than the second.03 79.00 80.04
1
2
To test for the equality of the means of the two examples.03 80.97 80.94
79. all “classical” tests including the ones used below are in package stats which is normally loaded.scan() 80. A <.03 80.03 80.490) Method A: 79. B) Welch Two Sample t-test data: A and B t = 3.98
80. p-value = 0.98 80. Note that in R.98 79.02 B <.027.02 79.and two-sample tests
So far we have compared a single sample to a normal distribution.00
80.97 Boxplots provide a simple graphical comparison of the two samples.97 80.04 80.05 80.95 79.02 80. A much more common operation is to compare aspects of two samples.04 79.02 80. df = 12.97 79.97 80.02 80. we can use an unpaired t-test by > t.02 79.98 79.94 79.03 80.97 boxplot(A.94 79.04 80.02
80.02 Method B: 80.00694 alternative hypothesis: true difference in means is not equal to 0
.3 One.96
79.
79.97 79. p.05 80.2499. Consider the following sets of data on the latent heat of the fusion of ice (cal/gm) from Rice (1995.02 80.03 80.scan() 79.95 79.04 80.04 79.98 80.97 80.

df = 19.07018320 sample estimates: mean of x mean of y 80.007497 alternative hypothesis: true location shift is not equal to 0
.5837. We can use the F test to test for equality in the variances. assuming normality. The two-sample Wilcoxon (or Mann-Whitney) test only assumes a common continuous distribution under the null hypothesis. B.equal=TRUE) Two Sample t-test data: A and B t = 3.test(A.test function).01385526 0.test(A.02077 79.06734788 sample estimates: mean of x mean of y 80.1052687 sample estimates: ratio of variances 0. p-value = 0. provided that the two samples are from normal populations. B) F test to compare two variances data: A and B F = 0.Chapter 8: Probability distributions
40
95 percent confidence interval: 0. > t. and so we can use the classical t-test that assumes equality of the variances. p-value = 0. > wilcox. num df = 12. var. > var.97875 All these tests assume normality of the two samples.4722. denom df = 7.02077 79.002551 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: 0.01669058 0.test(A. p-value = 0.3938 alternative hypothesis: true ratio of variances is not equal to 1 95 percent confidence interval: 0.97875 which does indicate a significant difference. B) Wilcoxon rank sum test with continuity correction data: A and B W = 89.5837405 which shows no evidence of a significant difference. By default the R function does not assume equality of variances in the two samples (in contrast to the similar S-Plus t.1251097 2.

points=FALSE.points=FALSE.05919 alternative hypothesis: two-sided Warning message: cannot compute correct p-values with ties in: ks.test(A. B) Two-sample Kolmogorov-Smirnov test data: A and B D = 0. There are several ways to compare graphically the two samples. add=TRUE) will show the two empirical CDFs. p-value = 0.test(A. B) Note the warning: there are several ties in each sample. verticals=TRUE. and qqplot will perform a Q-Q plot of the two samples.Chapter 8: Probability distributions
41
Warning message: Cannot compute exact p-value with ties in: wilcox. The following > plot(ecdf(A).5962. do. which suggests strongly that these data are from a discrete distribution (probably due to rounding). assuming a common continuous distribution: > ks.test(A. B)
. B)) > plot(ecdf(B). The Kolmogorov-Smirnov test is of the maximal vertical distance between the two ecdf’s. do. xlim=range(A. verticals=TRUE. We have already seen a pair of boxplots.

1 Grouped expressions
R is an expression language in the sense that its only command type is a function or expression which returns a result.1 Conditional execution: if statements
The language has available a conditional construction of the form > if (expr_1 ) expr_2 else expr_3 where expr 1 must evaluate to a single logical value and the result of the entire expression is then evident. ind) > for (i in 1:length(yc)) { plot(xc[[i]]. One possibility here is to use coplot(). Whereas & and | apply element-wise to vectors. loops and conditional execution
9.2 Repetitive execution: for loops. with elements a[i] if condition[i] is true. suppose ind is a vector of class indicators and we wish to produce separate plots of y versus x within classes.2 Control statements
9. and it may be used wherever any expression may be used. repeat and while
There is also a for loop construction which has the form > for (name in expr_1 ) expr_2 where name is the loop variable. . The “short-circuit” operators && and || are often used as part of the condition in an if statement. {expr_1 . be itself included in parentheses and used a part of an even larger expression. This has the form ifelse(condition. ind) > yc <. in particular multiple assignments are possible. or use xyplot from package lattice. for example.2. otherwise b[i]. the ifelse function. and only evaluate their second argument if necessary.1 which will produce an array of plots corresponding to each level of the factor. expr_m }. a. Since such a group is also an expression it may. There is a vectorized version of the if/else construct. now putting all plots on the one display. Commands may be grouped together in braces. && and || apply to vectors of length one.
9. and so on.split(y. b) and returns a vector of the length of its longest argument.. and expr 2 is often a grouped expression with its sub-expressions written in terms of the dummy name.. in which case the value of the group is the result of the last expression in the group evaluated. As an example.Chapter 9: Grouping. expr 1 is a vector expression.. Another way to do this. Even an assignment is an expression whose result is the value assigned.split(x. loops and conditional execution
42
9 Grouping. is as follows: > xc <. (often a sequence like 1:20).
.2. expr 2 is repeatedly evaluated as name ranges through the values in the vector result of expr 1. yc[[i]])
1
to be discussed later.
9.

Control statements are most often used in connection with functions which are discussed in Chapter 10 [Writing your own functions]. yc[[i]])) } (Note the function split() which produces a list of vectors obtained by splitting a larger vector according to the classes specified by a factor. and where more examples will emerge. page 44. loops and conditional execution
43
abline(lsfit(xc[[i]]. possibly abnormally. mostly used in connection with boxplots. See the help facility for further details. This is the only way to terminate repeat loops. The break statement can be used to terminate any loop. Code that takes a ‘whole object’ view is likely to be both clearer and faster in R. This is a useful function.Chapter 9: Grouping. The next statement can be used to discontinue one particular cycle and skip to the “next”.
. Other looping facilities include the > repeat expr statement and the > while (condition ) expr statement.) Warning: for() loops are used in R code much less often than in compiled languages.

qr(X) qr.function(y1. (This is ordinarily called the least squares estimate of the regression coefficients.yb2)/sqrt(s*(1/n1 + 1/n2)) tst } With this function defined. are themselves written in R and thus do not differ materially from user written functions. X. consider a function to emulate directly the Matlab backslash command. expr_2..) This would ordinarily be done with the qr() function.mean(y2) s1 <.((n1-1)*s1 + (n2-1)*s2)/(n1+n2-2) tst <. data$female).coef(X.1 Simple examples
As a first example. arg i. A function is defined by an assignment of the form > name <. This is an artificial example. consider a function to calculate the two sample t-statistic. that uses the arguments.function(X. to calculate a value. y2) { n1 <. convenience and elegance. y) { X <. where (X T X)− is a generalized inverse of X X. showing “all the steps”. you could perform two sample t-tests using a call such as > tstat <. which returns the coefficients of the orthogonal projection of the vector y onto the column space of the matrix. simpler ways of achieving the same end.. of course.var(y1).) and may occur anywhere a function call is legitimate. n2 <.. s2 <. however this is sometimes a bit tricky to use directly and it pays to have a simple function such as the following to use it safely. These are true R functions that are stored in a special internal form and may be used in further expressions and so on.twosam(data$male.length(y2) yb1 <. The value of the expression is the value returned for the function. yb2 <. In the process.(yb1 .mean(y1). since there are other. the R language allows the user to create objects of mode function.length(y1). var().Chapter 10: Writing your own functions
44
10 Writing your own functions
As we have seen informally along the way. A call to the function then usually takes the form name (expr_1. It should be emphasized that most of the functions supplied as part of the R system. . tstat As a second example. (usually a grouped expression)..function(arg_1. . postscript() and so on. such as mean(). arg_2.
10. The function is defined as follows: > twosam <. > bslash <. and learning to write useful functions is one of the main ways to make your use of R comfortable and productive.var(y2) s <. y)
.) expression The expression is an R expression. Thus given a n by 1 vector y and an n by p matrix X then X y is defined as (X T X)− X T y. the language gains enormously in power.

10. if arguments to called functions are given in the “name =object ” form.function(X.. data.2 Defining new binary operators
Had we given the bslash() function a different name. TRUE.function(data. Furthermore the argument sequence may begin in the unnamed. If so.) The function could then be used as X %!% y.) The matrix multiplication operator.frame. limit=20) > ans <. graph=TRUE. page 8.3 Named arguments and defaults
As first noted in Section 2. page 54
. The classical R function lsfit() does this job quite well. For example. data.fun1(d. y) { ..
10. we choose ! for the internal character. The function definition would then start as > "%!%" <. limit=20. and more1 . and specify named arguments after the positional arguments. Suppose. limit) { [function body omitted] } then the function may be invoked in several ways.frame. graph=TRUE. if fun1 were defined as > fun1 <.bslash(Xmat. they may be given in any order. 20) > ans <.coef() in the slightly counterintuitive way above to do this part of the calculation. namely one of the form %anything % it could have been used as a binary operator in expressions rather than in function form. graph.frame=df) are all equivalent. data. graph=TRUE. df. and the outer product matrix operator %o% are other examples of binary operators defined in this way. } it could be called as
1
See also the methods described in Chapter 11 [Statistical models in R]. %*%. positional form. limit=20) { . in which case they may be omitted altogether from the call when the defaults are appropriate. (The backslash symbol itself is not a convenient choice as it presents special problems in this context. Thus if there is a function fun1 defined by > fun1 <.. In many cases arguments can be given commonly appropriate default values.fun1(data=d. for example > ans <. } (Note the use of quote marks.function(data. yvar) and so on.. we may wish to make it a matrix binary operator for even more convenient use. It in turn uses the functions qr() and qr. for example. Hence there is probably some value in having just this part isolated in a simple to use function if it is going to be in frequent use.fun1(d.3 [Generating regular sequences]. df.Chapter 10: Writing your own functions
45
} After this object is created it may be used in statements such as > regcoeff <.

even involving other arguments to the same function. .has different semantics in R. though hardly difficult.
10. or as > ans <.qr(X) does not affect the value of the argument in the calling program. respectively. page 73.1 Efficiency factors in block designs
As a more complete.) [more omissions] }
10. and
. limit=10) which changes one of the defaults.frame.. If global and permanent assignments are intended within a function. An outline example is given below..) A block design is defined by two factors.) This can be done by including an extra argument.6 More advanced examples
10.fun1(d.1 [The par() function]. <<.. say blocks (b levels) and varieties (v levels).
10.’. graph=TRUE.4. which may then be passed on. fun1 <.Chapter 10: Writing your own functions
46
> ans <. if a little pedestrian. they are not restricted to be constants as in our simple example here. These are discussed further in Section 10. of the function. Thus the assignment X <.) { [omitted statements] if (graph) par(pch="*". df) which is now equivalent to the three cases above... This is a somewhat advanced. page 48. df. consider finding the efficiency factors for a block design. (See Section 12..7 [Scope]. (Some aspects of this problem have already been discussed in Section 5. See the help document for details..fun1(d. . example of a function. data. limit=20. If R and K are the v by v and b by b replications and block size matrices.or the function assign() can be used. then either the “superassignment” operator.function(data.4 The ‘.6. For example many graphics functions use the function par() and functions like plot() allow the user to pass on graphical parameters to par() to control the graphical output. S-Plus users should be aware that <<. It is important to note that defaults may be arbitrary expressions. literally ‘.’ argument
Another frequent requirement is to allow one function to pass on argument settings to another.. To understand completely the rules governing the scope of R assignments the reader needs to be familiar with the notion of an evaluation frame. for more details on the par() function. page 20. topic and is not covered further here.5 Assignments within functions
Note that any ordinary assignments done within the function are local and temporary and are lost after exit from the function.3 [Index matrices].

shown below.vector(table(blocks)) # remove dim attr R <. no. d <. no. where A = K −1/2 N R−1/2 .factor(blocks) # minor safety move b <. since sometimes these give additional useful qualitative information.list(rep("". varieties) { blocks <.length(levels(blocks)) varieties <.X > dimnames(temp) <.as.2 Dropping all names in a printed array
For printing purposes with large matrices or arrays. Removing the dimnames attribute will not achieve this effect. varieties) A <. v)) sv <.sv$d^2. varietycv=sv$v) } It is numerically slightly better to work with the singular value decomposition on this occasion rather than the eigenvalue routines.svd(A) list(eff=1 .factor(varieties) # minor safety move v <.vector(table(varieties)) # remove dim attr N <. rep(b.as.list() l <. but also the block and variety canonical contrasts.rep("". but rather the array must be given a dimnames attribute consisting of empty strings. > bdeff <. then the efficiency factors are defined as the eigenvalues of the matrix E = Iv − R−1/2 N T K −1 N R−1/2 = Iv − AT A.d a } With this function defined.1/sqrt(K) * N * rep(1/sqrt(R). It also illustrates how some effective and useful user functions can be quite short. ncol(X))) > temp.as. it is often useful to print them in close block form without the array names or numbers.dimnames <. i) } dimnames(a) <.
10.as. an array may be printed in close format using
. The result of the function is a list giving not only the efficiency factors as the first component.function(blocks. rep("".Chapter 10: Writing your own functions
47
N is the b by v incidence matrix. nrow(X)). X > temp <.length(levels(varieties)) K <.dimnames().0 for(i in dim(a)) { d[[l <. as a “wrap around” to achieve the same result. blockcv=sv$u.function(a) { ## Remove all dimension names from an array for compact printing.6. One way to write the function is given below.table(blocks.l + 1]] <. For example to print a matrix. rm(temp) This can be much more conveniently done using a function.

lim.0e-06.3 Recursive numerical integration
Functions may be recursive.6. that such functions. eps. it details one of the major differences between S-Plus and R. fun1) }
10. d. eps. lim . a.7 Scope
The discussion in this section is somewhat more technical than in other parts of this document. The integrand is evaluated at the end points of the range and in the middle.(b . However.function(f.function(f.dimnames(X) This is particularly useful for large integer arrays.((fa + fb) * (b . fb. eps = 1.a)/4 fd <. a1. eps.a1 .Chapter 10: Writing your own functions
48
> no. and may themselves define functions within themselves. The example is also given partly as a little puzzle in R programming. fun)) } } fa <.(a + b)/2 h <. fd. b. b. fd. Their values are determined by the process of binding the actual function arguments to the formal parameters.f(d) a1 <.h * (fa + fd) a2 <. The result is an adaptive integration process that concentrates function evaluations in regions where the integrand is farthest from linear.f(b) a0 <. b. a0. fun) { ## function ‘fun1’ is only visible inside ‘area’ d <. local variables and free variables. Otherwise the same process is recursively applied to each panel.a))/2 fun1(f. eps. fb. fb. The example below shows a naive way of performing one-dimensional numerical integration. a2. fa. b. however. where patterns are the real interest rather than the values. fa. are not inherited by called functions in higher evaluation frames as they would be if they were on the search path. a0. a.a2) < eps || lim == 0) return(a1 + a2) else { return(fun(f. a heavy overhead. Note. The formal parameters of a function are those occurring in the argument list of the function. If the one-panel trapezium rule answer is close enough to the two panel. There is. lim. fun) + fun(f. however. The symbols which occur in the body of a function can be divided into three classes. a. lim = 10) { fun1 <.1. or indeed variables. fa.f(a) fb <. then the latter is returned as the value. d. lim . and the function is only competitive with other algorithms when the integrand is both smooth and very difficult to evaluate. a. Local
.h * (fd + fb) if(abs(a0 . formal parameters. area <.
10.1.

Therefore it is a free variable and the scoping rules must be used to ascertain the value that is to be associated with it. Variables which are not formal parameters or local variables are called free variables. First we define a function called cube. cube <. The special assignment operator. The difference between evaluation in R and evaluation in S-Plus is that S-Plus looks for a global variable called n while R first looks for a variable called n in the environment created when cube was invoked. ## first evaluation in S S> cube(2) Error in sq(): Object "n" not found Dumped S> n <. <<-. they will have access to its value. In R the free variable bindings are resolved by first looking in the environment in which the function was created.2*x print(x) print(y) print(z) } In this function. Under lexical scope (R) it is the parameter to the function cube since that is the active binding for the variable n at the time the function sq was defined. Consider the following function definition.function(x) { y <. Free variables become local variables if they are assigned to.function() n*n n*sq() } The variable n in the function sq is not an argument to that function.3 S> cube(2) [1] 18 ## then the same function evaluated in R R> cube(2) [1] 8 Lexical scope can also be used to give functions mutable state. y is a local variable and z is a free variable. a function for making deposits and a function for stating the current balance. f <. We achieve this by creating the three functions within account and then returning a list containing them. In the following example we show how R can be used to mimic a bank account. Under static scope (S-Plus) the value is that associated with a global variable named n.Chapter 10: Writing your own functions
49
variables are those whose values are determined by the evaluation of expressions in the body of the functions. A functioning bank account needs to have a balance or total. When account is invoked it takes a numerical argument total and returns a list containing the three functions. Because these functions are defined in an environment which contains total.function(n) { sq <. This is called lexical scope. This operator looks back in enclosing environments for an environment that contains the
. x is a formal parameter. a function for making withdrawals. is used to change the value associated with total.

total . Your balance is". "\n\n") } ) } ross <.Last can be used. with the value of right hand side.account(100) robert <. withdraw = function(amount) { if(amount > total) stop("You don’t have that much money!\n") total <<. total. If that variable is unset. Finally. the file ‘Rprofile.amount cat(amount. total. total. open.
.site’ in the R home
2
In some sense this mimics the behavior in S-Plus since in S-Plus this operator always creates or assigns to a global variable.open. the special functions . Only when <<.8 Customizing the environment
Users can customize their environment in several different ways. For most users <<. The location of the site initialization file is taken from the value of the R_PROFILE environment variable.First and . "\n\n") }.account <.total + amount cat(amount.has been used in a function that was returned as the value of another function will the special behavior described here occur. Your balance is".account(200) ross$withdraw(30) ross$balance() robert$balance() ross$deposit(50) ross$balance() ross$withdraw(500)
10.function(total) { list( deposit = function(amount) { if(amount <= 0) stop("Deposits must be positive!\n") total <<.open.creates a global variable and assigns the value of the right hand side to it2 . balance = function() { cat("Your balance is".Chapter 10: Writing your own functions
50
symbol total and when it finds such an environment it replaces the value. "withdrawn. There is a site initialization file and every directory can have its own special initialization file. "\n\n") }. If the global or top-level environment is reached without finding the symbol total then that variable is created and assigned to there. in that environment. "deposited.

This file gives individual users control over their workspace and allows for different startup procedures in different working directories. If no ‘.function() { options(prompt="$ ". Thus.
.9 Classes. A second. If R is invoked in that directory then that file will be sourced. summary() for summarizing analyses of various types. profile file named ‘.Last <.function() { graphics. A definition in later files will mask definitions in earlier files.First() in either of the two profile files or in the ‘. An example makes things clearer. > . The class mechanism offers the user the facility of designing and writing generic functions for special purposes. For example."\nAdios\n")) } # a small safety measure. is (normally) executed at the very end of the session. "mystuff. This file should contain the commands that you want to execute every time R is started under your system. the definition in the example below alters the prompt to $ and sets up various other useful things that can then be taken for granted in the rest of the session.Rprofile’ files. Put the other way round. the file it points to is used instead of the ‘. > .Rprofile’3 can be placed in any directory. there is always a default action provided. # Is it time for lunch?
10. continue="+\t") # $ is the prompt options(digits=5.Rprofile’ file in the user’s home directory and uses that (if it exists). ‘Rprofile.RData’ and then . If the argument lacks any class attribute.Last(). if defined. Among the other generic functions are plot() for displaying objects graphically. length=999) # custom numbers and printout x11() # for graphics par(pch = "+") # plotting character source(file.
3
So it is hidden under UNIX.getenv("HOME").path(Sys. then R looks for a ‘. or has a class not catered for specifically by the generic function in question. It is automatically performed at the beginning of an R session and may be used to initialize the environment. the user profile.First <.Rprofile’ file is found in the startup directory. and anova() for comparing statistical models.R")) # my personal functions library(MASS) # attach a package } Similarly a function .First(). a generic function performs a task or action on its arguments specific to the class of the argument itself.off() cat(paste(date(). An example is given below. Any function named . generic functions and object orientation
The class of an object determines how it will be treated by what are known as generic functions. ‘. personal. the sequence in which files are executed is.Chapter 10: Writing your own functions
51
subdirectory ‘etc’ is used.site’.RData’ image has a special status. "R". If the environment variable R_PROFILE_USER is set.

the functions that can accommodate in some fashion objects of class "data.aov’ was found It was found in the following places registered S3 method for coef from namespace stats namespace:stats with value function (object.nls* coef. and more.frame" include [ [<[[<mean any plot as. "factor". for example > coef function (object.frame") Conversely the number of classes a generic function can handle can also be quite large. We can read these by either of > getAnywhere("coef.nls* coef.. .na(z)] } > getS3method("coef".) { z <..Chapter 10: Writing your own functions
52
The number of generic functions that can treat a class in a specific way can be quite large.matrix summary
A currently complete list can be got by using the methods() function: > methods(class="data.default* coef.. "aov") function (object.aov") A single object matching ’coef..frame".object$coef z[!is.Arima* coef. . "density".. For example.object$coef z[!is.aov* [5] coef. A complete list can be got again by using the methods() function: > methods(plot) For many generic functions the function body is quite short.) { z <. . For example the plot() function has a default method and variants for objects of classes "data. To see what methods are available we can use methods() > methods(coef) [1] coef.summary.) UseMethod("coef") The presence of UseMethod indicates this is a generic function.listof*
Non-visible functions are asterisked In this example there are six methods.. none of which can be seen by typing its name.na(z)] }
.

Chapter 10: Writing your own functions
53
The reader is referred to the R Language Definition for a more complete discussion of this mechanism.
.

y~0+x y ~ -1 + x y ~ x . Suppose y.1 Defining statistical models. Later we make some rather more ambitious presumptions. in particular with regression analysis and the analysis of variance. and the second an explicit one. . x2. . log(y). . homoscedastic errors
p
yi =
j=0
βj xij + ei . . . y~x y~1+x Both imply the same simple linear regression model of y on x. . The requirements for fitting statistical models are sufficiently well defined to make it possible to construct general tools that apply in a broad spectrum of problems. As we mention in the introduction. and one needs to ask for the details by calling extractor functions. x1.
ei ∼ NID(0. The first has an implicit intercept term.Chapter 11: Statistical models in R
54
11 Statistical models in R
This section presumes the reader has some familiarity with statistical methodology.
. the basic output is minimal. The following formulae on the left side below specify statistical models as described on the right. . on x1 and x2 (with an implicit intercept term). without an intercept term). X is the model matrix or design matrix and has columns x0 . . the determining variables.
Examples
Before giving a formal specification.
i = 1. B. formulae
The template for a statistical model is a linear regression model with independent. . X is a matrix and A. . are numeric variables. n
In matrix terms this would be written y = Xβ + e where the y is the response vector. xp . . x. . x1 . are factors. x0.
11. . Very often x0 will be a column of ones defining an intercept term. . C.1 Simple linear regression of y on x through the origin (that is. log(y) ~ x1 + x2 Multiple regression of the transformed variable. R provides an interlocking suite of facilities that make fitting statistical models very simple. σ 2 ). a few examples may usefully set the picture. namely that something is known about generalized linear models and nonlinear regression.

determined by factor C.. and error strata determined by factor C.1 Separate simple linear regression models of y on x within the levels of A. or 1.
A*B A + B + A:B B %in% A A/B Two factor non-additive model of y on A and B. for an ordinary linear model. with classes determined by A. The form. A and B. Both formulae specify the same model. is either • a vector or matrix expression. and the second uses explicit powers. (the first is optional). implying the inclusion or exclusion of a term in the model.
. For example a split plot experiment. and with covariate x.2) y ~ 1 + x + I(x^2) Polynomial regression of y on x of degree 2. y ~ A*B + Error(C) An experiment with two treatment factors.2) Multiple regression y with model matrix consisting of the matrix X as well as polynomial terms in x to degree 2. (or expression evaluating to a vector or matrix) defining the response variable(s). The operator ~ is used to define a model formula in R.A:B:C Three factor experiment but with a model containing main effects and two factor interactions only. In abstract terms all four specify the same model subspace. is response ~ op_1 term_1 op_2 term_2 op_3 term_3 . The first two specify the same crossed classification and the second two specify the same nested classification. is an operator. Single classification analysis of covariance model of y. with different codings. y ~ X + poly(x. with whole plots (and hence also subplots). The first form uses orthogonal polynomials.. with classes determined by A.Chapter 11: Statistical models in R
55
y ~ poly(x. where response op i term i is a vector or matrix. y~A y~A+x y y y y ~ ~ ~ ~ Single classification analysis of variance model of y. either + or -. as basis. The last form produces explicit estimates of as many different intercepts and slopes as there are levels in A. y~A*x y ~ A/x y ~ A/(1 + x) .
y ~ (A + B + C)^2 y ~ A*B*C .

) For ordered factors the k − 1 columns are the orthogonal polynomials on 1. the whole behavior can be changed by the options setting for contrasts. the specification of the parameters being implicit.29): Y ~M Y is modeled as M. M_1 * M_2 M_1 + M_2 + M_1 :M_2 . . . The default setting in R is
. if the intercept is omitted in a model that contains a factor term. (Thus the implicit parameterization is to contrast the response at each level with that at the first. but with a different coding. Inside M all operators have their normal arithmetic meaning. M_1 + M_2 Include M 1 and M 2. This is easy if we have continuous variables. Note particularly that the model formulae specify the columns of the model matrix. . and that term appears in the model matrix. omitting the constant term. then the “subclasses” factor. . the first such term is encoded into k columns giving the indicators for all the levels.1 Contrasts
We need at least some idea how the model formulae specify the columns of the model matrix. . If both terms are factors. First. This is not the case in other contexts. What about a k-level factor A? The answer differs for unordered and ordered factors.’ becomes ‘:’ since the period is a valid name character in R. For unordered factors k − 1 columns are generated for the indicators of the second. vectors or matrices connected by formula operators. as each provides one column of the model matrix (and the intercept will provide a column of ones if included in the model). Although the answer is already complicated.
Note that inside the parentheses that usually enclose function arguments all operators have their normal arithmetic meaning.Chapter 11: Statistical models in R
56
• a factor. M_1 / M_2 M_1 + M_2 %in% M_1 . p. The notation is summarized below (based on Chambers & Hastie. M_1 : M_2 The tensor product of M 1 and M 2. One inevitable change is that the operator ‘.1. A 1 stands for an intercept column and is by default included in the model matrix unless explicitly removed. k. 1992. The formula operators are similar in effect to the Wilkinson and Rogers notation used by such programs as Glim and Genstat.M_2 Include M 1 leaving out terms of M 2.
11. M_1 %in% M_2 Similar to M_1 :M_2 . The function I() is an identity function used to allow terms in model formulae to be defined using arithmetic operators. . kth levels of the factor. Second. . or • a formula expression consisting of factors. M ^n I(M ) All terms in M together with “interactions” up to order n Insulate M. . it is not the whole story. for example in specifying nonlinear models. In all cases each term defines a collection of columns either to be added to or removed from the model matrix. M_1 .

So if you need to compare your results to those of a textbook or paper which used S-Plus. and a streamlined version of the call is as follows: > fitted. plotted and so on by using generic functions that orient themselves to objects of class "lm".
11.
.3 Generic functions for extracting model information
The value of lm() is a fitted model object. you will need to set options(contrasts = c("contr. as the contrast scheme to be used can be set for each term in the model using the functions contrasts and C. Although the details are complicated. data = production) would fit a multiple regression model of y on x1 and x2 (with implicit intercept term). S using Helmert contrasts. provided that marginality is preserved. "contr. Fitting.poly")) This is a deliberate difference.2 Linear models
The basic function for fitting ordinary multiple models is lm().
11. This is the case regardless of whether data frame production has been attached on the search path or not. and is for experts only. The important (but technically optional) parameter data = production specifies that any variables needed to construct the model should come first from the production data frame. Information about the fitted model can then be displayed. We have not yet considered interaction terms: these generate the products of the columns introduced for their component terms.helmert".Chapter 11: Statistical models in R
57
options(contrasts = c("contr. Long form: coefficients(object ).lm(y ~ x1 + x2. as treatment contrasts (R’s default) are thought easier for newcomers to interpret. model formulae in R will normally generate the models that an expert statistician would expect. coef(object ) Extract the regression coefficient (matrix).frame ) For example > fm2 <.model <. These include add1 deviance formula predict step alias drop1 kappa print summary anova effects labels proj vcov coef family plot residuals A brief description of the most commonly used ones is given below. data = data.lm(formula.treatment". extracted. object_2 ) Compare a submodel with an outer model and produce an analysis of variance table. "contr. technically a list of results of class "lm". a model with an interaction but not the corresponding main effects will in general lead to surprising results.poly")) The main reason for mentioning this is that R and S have different defaults for unordered factors. We have still not finished. anova(object_1. for example.

formula(object ) Extract the model formula. It should be noted that in addition aov() allows an analysis of models with multiple error strata such as split plot experiments. weighted if appropriate.data) would typically be used to describe an experiment with mean model v + n*p*k and three error strata.frame ) The data frame supplied must have variables specified with the same labels as the original. The model with the smallest value of AIC (Akaike’s An Information Criterion) discovered in the stepwise search is returned.frame ) operates at the simplest level in a very similar way to the function lm(). data=farm. namely “between farms”. newdata=data. The model formula response ~ mean. The value is a vector or matrix of predicted values corresponding to the determining variable values in data.Chapter 11: Statistical models in R
58
deviance(object ) Residual sum of squares.formula ) specifies a multi-stratum experiment with error strata defined by the strata.formula + Error(strata.4 Analysis of variance and model comparison
The model fitting function aov(formula. or balanced incomplete block designs with recovery of inter-block information. fitted values and some diagnostics. showing residuals.formula. predict(object. plot(object ) Produce four plots.frame. “within farms. For example. and most of the generic functions listed in the table in Section 11. Short form: resid(object ). vcov(object ) Returns the variance-covariance matrix of the main parameters of a fitted model object. when it defines a two strata experiment. In the simplest case. with all determining variables factors. namely between and within the levels of the factor.formula is simply a factor. summary(object ) Print a comprehensive summary of the results of the regression analysis.3 [Generic functions for extracting model information]. data=data.
. weighted as appropriate. Most often used implicitly.
11. strata. page 57 apply. print(object ) Print a concise version of the object. residuals(object ) Extract the (matrix of) residuals. step(object ) Select a suitable model by adding or dropping terms and preserving hierarchies.aov(yield ~ v + n*p*k + Error(farms/blocks). a model formula such as that in: > fm <. between blocks” and “within blocks”.

new. For further details. For multistratum experiments the procedure is first to project the response onto the error strata. fit an additional model including a sixth regressor variable. This does not give different information to the default. The name ‘.4.2.) The display is then an ANOVA table showing the differences between the fitted models when fitted in sequence.update(fm05. + x6) > smf6 <. Hence only for orthogonal experiments will the order of inclusion be inconsequential.) ~ .formula the special name consisting of a period.
. The sums of squares shown are the decrease in the residual sums of squares resulting from an inclusion of that term in the model at that place in the sequence. The names of these give a good clue to their purpose.model. drop1() and step().) would fit a five variate multiple regression with variables (presumably) from the data frame production. ~ .1. A more flexible alternative to the default full ANOVA table is to compare two or more models directly using the anova() function. but rather makes it easier to comprehend and control.model <. The fitted models being compared would usually be an hierarchical sequence. but with slightly different meaning. > anova(fitted. data = production) would fit a model with response y and regressor variables all other variables in the data frame production.formula ) In the new.Chapter 11: Statistical models in R
59
11. can be used to stand for “the corresponding part of the old model formula”.model. ‘. see Chambers & Hastie (1992). sqrt(. Note especially that if the data= argument is specified on the original call to the model fitting function. again in sequence.lm(y ~ x1 + x2 + x3 + x4 + x5.. For example > fmfull <. of course.’ can also be used in other contexts.
11. Its form is > new. and to fit the mean model to each projection. data = production) > fm6 <. . fitted. Other functions for exploring incremental sequences of models are add1().lm(y ~ .update(fm6.5 Updating fitted models
The update() function is largely a convenience function that allows a model to be fitted that differs from one previously fitted usually by just a few additional or removed terms.update(old.. this information is passed on through the fitted model object to update() and its allies.model. > fm05 <. For example.’. but for full details see the on-line help. and fit a variant on the model where the response had a square root transform applied. . . only.1 ANOVA tables
Note also that the analysis of variance table (or tables) are for a sequence of fitted models.

. So it is assumed that the distribution of y is determined by its mean and possibly a scale parameter as well. binomial. • The distribution of y is of the form fY (y. A generalized linear model may be described in terms of the following sequence of assumptions: • There is a response. . Those automatically available are shown in the following table: Family name binomial gaussian Gamma Link functions logit. ϕ) = exp A {yλ(µ) − γ (λ(µ))} + τ (y.6. • The mean. µ. poisson. inverse identity. is a smooth invertible function of the linear predictor: µ = m(η). These assumptions are loose enough to encompass a wide class of models useful in statistical practice. A represents a prior weight. inverse.6 Generalized linear models
Generalized linear modeling is a development of linear models to accommodate both nonnormal response distributions and transformations to linearity in a clean and straightforward way. Each response distribution admits a variety of link functions to connect the mean with the linear predictor. ϕ) ϕ
where ϕ is a scale parameter (possibly known). but tight enough to allow the development of a unified methodology of estimation and inference. log
. . This linear function is called the linear predictor. only. x2 . such as McCullagh & Nelder (1989) or Dobson (1990). In the latter case the variance function must be specified as a function of the mean. but in other cases this function is implied by the response distribution. µ. and is constant for all observations.Chapter 11: Statistical models in R
60
11. whose values influence the distribution of the response.
11. log. y. hence xi has no influence on the distribution of y if and only if βi = 0. log.1 Families
The class of generalized linear models handled by facilities supplied in R includes gaussian. cloglog identity. at least approximately. is called the link function. The reader is referred to any of the current reference works on the subject for full details. and is usually written η = β1 x1 + β2 x2 + · · · + βp xp . • The stimulus variables influence the distribution of y through a single linear function. η = m−1 (µ) = (µ)
and this inverse function. assumed known but possibly varying with the observations. probit. inverse gaussian and gamma response distributions and also quasilikelihood models where the response distribution is not explicitly specified. of interest and stimulus variables x1 . (). and µ is the mean of y. .

Note how the gaussian family is not automatically provided with a choice of links. the name of the link may also be supplied with the family name.model <.Chapter 11: Statistical models in R
61
1/mu^2. so no parameter is allowed. 1/mu^2. log. data=data. data=sales) but much less efficiently. In the case of the quasi family. The data is shown below: Age: 20 35 45 55 70 No.1 [Families]. The family has to be specified in a different way. which is the instrument by which the family is described. probit. log.
The binomial family
Consider a small. family = gaussian. On the Aegean island of Kalythos the male inhabitants suffer from a congenital eye disease. Where there is a choice of links. family=family. sqrt The combination of a response distribution.generator. as we shall see later. tested: 50 50 50 50 50 No. the same mechanism as was used for linear models can still be used to specify the linear part of a generalized model. The R function to fit a generalized linear model is glm() which uses the form > fitted.
The gaussian family
A call such as > fm <. in parentheses as a parameter.glm(y ~ x1 + x2. the effects of which become more marked with increasing age. inverse. supplied family generators are given under “Family Name” in the table in Section 11.gaussian poisson quasi
11.frame ) The only new feature is the family.lm(y ~ x1+x2. cloglog.glm(formula. its use is quite simple.
inverse. It is the name of a function that generates a list of functions and expressions that together define and control the model and estimation process. the variance function may also be specified in this way. a link function and various other pieces of information that are needed to carry out the modeling exercise is called the family of the generalized linear model. Samples of islander males of various ages were tested for blindness and the results recorded. Some examples make the process clear.6. If a problem requires a gaussian family with a nonstandard link. identity. this can usually be achieved through the quasi family.generator. log identity. identity. Although this may seem a little complicated at first sight. page 60.2 The glm() function
Since the distribution of the response depends on the stimulus variables through a single linear function only. artificial example. from Silvey (1970).6. sqrt logit. inverse. data = sales) achieves the same result as > fm <. The names of the standard. blind: 6 17 26 37 44
.

In both cases the LD50 is LD50 = −β0 /β1 that is.data.ld50(coef(fmp)). This is a large and important subject we will not discuss further here. family = binomial.function(b) -b[1]/b[2] > ldp <. c(ldp.kalythos$y) To fit the models we use > fmp <. and to estimate for each model the LD50. F (z) = ez /(1 + ez ). y = c(6. ldl <. F (z) = Φ(z) is the standard normal distribution function. both models have the form y ∼ B(n. kalythos$n .frame(x = c(20.cbind(kalythos$y. and so must be a 0/1 vector. It even forms a major part of the use of non-gaussian generalized models overall. If y is the number of blind at age x and n the number tested. its first level is taken as failure (0) and all other levels as ‘success’ (1). To find the LD50 estimate we can use a simple function: > ld50 <.Chapter 11: Statistical models in R
62
The problem we consider is to fit both logistic and probit models to this data.601 years respectively.
Poisson models
With the Poisson family the default link is the log.glm(Ymat ~ x. Here we need the second of these conventions. and in the logit case (the default). • If the response is a factor. so we add a matrix to our data frame: > kalythos$Ymat <. data = kalythos) Since the logit link is the default the parameter may be omitted on the second call. that is the age at which the chance of blindness for a male inhabitant is 50%.45.35.5).17. F (β0 + β1 x)) where for the probit case.44)) To fit a binomial model using glm() there are three possibilities for the response: • If the response is a vector it is assumed to hold binary data.ld50(coef(fml)). the point at which the argument of the distribution function is zero.70). whose actual distribution is often multinomial. To see the results of each fit we could use > summary(fmp) > summary(fml) Both models fit (all too) well. • If the response is a two-column matrix it is assumed that the first column holds the number of successes for the trial and the second holds the number of failures.
. family = binomial(link=probit). The first step is to set the data up as a data frame > kalythos <.37. ldl) The actual estimates from this data are 43. data = kalythos) > fml <.663 years and 43. n = rep(50.26. and in practice the major use of this family is to fit surrogate Poisson log-linear models to frequency data.glm(Ymat ~ x.55.

2. as needed. which provide the functionality (and more) of S-Plus’s ms() and nlminb(). family = quasi(link=inverse.counts)
Quasi-likelihood models
For all families the variance of the response will depend on the mean and will have the scale parameter as a multiplier.1. incidentally. All the methods require initial guesses about what parameter values to try. x2 = −1/z1 . a Poisson generalized linear model may be fitted as in the following example: > fmod <. family = poisson(link=sqrt). data = worm. for example for the poisson distribution Var[y] = µ. variance=constant).7 Nonlinear least squares and maximum likelihood models
Certain forms of nonlinear model can be fitted by Generalized Linear Models (glm()). The form of dependence of the variance on the mean is a characteristic of the response distribution. data = biochem) The reader is referred to the manual and the help document for further information. consider fitting the non-linear regression y= which may be written alternatively as y= 1 +e β1 x1 + β2 x2 θ1 z1 +e z2 − θ2
where x1 = z2 /z1 .glm(y ~ x1 + x2 . For example. nlm() and (from R 2. but rather only a link function and the form of the variance function as it depends on the mean. For quasi-likelihood estimation and inference the precise response distribution is not specified.glm(y ~ A + B + x.
. But in the majority of cases we have to approach the nonlinear curve fitting problem as one of nonlinear optimization. Unlike linear regression for example. and they do this by trying out various parameter values iteratively. As a graceful alternative to the latter.Chapter 11: Statistical models in R
63
Occasionally genuinely Poisson data arises in practice and in the past it was often analyzed as gaussian data after either a log or a square-root transformation. Since quasi-likelihood estimation uses formally identical techniques to those for the gaussian distribution. Supposing a suitable data frame to be set up we could fit this non-linear regression as > nlfit <.
11. R’s nonlinear optimization routines are optim(). this family provides a way of fitting gaussian models with non-standard link functions or variance functions. there is no guarantee that the procedure will converge on satisfactory estimates. β1 = 1/θ1 and β2 = θ2 /θ1 . We seek the parameter values that minimize some index of lack-of-fit.0) nlminb(). and convergence may depend critically upon the quality of the starting values.

06. 1.1). 107. 97. out$minimum is the SSE. but these starting values of 200 and 0.02. 139. K) data: df Vm K 212. and out$estimate are the least squares estimates of the parameters. This method makes sense if the observed errors could have plausibly arisen from a normal distribution.06412123 residual sum-of-squares: 1195. y=y) > fit <. 0. A 95% confidence interval would be the parameter estimate ± 1. 159.data. One way to find sensible starting values is to plot the data.1 + xfit) lines(spline(xfit. 0.7. 200) The fit criterion to be minimized is: > fn <. .11. To obtain the approximate standard errors (SE) of the estimates we do: > sqrt(diag(2*out$minimum/(length(y) .02. 0.06412146 + xfit) lines(spline(xfit. guess some parameter values.212. 191. Vm.1.(p[1] * x)/(p[2] + x))^2) In order to do the fit we need initial estimates of the parameters.seq(. 0.200 * xfit/(0. page 51.22.nls(y ~ SSmicmen(x.seq(.2) * solve(out$hessian))) The 2 in the line above represents the number of parameters. 1.Chapter 11: Statistical models in R
64
11.22. so we can use > df <.02. 0. 152.02.1 Least squares
One way to fit a nonlinear model is by minimizing the sum of the squared errors (SSE) or residuals.1.96 SE.c(0. Now do the
The standard package stats provides much more extensive facilities for fitting non-linear models by least squares.function(p) sum((y . hessian = TRUE) After the fitting. The data are: > x <. df) > fit Nonlinear regression model model: y ~ SSmicmen(x.68384222 * xfit/(0.05) yfit <.449
. 0. 0. We can superimpose the least squares fit on a new plot: > > > > plot(x. y) xfit <.11.1 seem adequate. yfit)) plot(x. yfit))
We could do better. 207. 0. and superimpose the model curve using those values. 201.nlm(fn. .68370711 0. 47.05) yfit <. > > > > fit: > out <.10.56.10) > y <. 0. The model we have just fitted is the Michaelis-Menten model. Here is an example from Bates & Watts (1988). 1. 0.56.frame(x=x.06. y) xfit <. p = c(200. Vm. 123. K).c(76. 1.

and out$estimate are the maximum likelihood estimates of the parameters.Chapter 11: Statistical models in R
65
> summary(fit) Formula: y ~ SSmicmen(x. pp.57e-05 Residual standard error: 10. y)) )) We pick sensible starting values and do the fit: > out <. Error t value Pr(>|t|) Vm 2. hessian = TRUE) After the fitting. 60.412e-02 8. 56. 1. These functions make heavy use of formulae to specify the models. > n <.20). 1.127e+02 6. The recommended nlme package provides functions lme() and nlme() for linear and non-linear mixed-effects models.7842.7242.743 1.281e-03 7.947e+00 30. 59. 108–111. • Mixed models. Vm. 1. The method finds the parameter values which maximize the log likelihood.function(p) sum( . 60) 63. that is linear and non-linear regressions in which some of the coefficients correspond to random effects.93 on 10 degrees of freedom Correlation of Parameter Estimates: Vm K 0. 1.96 SE. The data are: > x <.7651
11.8610. 61.615 3. 18.24e-11 K 6. 1.c(1.8113. K) Parameters: Estimate Std. 53. Here is an example from Dobson (1990). p = c(-50.c(59.7.8839) 52. This example fits a logistic model to dose-response data.7552.6907.c( 6.8 Some non-standard models
We conclude this chapter with just a brief mention of some of the other facilities available in R for special regression and data analysis problems. 1. or equivalently which minimize the negative log-likelihood.(y*(p[1]+p[2]*x) . 1. which clearly could also be fit by glm(). To obtain the approximate SEs of the estimates we do: > sqrt(diag(solve(out$hessian))) A 95% confidence interval would be the parameter estimate ± 1. 62.nlm(fn. > y <.2 Maximum likelihood
Maximum likelihood is a method of nonlinear model fitting that applies even if the errors are not normal.
. 60)
The negative log-likelihood to minimize is: > fn <. 28.
11. 13.n*log(1+exp(p[1]+p[2]*x)) + log(choose(n. out$minimum is the negative log-likelihood.8369. 62.

but many other generic functions such as plot() and text() are well adapted to displaying the results of a tree-based model fit in a graphical way. Function loess is in the standard package stats. recursively. usually one for each determining variable.
. tree-based models seek to bifurcate the data.Chapter 11: Statistical models in R
66
• Local approximating regressions. There are several functions available for fitting regression models in a way resistant to the influence of extreme outliers in the data. Models are again specified in the ordinary linear model form. Tree models are available in R via the user-contributed packages rpart and tree. and as heterogeneous as possible between. implemented in user-contributed packages gam and mgcv. • Additive models. Rather than seek an explicit global linear model for prediction or interpretation. together with code for projection pursuit regression. Such regressions are useful for highlighting a trend in messy data or for data reduction to give some insight into a large data set. The loess() function fits a nonparametric regression by using a locally weighted regression. for example function rlm in package MASS. The model fitting function is tree(). at critical points of the determining variables in order to partition the data ultimately into groups that are as homogeneous as possible within. • Tree-based models. This technique aims to construct a regression function from smooth additive functions of the determining variables. Function lqs in the recommended package MASS provides state-of-art algorithms for highly-resistant fits. • Robust regression. Functions avas and ace in package acepack and functions bruto and mars in package mda provide some examples of these techniques in user-contributed packages to R. An extension is Generalized Additive Models. Less resistant but statistically more efficient methods are available in packages. The results often lead to insights that other data analysis methods tend not to yield.

it is useful to know that the command used is X11() under UNIX. it produces a plot of the values in the vector against their index in the vector. an existing plot. or extract information from.
12. There is a recommended package lattice which builds on grid and provides ways to produce multi-panel plots akin to those in the Trellis system in S. titles and so on. plot(x. using a pointing device such as a mouse.1 High-level plotting commands
High-level plotting functions are designed to generate a complete plot of the data passed as arguments to the function.Chapter 12: Graphical procedures
67
12 Graphical procedures
Graphical facilities are an important and extremely versatile component of the R environment. but in most cases. windows() under Windows and quartz() under Mac OS X. lines and labels. In addition. erasing the current plot if necessary. This manual only describes what are known as ‘base’ graphics. A separate graphics sub-system in package grid coexists with base – it is more powerful but harder to use. Plotting commands are divided into three basic groups: • High-level plotting functions create a new plot on the graphics device. Where appropriate. R maintains a list of graphical parameters which can be manipulated to customize your plots. The graphics facilities can be used in both interactive and batch modes. This is a generic function: the type of plot produced is dependent on the type or class of the first argument. Although this is done automatically. this produces a time-series plot. R plotting commands can be used to produce a variety of graphical displays and to create entirely new kinds of display. It is possible to use the facilities to display a wide variety of statistical graphs and also to build entirely new types of graph.1. labels.1 The plot() function
One of the most frequently used plotting functions in R is the plot() function. If x is a numeric vector. axes. • Low-level plotting functions add more information to an existing plot. labels and titles are automatically generated (unless you request otherwise. If
. interactive use is more productive.) High-level plotting commands always start a new plot. y ) plot(xy ) If x and y are vectors. such as extra points. plot(x. Interactive use is also easy because at startup time R initiates a graphics device driver which opens a special graphics window for the display of interactive graphics. The same effect can be produced by supplying one argument (second form) as either a list containing two elements x and y or a two-column matrix. plot(x ) If x is a time series. y ) produces a scatterplot of y against x. Once the device driver is running. • Interactive graphics functions allow you interactively add information to. possibly with axes.
12.

1. then the command > coplot(a ~ b | c) produces a number of scatterplots of a against b for given values of c.3 Display graphics
Other high-level graphics functions produce different types of plots. it produces a plot of imaginary versus real parts of the vector elements. When three or four variables are involved a coplot may be more enlightening. it is divided into a number of conditioning intervals and for each interval a is plotted against b for values of c within the interval. y is a numeric vector. plot(df ) plot(~ expr ) plot(y ~ expr ) df is a data frame. If c is a factor. that is.
12. You can also use two given variables with a command like > coplot(a ~ b | c + d) which produces scatterplots of a against b for every joint conditioning interval of c and d..
12.g.Chapter 12: Graphical procedures
68
x is a complex vector. this simply means that a is plotted against b for every level of c.smooth(). If X is a numeric matrix or data frame. Some examples are:
. every column of X is plotted against every other column of X and the resulting n(n − 1) plots are arranged in a matrix with plot scales constant over the rows and columns of the matrix.2 Displaying multivariate data
R provides two very useful functions for representing multivariate data. When c is numeric. An example panel function useful for coplots is panel. The third form plots y against every object named in expr. the second form produces boxplots of y for each level of f. The coplot() and pairs() function both take an argument panel= which can be used to customize the type of plot which appears in each panel. expr is a list of object names separated by ‘+’ (e. If a and b are numeric vectors and c is a numeric vector or factor object (all of the same length). y is any object. The number and position of intervals can be controlled with given. The default is points() to produce a scatterplot but by supplying some other low-level graphics function of two vectors x and y as the value of panel= you can produce any type of plot you wish. a + b + c).intervals() is useful for selecting intervals. plot(f ) plot(f. The first two forms produce distributional plots of the variables in a data frame (first form) or of a number of named objects (second form). The first form generates a bar plot of f . y ) f is a factor object.1. the command > pairs(X) produces a pairwise scatterplot matrix of the variables defined by the columns of X.values= argument to coplot()—the function co.

) Constructs a dotchart of the data in x. z. The type= argument controls the type of plot produced. The default.) Plots of three variables. image(x. . nclass=n ) hist(x.4 Arguments to high-level plotting functions
There are a number of arguments which may be passed to high-level graphics functions. . y or both axes to be logarithmic...
Causes the x. This will work for many. . means include axes. hist(x) hist(x.1. log="x" log="y" log="xy" type= Forces the function to act as a low-level graphics function. For example it allows easy visual selection of all data entries with values lying in specified ranges.) Produces a histogram of the numeric vector x.Chapter 12: Graphical procedures
69
qqnorm(x) qqline(x) qqplot(x.. y. The first form plots the numeric vector x against the expected Normal order scores (a normal scores plot) and the second adds a straight line to such a plot by drawing a line through the distribution and data quartiles. the bars represent relative frequencies divided by bin width instead of counts. z. dotchart(x. the contour plot draws contour lines to represent the value of z. and the persp plot draws a 3D surface. . . The image plot draws a grid of rectangles using different colours to represent the value of z. but a recommendation can be given with the nclass= argument. superimposing the plot on the current plot (some functions only)...) contour(x. the breakpoints can be specified exactly with the breaks= argument. y) Distribution-comparison plots. axes=TRUE. as follows: add=TRUE axes=FALSE Suppresses generation of axes—useful for adding your own custom axes with the axis() function. z. Alternatively.) persp(x. as follows: type="p" type="l" type="b" Plot individual points (the default) Plot lines Plot points connected by lines (both)
.. The third form plots the quantiles of x against those of y to compare their respective distributions.. y.
12. In a dotchart the y-axis gives a labelling of the data in x and the x-axis gives its value. A sensible number of classes is usually chosen.. If the probability=TRUE argument is given. but not all... y. breaks=b. types of plot.

labels. and
. The default is 1:length(x). type="n"). low-level plotting commands can be used to add extra information (such as points. Some of the more useful low-level plotting functions are: points(x. Ideal for creating plots with subsequent low-level graphics functions. placed at the top of the plot in a large font. However axes are still drawn (by default) and the coordinate system is set up according to the data.) text(x. sub=string Sub-title.) Add text to a plot at points given by x. in the second.
xlab=string ylab=string Axis labels for the x and y axes. y) lines(x. plot()’s type= argument can also be passed to these functions (and defaults to "p" for points() and "l" for lines().obj ) Adds a line of slope b and intercept a to the current plot. the bottom.. y[i]). h=y may be used to specify y-coordinates for the heights of horizontal lines to go across a plot. lines or text) to the current plot.. y. y. No plotting at all. Use these arguments to change the default labels. . usually the names of the objects used in the call to the high-level plotting function. text(x. names) The graphics parameter type="n" suppresses the points but sets up the axes. y) Adds points or connected lines to the current plot. y. the top of the vertical defines the point. as specified by the character vector names for the points. In the first form.Chapter 12: Graphical procedures
70
type="o" type="h" type="s" type="S" type="n"
Plot points overlaid by lines Plot vertical lines from points to the zero axis (high-density) Step-function plots.
12. Normally labels is an integer or character vector in which case labels[i] is plotted at point (x[i]. Note: This function is often used in the sequence > plot(x. main=string Figure title.2 Low-level plotting commands
Sometimes the high-level plotting functions don’t produce exactly the kind of plot you desire. In this case. b) abline(h=y ) abline(v=x ) abline(lm. placed just below the x-axis in a smaller font. and the text() function supplies special characters. y. abline(a.

For example. Where x and y arguments are required.. or title. This can be achieved in R by specifying an expression rather than a character string in any one of text. Low-level plotting functions usually require some positioning information (e. or fill it if the graphics device allows the filling of figures. counting clockwise from the bottom. polygon(x. mtext. y) and (optionally) shade it in with hatch lines. .. .g.) Draws a polygon defined by the ordered vertices in (x.1 Mathematical annotation
In some cases. it is also sufficient to supply a single argument being a list with elements named x and y. fill=v ) Colors for filled boxes legend( . col=v ) Colors in which points or lines will be drawn legend( .) Adds a legend to the current plot at the specified position. the following code draws the formula for the Binomial probability function:
.2.) Adds an axis to the current plot on the side given by the first argument (1 to 4.
12.. legend. Also lm. in that order. line styles. axis. Similarly a matrix with two columns is also valid input.. axis(side.obj may be list with a coefficients component of length 2 (such as the result of model-fitting functions... Coordinates are given in terms of user coordinates which are defined by the previous high-level graphics command and are chosen based on the supplied data. sub) Adds a title main to the top of the current plot in a large font and (optionally) a sub-title sub at the bottom in a smaller font. Useful for adding custom axes after calling plot() with the axes=FALSE argument. lwd=v ) Line widths legend( .. lty=v ) Line styles legend( . In this way functions such as locator() (see below) may be used to specify positions on a plot interactively.. as follows: legend( . .) Other arguments control the positioning of the axis within or beside the plot. pch=v ) Plotting characters (character vector) title(main.Chapter 12: Graphical procedures
71
v=x similarly for the x-coordinates for vertical lines. and tick positions and labels. are identified with the labels in the character vector legend. y. legend(x. Plotting characters. x and y coordinates) to determine where to place the new plot elements. y. At least one other argument v (a vector the same length as legend) with the corresponding values of the plotting unit must also be given.) which are taken as an intercept and slope. colors etc. it is useful to add mathematical symbols and formulae to a plot.

2 Hershey vector fonts
It is possible to specify Hershey vector fonts for rendering text when using the text and contour functions. locator() is usually called with no arguments. There are three reasons for using the Hershey fonts: • Hershey fonts can produce better output. The type argument allows for plotting at the selected points and has the same effect as for high-level graphics commands.3 Interacting with graphics
R also provides functions which allow users to extract or add information to a plot using a mouse. locator() returns the locations of the points selected as a list with two components x and y. In particular.Chapter 12: Graphical procedures
72
> text(x. It is particularly useful for interactively selecting positions for graphic elements such as legends or labels when it is difficult to calculate in advance where the graphic should be placed. the default is no plotting. to place some informative text near an outlying point. cartographic symbols and astronomical symbols. the command > text(locator(1).) identify(x. atop(n. expression(paste(bgroup("(". The simplest of these is the locator() function: locator(n. (locator() will be ignored if the current device. More information. such as postscript does not support interactive pointing. p^x. including a full listing of the features available can obtained from within R using the commands: > help(plotmath) > example(plotmath) > demo(plotmath)
12. adj=0) may be useful. • Hershey fonts provide cyrillic and japanese (Kana and Kanji) characters. This continues until n (default 512) points have been selected. including tables of Hershey characters can be obtained from within R using the commands: > > > > help(Hershey) demo(Hershey) help(Japanese) demo(Japanese)
12. y. labels) Allow the user to highlight any of the points defined by x and y (using the left mouse button) by plotting the corresponding component of labels nearby (or
. x). or another mouse button is pressed. For example. especially on a computer screen. • Hershey fonts provide certain symbols that may not be available in the standard fonts. ")").2. there are zodiac signs. q^{n-x}))) More information. for rotated and/or small text. "Outlier". y. type) Waits for the user to select locations on the current plot using the left mouse button.

When the process is terminated (see above). "lty")) With a character vector argument. lty=2) With named arguments (or a single list argument). If there is a point near the mouse pointer it will be marked with its index number (that is. R maintains a list of a large number of graphics parameters which control things such as line style. affecting only a single graphics function call. or temporarily. as a list. par() Without arguments.) and a value (a color number. Returns the indices of the selected points when another button is pressed. affecting all graphics functions which access the current device. we may wish the user to select some observation of interest from a graphical display and then manipulate that observation in some way. which controls colors. sets the values of the named graphics parameters. rather than their positions. particularly for presentation or publication purposes.
12. Setting graphics parameters with the par() function changes the value of the parameters permanently. and returns the original values of the parameters as a list. Every graphics parameter has a name (such as ‘col’. For example. Given a number of (x. colors. returns only the named graphics parameters (again. or disable marking altogether with the plot = FALSE argument. however. Graphics parameters can be set in two ways: either permanently.4 Using graphics parameters
When creating graphics. its position in the x/y vectors) plotted nearby. you could use some informative string (such as a case name) as a highlight by using the labels argument to identify(). You can.
par(c("col".1 Permanent changes: The par() function
The par() function is used to access and modify the list of graphics parameters for the current graphics device. customize almost every aspect of the display using graphics parameters. R’s defaults do not always produce exactly that which is required. and each device has a default set of parameters when initialized. y) coordinates in two numeric vectors x and y. Sometimes we want to identify particular points on a plot. y) > identify(x.Chapter 12: Graphical procedures
73
the index number of the point if labels is absent).) par(col=4. we could use the identify() function as follows: > plot(x. for example. figure arrangement and text justification among many others. identify() returns the indices of the selected points. Alternatively. you can use these indices to extract the selected points from the original vectors x and y.4. but simply allows the user to move the mouse pointer and click the left mouse button near a point. returns a list of all graphics parameters and their values for the current device. y) The identify() functions performs no plotting itself. in the sense that all future calls to graphics functions (on the current device)
.) A separate list of graphics parameters is maintained for each active device.
12.

12. You can think of setting graphics parameters in this way as setting “default” values for the parameters.2 Temporary changes: Arguments to graphics functions
Graphics parameters may also be passed to (almost) any graphics function as named arguments. Unfortunately. . You can restore the initial values by saving the result of par() when making changes. plotting commands . which will be used by all graphics functions unless an alternative value is given. This has the same effect as passing the arguments to the par() function. Graphics parameters will be presented in the following form: name =value A description of the parameter’s effect. . . For example: > plot(x.5 Graphics parameters list
The following sections detail many of the commonly-used graphical parameters. . plotting commands . > oldpar <. name is the name of the parameter. that is. and restoring the initial values when plotting is complete.
. without changing the default plotting character for future plots. . pch="+") produces a scatterplot using a plus sign as the plotting character. The R help documentation for the par() function provides a more concise summary. and then restore the original values so as not to affect the user’s R session. this is provided as a somewhat more detailed alternative. except that the changes only last for the duration of the function call. y. This is often undesirable behavior—usually we want to set some graphics parameters. even when par() is called from within a function. the argument name to use in calls to par() or a graphics function. lty=2) .
1
Some graphics parameters such as the size of the current device are for information only.par(col=4. . . > par(oldpar)
12.readonly=TRUE) .Chapter 12: Graphical procedures
74
will be affected by the new value. Note that axes is not a graphics parameter but an argument to a few plot methods: see xaxt and yaxt. Note that calls to par() always affect the global values of graphics parameters. do some plotting.4. value is a typical value you might use when setting the parameter. this is not implemented entirely consistently and it is sometimes necessary to set and reset graphics parameters using par(). > par(oldpar) To save and restore all settable1 graphical parameters use > oldpar <. .par(no.

but it is usually ‘◦’. adj=-0.1 Justification of text relative to the plotting position.sub font=2
The color to be used for axis annotation. 4 to bold italic and 5 to a symbol font (which include Greek letters). When pch is given as an integer between 0 and 25 inclusive. Colors to be used for points.main col. respectively. or some combination of both. 0 means left justify. Not all devices support this." as the plotting character. filled regions and images. use the command > legend(locator(1). In addition.axis col. Alternative line styles are not supported on all graphics devices (and vary on those that do) but line type 1 is always a solid line. Desired width of lines. so a value of -0. and line types 2 and onwards are dotted or dashed lines.1 Graphical elements
R plots are made up of points.lab font. and some have restrictions on the widths that can be used.1 leaves a gap of 10% of the text width between the text and the plotting position. The actual value is the proportion of text that appears to the left of the plotting position. 1 means right justify and 0. To see what the symbols are. 2 to bold face.
. lines. Plotted points tend to appear slightly above or below the appropriate position unless you use ". as follows: pch="+" Character to be used for plotting points. main and sub-titles. as.lab col.5 means to center horizontally about the plotting position. in multiples of the “standard” line width.
pch=4
lwd=2
col=2 col. device drivers arrange so that 1 corresponds to plain text. line type 0 is always invisible.main font.) Graphical parameters exist which control how these graphical elements are drawn. which produces centered points.axis font.character(0:25). Line widths. text. a specialized plotting symbol is produced. but can be coloured in different ways: see the help on points and its examples. main and sub-titles. lty=2 Line types. If possible.5.sub The font to be used for axis annotation. x and y labels.
font. The default varies with graphics drivers. lines. x and y labels. Affects axis lines as well as lines drawn with lines(). pch can be a character or a number in the range 32:255 representing a character in the current font. text and polygons (filled regions. 3 to italic. An integer which specifies which font to use for text. etc. pch = 0:25) Those from 21 to 25 may appear to duplicate earlier symbols. respectively.Chapter 12: Graphical procedures
75
12. A number from the current palette (see ?palette) or a named colour.

3 Figure margins
A single plot in R is known as a figure and comprises a plot region surrounded by margins (possibly containing axis labels. Negative values give tick marks outside the plotting region.2 Axes and tick marks
Many of R’s high-level plots have axes.
. as a fraction of the size of the plotting region. main and sub-titles.01 Length of tick marks.
12. 1 means always horizontal.5.)
xaxs="r" yaxs="i"
12. and 2 means always perpendicular to the axis. The value is the desired size of text characters (including plotting characters) relative to the default text size. respectively. (S has other styles not implemented in R. Axes have three main components: the axis line (line style controlled by the lty graphics parameter).5.) These components can be customized with the following graphics parameters. Positive numbers measure outside the plot region. and the final component is the distance from the axis position to the axis line (usually zero).main cex.) Choosing a too-small value for this parameter may result in all tick labels being rounded to the same number! las=1 Orientation of axis labels. When tck is small (less than 0.5
Character expansion.Chapter 12: Graphical procedures
76
cex=1. Use tck=0. tck=0. and you can construct axes yourself with the lowlevel axis() graphics function.-1. titles. in text lines. however style "r" leaves a small amount of space at the edges.01 and mgp=c(1.axis cex. The third number is the desired length of axis labels. x and y labels. etc.sub
The character expansion to be used for axis annotation.) and (usually) bounded by the axes themselves.5) the tick marks on the x and y axes are forced to be the same size.0) for internal tick marks. respectively.
mgp=c(3. Axis styles for the x and y axes. 0) Positions of axis components. 0 means always parallel to axis.5. The second component is the distance to the tick labels. The first component is the distance from the axis label to the axis position. in characters (including the decimal point.lab cex. 7. 12) The first two numbers are the desired number of tick intervals on the x and y axes respectively. lab=c(5. the tick marks (which mark off unit divisions along the axis line) and the tick labels (which mark the units.
cex. negative numbers inside. A value of 1 gives grid lines. 1. With styles "i" (internal) and "r" (the default) tick marks always fall within the range of the data.

0
0.Chapter 12: Graphical procedures
77
A typical figure is
−−−−−−−−−−−−−−−−−− −−−−−−−−−−−−−−−−−− −−−−−−−−−−−−−−−−−− −−−−−−−−−−−−−−−−−− −−−−−−−−−−−−−−−−−− −−−−−−−−−−−−−−−−−−
3. 0.
.5. top and right margins. 2. left. respectively.0
y
−1.
mar=c(4. however this may not be enough when many figures share the same page.5 −3.5
mai[2]
−1. Furthermore. 0) Widths of the bottom. measured in inches. 1) Similar to mai. except the measurement unit is text lines. 2.0 x
1. using the postscript() driver with the height=4 argument will result in a plot which is about 50% margin unless mar or mai are set explicitly. When multiple figures are in use (see below) the margins are reduced.0
mar[3]
Plot region
1. The bottom and left margins must be large enough to accommodate the axis and tick labels. The default values chosen for this parameter are often too large. the right-hand margin is rarely needed. 0.
mar and mai are equivalent in the sense that setting one changes the value of the other.5.5
3. the default is chosen without regard to the size of the device surface: for example. and neither is the top margin if no title is being used.0 −3.5
0.0
mai[1]
Margin
Graphics parameters controlling figure layout include:
mai=c(1.

4 Multiple figure environment
R allows you to create an n by m array of figures on a single page. right.5. Values are the positions of the left. the last two are the number of rows and columns in the multiple figure array. 3. The layout in the Figure could have been created by setting mfrow=c(3. 4)/10 Position of the current figure on the page.3. Setting either of these can reduce the base size of symbols and text (controlled by par("cex") and the pointsize of the device). fig=c(4.2)
mfrow=c(3. 1. as a percentage of the page measured
. Each figure has its own margins. 2) Position of the current figure in a multiple figure environment.2. the figure shows the page after four plots have been drawn. mfg=c(2. as shown in the following figure. You can even use different values for the last two numbers than the true values for unequally-sized figures on the same page. the reduction factor is 0.66. In a layout with exactly two rows and columns the base size is reduced by a factor of 0.83: if there are three or more of either rows or columns. 4) Set the size of a multiple figure array. 9. and the array of figures is optionally surrounded by an outer margin.2). bottom and top edges respectively. Set this parameter to jump between figures in the array.
−−−−−−−−−−−−−−− −−−−−−−−−−−−−−− −−−−−−−−−−−−−−− −−−−−−−−−−−−−−− −−−−−−−−−−−−−−− oma[3]
omi[4]
mfg=c(3. the second is the number of columns. mfrow fills by rows.2)
omi[1]
The graphical parameters relating to multiple figures are as follows: mfcol=c(3. The only difference between these two parameters is that setting mfcol causes figures to be filled by column. The first two numbers are the row and column of the current figure.Chapter 12: Graphical procedures
78
12. The first value is the number of rows. 2) mfrow=c(2. 2.

0) Size of outer margins. so you must create them explicitly using oma or omi.off() This ensures that the device finishes cleanly. Text can be added to the outer margins with the mtext() function with argument outer=TRUE. starting with the bottom margin and working clockwise. 0. (Not always available: see its help page. There are no outer margins by default. oma=c(2. For example. for example in the case of hardcopy devices this ensures that every page is completed and has been sent to the printer.” for example) into a form that the particular device can understand. Outer margins are particularly useful for page-wise titles. Before this can begin. however. 0) omi=c(0. however. the first measures in text lines and the second in inches. Produces a bitmap PNG file. best used for image plots.screen() and layout() functions. 3.) Produces a bitmap JPEG file.8. Set this parameter for arbitrary positioning of figures within a page. 0. The purpose of a device driver is to convert graphical instructions from R (“draw a line.) For use with the X11 window system on Unix-alikes
When you have finished with a device.Chapter 12: Graphical procedures
79
from the bottom left corner.6 Device drivers
R can generate graphics (of varying levels of quality) on almost any type of display or printing device. as well as by the grid and lattice packages. use new=TRUE as well (unlike S). (Not always available: see its help page. This is done by starting a device driver. Some commonly-used device drivers are: X11() windows() For use on Windows quartz() For use on Mac OS X postscript() For printing on PostScript printers. There is one such function for every device driver: type help(Devices) for a list of them all. be sure to terminate the device driver by issuing the command > dev.)
. or creating PostScript graphics files. The example value would be for a figure in the bottom right of the page. If you want to add a figure to a current page.
12. issuing the command > postscript() causes all future graphics output to be sent to the printer in PostScript format. pdf() png() jpeg() Produces a PDF file. R needs to be informed what type of device it is dealing with. which can also be included into PDF files. Device drivers are started by calling a device driver function. 0. etc. (This will happen automatically at the normal end of a session. More complicated arrangements of multiple figures can be produced by the split. Like mar and mai.

onefile=FALSE.metafile() [Windows] quartz() [Mac OS X]
postscript() pdf() png() jpeg() tiff() bitmap() . The plot will be in landscape orientation unless the horizontal=FALSE argument is given. This works best when encapsulated PostScript is produced: R always produces conformant output. Many usages of PostScript output will be to incorporate the figure in another document. perhaps for inclusion in a document.2 Multiple graphics devices
In advanced use of R it is often useful to have several graphics devices in use at the same time. they form a numbered sequence with names giving the kind of device at any position.Chapter 12: Graphical procedures
80
12.printer() win..ps". horizontal=FALSE. the command > postscript("file.6.1 PostScript diagrams for typeset documents
By passing the file argument to the postscript() device driver function. width=6. It is important to note that if the file named in the command already exists. you may store the graphics in PostScript format in a file of your choice. When multiple devices are open.6. height=5. pointsize=10)
12. Thus to produce a plot for inclusion use something like > postscript("plot1. The main commands used for operating with multiple devices. This unusual notation stems from S-compatibility: it really means that the output will be a single page (which is part of the EPSF specification). pointsize=10) will produce a file containing PostScript code for a figure five inches high. and you can control the size of the graphic with the width and height arguments (the plot will be scaled as appropriate to fit these dimensions. height=8. and their meanings are as follows: X11() [UNIX]
windows() win. thus extending by one the device list.
. to which graphics output will be sent. This is the case even if the file was only created earlier in the same R session.) For example. but only marks the output as such when the onefile=FALSE argument is supplied. and this is known as the current device. horizontal=FALSE. Each new call to a device driver function opens a new graphics device.eps". This device becomes the current device.. it will be overwritten. Of course only one graphics device can accept graphics commands at any one time.

.print is similar. graphics. such as printing hardcopies...off(k ) Terminate the graphics device at point k of the device list.7 Dynamic graphics
R does not have builtin capabilities for dynamic or interactive graphics.’. such as postscript devices. e.
.list() Returns the number and name of all active devices.. are immediately performed.set(which=k ) Can be used to change the current graphics device to the one at position k of the device list.org/rggobi. Returns the number and label of the device. so that end actions.ggobi.. which=k ) dev. or previous to the current device. described at http://www. For some devices. extensive dynamic graphics facilities are available in the system GGobi by Swayne. package rgl provides ways to interact with 3D plots. for example of surfaces. .ggobi. dev.
12. but the copied device is immediately closed. with extra arguments. depending on how the device was initiated... which=k ) Make a copy of the device k.print(device.. dev. Here device is a device function. if needed. Cook and Buja available from http://www.copy(device.next() dev. dev. Also.. this will either print the file immediately or correctly complete the file for later printing. dev. rotating point clouds or to “brushing” (interactively highlighting) points.off() Terminate all graphics devices on the list. However.g. specified by ‘. except the null device.org/ and these can be accessed from R via the package rggobi. dev. The device at position 1 on the list is always the null device which does not accept graphics commands at all. respectively.prev() Returns the number and name of the graphics device next to. such as postscript.Chapter 12: Graphical procedures
81
dev.

the boot package containing functions from Davison & Hinkley (1997)). See Section “R packages” in R FAQ. for a complete list. The R FAQ contains a list that was current at the time of release. To load a particular package (e. The process of developing packages is described in Section “Creating R packages” in Writing R Extensions. They should be automatically available in any R installation.org/ and its mirrors). and to aid package developers.g. Some (the recommended packages) are distributed with every binary distribution of R. Some packages may be loaded but not available on the search list (see Section 13. They contain the basic functions that allow R to work.. Here. use a command like > library(boot) Users connected to the Internet can use the install.
13. written by many different authors. use > help.packages() functions (available through the Packages menu in the Windows and RAqua GUIs. and the datasets and standard statistical and graphical functions that are described in this manual.packages() and update. we will describe them from a user’s point of view.
.2 Contributed packages and CRAN
There are hundreds of contributed packages for R. To see which packages are currently loaded. Some of these packages implement specialized statistical methods. but the collection of available packages changes frequently. use > search() to display the search list. and others are designed to complement textbooks. To see which packages are installed at your site. and then navigate to the package listing in the Reference section.1 Standard packages
The standard (or base) packages are considered part of the R source code.3 [Namespaces]. Most are available for download from CRAN (http://CRAN. see Section “Installing packages” in R Installation and Administration) to install and update packages. Only when a package is loaded are its contents available.start() to start the HTML help system.R-project. and other repositories such as Bioconductor (http://www. others give access to data or hardware. who are protected from name clashes with other code.
13. page 83): these will be included in the list given by > loadedNamespaces() To see a list of all available help topics in an installed package. This is done both for efficiency (the full list would take more memory and would take longer to search than a subset).bioconductor.Chapter 13: Packages
82
13 Packages
All R functions and datasets are stored in packages.org/). issue the command > library() with no arguments.

The colon operators described above will also cause automatic loading of the associated package. In the example above. Users are more likely to use the getAnywhere() function. and currently all of the base and recommended packages do except the datasets package. Only functions that are exported from the package can be retrieved in this way. There are two operators that work with namespaces. Packages are often inter-dependent.Chapter 13: Packages
83
13.3 Namespaces
Packages can have namespaces. and they provide a way to refer to an object within a particular package. and breaking every function that tries to transpose a matrix. they prevent functions from breaking when a user (or other package writer) picks a name that clashes with one in the package. The double-colon operator :: selects definitions from a particular namespace. which searches multiple packages. because it is defined in the base package. When packages with namespaces are loaded automatically they are not added to the search list. but users might define their own function named t. The triple-colon operator ::: may be seen in a few places in R code: it acts like the double-colon operator but also allows access to hidden objects. For example. Namespaces do three things: they allow the package writer to hide functions and data that are meant only for internal use.
. Namespaces prevent the user’s definition from taking precedence. t() is the transpose function in R. and loading one may cause others to be automatically loaded. the transpose function will always be available as base::t.

Iconify the help window and move on to the next part. A graphics window will appear automatically. y) Standard point plot. The R program begins. $R Start R as appropriate for your platform. (Clean up). data=dummy. (Within R.rnorm(50) y <. start your windowing system. Login.lm(y ~ x. with a banner. Many features of the system will be unfamiliar and puzzling at first. but this puzzlement will soon disappear. dummy <.
.data.lowess(x. weight=1/w^2) summary(fm1) Since we know the standard deviations. y) Plot the points in the plane. x and y.
x <.lm(y ~ x. 2. w <. we can do a weighted regression. fm <.1 + sqrt(x)/2 A ‘weight’ vector of standard deviations. y) Make a nonparametric local regression function. You should briefly explore the features of this facility with the mouse. ls() rm(x. x <. we are modelling y dependent on x. fm1 <.)
help. . y) See which R objects are now in the R workspace. plot(x.start() Start the HTML interface to on-line help (using a web browser available at your machine).Appendix A: A sample session
84
Appendix A A sample session
The following session is intended to introduce to you some features of the R environment by using them. the prompt on the left hand side will not be shown to avoid confusion. . . data=dummy) summary(fm) Fit a simple linear regression and look at the analysis. attach(dummy) Make the columns in the data frame visible as variables. y= x + rnorm(x)*w) dummy Make a data frame of two columns. . Remove objects no longer needed.frame(x=x. plot(x. lrf <.and y-coordinates. and look at it. 20).rnorm(x) Generate two pseudo-random normal vectors of x. With y to the left of the tilde.1:20 Make x = (1.

file("data". The next section will look at data from the classical experiment of Michaelson and Morley to measure the speed of light. abline(coef(fm1). This dataset is available in the morley object. lty=3) The true regression line: (intercept 0. filepath <. x. plot(Expt. xlab="Experiment No. main="Residuals vs Fitted") A standard regression diagnostic plot to check for heteroscedasticity.
plot(fitted(fm). suitably coded. detach() Remove data frame from the search path.factor(mm$Run) Change Expt and Run into factors. file. ylab="Residuals". mm <. Can you see it? qqnorm(resid(fm).Appendix A: A sample session
85
lines(x.read. abline(coef(fm)) Unweighted regression line. slope 1).table(filepath) mm Read in the Michaelson and Morley data as a data frame. data=mm) summary(fm) Analyze as a randomized block.tab" .system. xlab="Fitted values".aov(Speed ~ Run + Expt. main="Residuals Rankit Plot") A normal scores plot to check for skewness. and look at it. package="datasets") filepath Get the path to the data file. fm1.) rm(fm. mm$Expt <. dummy) Clean up again. "morley. resid(fm). There are five experiments (column Expt) and each has 20 runs (column Run) and sl is the recorded speed of light. with ‘runs’ and ‘experiments’ as factors.show(filepath) Optional. kurtosis and outliers. main="Speed of Light Data".factor(mm$Expt) mm$Run <. col = "red") Weighted regression line.table function. fm <. (Not very useful here. Speed. abline(0. lrf. Look at the file. attach(mm) Make the data frame visible at position 3 (the default). 1.") Compare the five experiments with simple boxplots. but we will read it to illustrate the read. lrf$y) Add in the local regression.
.

One method would be to take complex numbers with standard normal real and imaginary parts .1). w <. xlab="x". w. and to map any outside the circle onto their reciprocal. rm(th. w) . xlim=c(-1. 1/w. . ylab="y") lines(z) The second method uses the uniform distribution. q() Quit the R program.Appendix A: A sample session
87
w <.1).1).1).
.sqrt(runif(100))*exp(2*pi*runif(100)*1i) plot(w. You will be asked if you want to save the R workspace. you probably do not want to save it.xlab="x".ifelse(Mod(w) > 1. plot(w. . pch="+".rnorm(100) + rnorm(100)*1i Suppose we want to sample points within the unit circle. The points should now look more evenly spaced over the disc. ylim=c(-1. . z) Clean up again. and for an exploratory session like this. ylim=c(-1. xlim=c(-1. ylab="y") lines(z) All points are inside the unit circle. pch="+". but the distribution is not uniform. . w <.

site’ is used (if it exists). (See help("Startup") for a precise description. or reside in ‘. • Then R searches for the site-wide startup profile unless the command line option ‘--no-site-file’ was given. The startup mechanism is as follows (see also the on-line help for topic ‘Startup’ for more information. Most options control what happens at the beginning and at the end of an R session. as a wrapper to various R tools (e. The name of this file is taken from the environment variable R_PROFILE_USER.RData’. it is executed. a file called ‘. The name of this file is taken from the value of the R_PROFILE environment variable.Renviron’ in the current or in the user’s home directory (in that order) are searched for. The user file is the one pointed to by the environment variable R_ ENVIRON_USER if this is set. there are options for controlling the memory available to the R process (see the on-line help for topic ‘Memory’ for more information). This function (as well as . R searches for user and site files to process for setting environment variables. via the R CMD interface. In addition. or.Appendix B: Invoking R
88
Appendix B Invoking R
B. • It also loads a saved image from ‘. Users will not normally need to use these unless they are trying to limit the amount of memory used by R. R accepts the following command-line options. for processing files in R documentation format or manipulating add-on packages) which are not intended to be called “directly”. • Then. files ‘.g. R searches for a user profile and sources it. R_PRINTCMD (the default print command) and R_LIBS (specifies the list of R library trees searched for add-on packages). You need to ensure that either the environment variable TMPDIR is unset or it points to a valid place to create temporary files and directories. if this is unset.First exists. the default ‘R_HOME /etc/Rprofile. otherwise. if a function .1 Invoking R from the command line
When working in UNIX or at a command line in Windows.Last which is executed at the end of the R session) can be defined in the appropriate startup profiles.site’ is used if this exists.Rprofile’ in the current directory or in the user’s home directory (in that order) is searched for. ‘R_HOME /etc/Renviron. unless ‘--no-init-file’ was given.) Variables you might want to set include R_PAPERSIZE (the default paper size). and the section below for some Windows-specific details).
. ‘--help’ ‘-h’ Print short help message to standard output and exit successfully. the command ‘R’ can be used both for starting the main R program in the form R [options] [<infile] [>outfile]. • Finally.. if unset. These files should contain lines of the form ‘name =value ’. The name of the site file is the one pointed to by the environment variable R_ENVIRON. If that variable is unset. • Unless ‘--no-environ’ was given.RData’ if there is one (unless ‘--no-restore’ or ‘--no-restore-data’ was specified).

Apart from the front-end shell script and the man page. this also includes ‘--no-Rconsole’. packages. ‘--no-environ’ Do not read any user file to set environment variables. ‘--restore’ ‘--no-restore’ ‘--no-restore-data’ Control whether saved images (file ‘. ‘-e expression ’ Use expression as an input line.
‘--save’ ‘--no-save’ Control whether data sets should be saved or not at the end of the R session. in non-interactive use one of these must be specified or implied by some other option (see below). Implies ‘--no-save’ unless ‘--save’ has been set. ‘--no-init-file’ Do not read the user’s profile at startup.Appendix B: Invoking R
89
‘--version’ Print version information to standard output and exit successfully. the user is asked for the desired behavior when ending the session with q().) into this directory. ‘--no-init-file’ and ‘--no-restore’. but can be set by the environment variable R_HISTFILE) should be restored at startup or not. ‘RHOME’ Print the path to the R “home directory” to standard output and exit successfully. ‘-f file ’ ‘--file=file ’ Take input from file: ‘-’ means stdin. One or more ‘-e’ options can be used.RData’ in the directory where R was started) should be restored at startup or not. but not together with ‘-f’ or ‘--file’. Implies ‘--no-save’ unless ‘--save’ has been
.) ‘--no-restore-history’ Control whether the history file (normally file ‘. ‘--no-site-file’ Do not read the site-wide profile at startup. This needs to be an encoding known to iconv: see its help page. The default is to restore. Under Windows. ‘--encoding=enc ’ Specify the encoding to be assumed for input from the console or stdin. The default is to restore. ‘--no-environ’. ‘--vanilla’ Combine ‘--no-save’. ‘--no-Rconsole’ (Windows only) Prevent loading the ‘Rconsole’ file at startup.Rhistory’ in the directory where R was started. (‘--no-restore’ implies all the specific ‘--no-restore-*’ options. If neither is given in an interactive session. etc. R installation puts everything (executables. ‘--no-site-file’.

(There is a limit of 10.) ‘--no-readline’ (UNIX only) Turn off command-line editing via readline. See Appendix C [The command-line editor].9.expand. ‘-f’ or ‘--file’ asserts non-interactive use even if ‘--interactive’ is given. or ‘k’. ‘--quiet’ ‘--silent’ ‘-q’ Do not print out the initial copyright and welcome messages.5Gb on some 64-bit versions of Windows. ‘--min-vsize=N ’ ‘--max-vsize=N ’ Specify the minimum or maximum amount of memory used for variable size objects by setting the “vector heap” size to N bytes. A cons cell takes 28 bytes on a 32-bit machine. ‘M’. Command-line editing is enabled by default interactive use (see ‘--interactive’). for more information. (The default is to deduce that R is being run interactively if and only if ‘stdin’ is connected to a terminal or pty. Currently the maximum value accepted is 100000.000 bytes on the total length of expressions used in this way.Appendix B: Invoking R
90
set. See the previous option for details on N. ‘Mega’ (2^20). ‘K’. This is useful when running R from within Emacs using the ESS (“Emacs Speaks Statistics”) package. or regular ‘kilo’ (1000). ‘--min-nsize=N ’ ‘--max-nsize=N ’ Specify the amount of memory used for fixed size objects by setting the number of “cons cells” to N. ‘--max-mem-size=N ’ (Windows only) Specify a limit for the amount of memory to be used both for R objects and working areas. This is set by default to the smaller of 1. ‘--slave’ Make R run as quietly as possible. This defaults to 10000.) Using ‘-e’. and usually 56 bytes on a 64-bit machine. Here. This option is intended to support programs which use R to compute results for them. but can be increased to allow large and complicated calculations to be done.5Gb on versions of Windows that support 3Gb per process and have the support enabled: see the ‘rw-FAQ’ Q2. (computer) ‘Kilo’ (2^10). 3.
1
2. and must be between 32Mb and the maximum allowed on that version of Windows.
. N must either be an integer or an integer ending with ‘G’.
‘--interactive’ (UNIX only) Assert that R really is being run interactively even if input has been redirected: use if input is from a FIFO or pipe and fed from an interactive program. ‘--max-ppsize=N ’ Specify the maximum size of the pointer protection stack as N locations. This option also affects tilde-epansion: see the help for path. page 95. meaning ‘Giga’ (2^30).5Gb1 and the amount of physical RAM in the machine. It implies ‘--quiet’ and ‘--no-save’.

provided that ‘Tcl/Tk’ support is available. but not intended to be called “directly”. Currently. Warning and error messages are sent to the error channel (stderr).
A Convert Rd format to various other formats. ‘--debugger=name ’ ‘-d name ’ (UNIX only) Run R through debugger name. Post-process R profiling files. Check add-on packages. Currently. possible values for type are ‘X11’ (the default) and. and should instead be given when starting the R executable from inside the debugger. further command line options are disregarded. (For back-compatibility. including asserting interactive use without the command-line editor. ‘x11’ and ‘tk’ are accepted. The general form is R CMD command args where command is the name of the tool and args the arguments passed on to it. (UNIX only) Compile files for use with R. Build (that is. R code uses this option to control the printing of diagnostic messages. package) add-on packages. (UNIX only) Front-end for creating executable programs. including HTML.
Note that input and output can be redirected in the usual way (using ‘<’ and ‘>’). Build shared library for dynamic loading. ‘Tk’. and in particular set R’s option verbose to TRUE.
.
‘--verbose’ Print more information about progress. but the line length limit of 4095 bytes still applies. Remove add-on packages. plain text. the following tools are available. ‘--gui=type ’ ‘-g type ’ (UNIX only) Use type as graphical user interface (note that this also includes interactive graphics). The command R CMD allows the invocation of various tools which are useful in conjunction with R.Appendix B: Invoking R
91
‘--ess’
(Windows only) Set Rterm up for use by R-inferior-mode in ESS. For most debuggers (the exceptions are valgrind and recent versions of gdb). Install add-on packages. Rd2txt can be used as shorthand for Rd2txt -t txt. L TEX. and extracting the examples.) ‘--args’ This flag does nothing except cause the rest of the command line to be skipped: this can be useful to retrieve values from it with commandArgs(TRUE). BATCH COMPILE SHLIB INSTALL REMOVE build check LINK Rprof Rdconv Rd2txt Run R in batch mode.

that gives the home directory. Rd2pdf can be used as shorthand for Rd2dvi --pdf. (This mechanism is used for drag-and-drop and file association with RGui.exe). but also works for Rterm. Environment variables can be supplied as ‘name =value ’ pairs on the command line.2 Invoking R under Windows
There are two ways to run R under Windows. if the environment variable HOME is defined.Appendix B: Invoking R
92
Rd2dvi Rd2pdf Sd2Rd Stangle Sweave Rdiff config javareconf
Convert Rd format to DVI/PDF. If the named file does not exist it sets the working directory if the parent directory exists. You need to ensure that either the environment variables TMPDIR. there is a console-based GUI (Rgui. In addition. It first tries to use the Windows "personal" directory (typically C:\Documents and Settings\username\My Documents in Windows XP). Convert S documentation to Rd format.com or a more capable shell).exe.
B. After those two user-controllable settings. you can use R CMD cmd args for any other executable cmd on the path: this is useful to have the same environment as R or the specific commands run under. Under a Unix-alike.)
.) For interactive use. If there is an argument ending ‘. for example to run ldd or pdflatex. If the environment variable R_USER is defined. invoking by R. Failing all those.exe or command. (Unix only) Update the Java configuration variables
Use R CMD command --help to obtain usage information for each of the tools accessible via the R CMD interface. the home directory is taken to be the starting directory. the methods described in the previous section may be used. The startup procedure under Windows is very similar to that under UNIX. cmd. that gives the home directory.exe or more directly by Rterm. as this is not always defined on Windows. Within a terminal window (e. (These are principally intended for batch use.g. and environment variables HOMEDRIVE and HOMEPATH are defined (and they normally are) these define the home directory. TMP and TEMP are either unset or one of them points to a valid place to create temporary files and directories.RData’ (in any case) it is interpreted as the path to the workspace to be restored: it implies ‘--restore’ and sets the working directory to the parent of the named file. Next.exe. if cmd is ‘perl’ or ‘awk’ it is replaced by the full path to the Perl or AWK command found when R was configured. Extract S/R code from Sweave documentation Process Sweave documentation Diff R output ignoring headers etc Obtain configuration information. R tries to find system defined home directories. but references to the ‘home directory’ need to be clarified. If that fails.exe.

pl’) with several environment variables set appropriately.tex’. ‘--debug’ Enable the “Break to debugger” menu item in Rgui.R. but the startup and current working directory are set as the user’s home directory unless a different startup directory is given in the Preferences window accessible from within the GUI. the methods described in the first subsection apply. which can be invoked by Rscript foo.exe’.exe’ on your path. with the path to R’s ‘share/texmf’ macros appended to A TEXINPUTS. The ‘home directory’ is the one inside the R.
In Windows with R CMD you may also specify your own ‘. R_OSTYPE. You can pass parameters to scripts via additional arguments on the command line: for example R CMD BATCH --args arg1 arg2 foo. It is a standard double-clickable Mac OS X application.R arg1 arg2 and this can also be used to write executable script files like (at least on Unix-alikes. the recommended way is to use R CMD BATCH foo.R & runs a background job.
B.bat’.exe mydoc A will run L TEX on ‘mydoc.sh’ or ‘. PATH.framework. ‘. if you already have ‘latex.pl’ file.app window by invoking R. There is also console-based GUI (R. and trigger a break to the debugger during command line processing. (Unfortunately. PERL5LIB.commandArgs(TRUE) This is made simpler by the alternative front-end Rscript. If you want to run this in the background or as a batch job use OS-specific facilities to do so: for example in most shells on Unix-alike OSes R CMD BATCH foo. ‘--mdi’ ‘--sdi’ ‘--no-mdi’ Control whether Rgui will operate as an MDI program (with multiple child windows within one main window) or an SDI application (with multiple toplevel windows for the console. ‘.R’ of R commands.exe. including R_HOME. It will be run under the appropriate interpreter (Perl for ‘.)
B. The command-line setting overrides the setting in the user’s ‘Rconsole’ file.app) that by default is installed in the Applications folder on your system. The startup procedure under Mac OS X is very similar to that under UNIX.4 Scripting with R
If you just want to run a file ‘foo.Appendix B: Invoking R
93
The following additional command-line options are available when invoking RGui. and in some Windows shells)
. R_CMD. this does not help with the MiKTeX build of L TEX. R_VERSION. then R CMD latex. For example. Within a Terminal. graphics and pager). and TEXINPUTS.R & will pass arguments to a script which can be retrieved as a character vector by args <.3 Invoking R under Mac OS X
There are two ways to run R under Mac OS X.

77 3.03 28.70 3.50 3. .70 3.scan(n=24) 2.. e..90 3. If you do not wish to hardcode the path to Rscript but have it in your path (which is normally the case for an installed R except on Windows).. It is commonplace to write R scripts with segments like chem <..40 2. q(status=<exit status code>) If this is entered into a text file ‘runfoo’ and this is made executable (by chmod 755 runfoo).40 2.
. Another way to write executable script files (suggested by Fran¸ois Pinard) is to use a c here document like #!/bin/sh [environment variables can be set here] R --slave [other options] <<EOF R program goes here.Appendix B: Invoking R
94
#! /path/to/Rscript args <.).40 3.37 3. use "stdin" as a file connection. This writes R output to ‘stdout’ and ‘stderr’.60 3.20 5.28 3.40 3.70 2. it can be invoked for different arguments by runfoo arg1 arg2 For further options see help("Rscript").95 3. At least in Bourne and bash shells.20 3..50 2.10 3.70 and stdin() refers to the script file to allow such traditional usage. Very short scripts can be passed to Rscript on the command-line via the ‘-e’ flag. and this can be redirected in the usual way for the shell running the command.70 2. use #! /usr/bin/env Rscript . One thing to consider is what stdin() refers to.40 2.commandArgs(TRUE) . EOF but here stdin() refers to the program source and "stdin" will not be usable..80 2. If you want to refer to the process’s ‘stdin’.g.. the #! mechanism does not allow extra arguments like #! /usr/bin/env Rscript --vanilla.03 3. scan("stdin"..

including the erroneous lines. In vi mode character insertion mode is started by M-i or M-a. If your terminal does not have a META key. changed if necessary. The ESC character sequences are also allowed on terminals with real Meta keys. and commands in your history may be recalled. Meta characters. displacing any characters to the right of the cursor.1 Preliminaries
When the GNU readline library is available at the time R is configured for compilation under UNIX. an inbuilt command line editor allowing recall. In Emacs-style command-line editing any straight typing you do while in this editing phase causes the characters to be inserted in the command you are editing. occasionally the ‘Windows’ key. Control characters. Many of these use either Control or Meta characters. Thus. such as Control-m. and written as M-b in the following. editing and re-submission of prior commands is used. see the URL http://ESS.Appendix C: The command-line editor
95
Appendix C The command-line editor
C.R-project. you could type ESCB.2 Editing actions
The R program keeps a history of the command lines you type. characters are typed and insertion mode is finished by typing a further ESC. and the file ‘README. Other editing actions are summarized in the following table. It can be disabled (useful for usage with ESS1 ) using the startup option ‘--no-readline’.
C. you can also use the up and down arrow keys instead of C-p and C-n. Go to the next command (forwards in the history). When using R with readline capabilities. Windows versions of R have somewhat simpler command-line editing: see ‘Console’ under the ‘Help’ menu of the GUI.exe. are typed by holding down META2 and pressing B.
. are obtained by holding the CTRL down while you press the M key. respectively. and re-submitted as new commands.
C. such as Meta-b.
On most terminals.org On a PC keyboard this is usually the Alt key.
1 2
The ‘Emacs Speaks Statistics’ package. you can still type Meta characters using two-character sequences starting with ESC. the functions described below are available.Rterm’ for command-line editing under Rterm. Find the last command with the text string in it. to enter M-b. Note that other versions of readline exist and may be used by the inbuilt command line editor: this may happen on Mac OS X. Note that case is significant for Meta characters. Pressing the RET command at any time causes the command to be re-submitted.3 Command-line editor summary
Command recall and vertical motion
C-p C-n C-r text Go to the previous command (backwards in the history). and are written as C-m below.

and “save” it. respectively.Appendix C: The command-line editor
96
Horizontal motion of the cursor
C-a C-e M-b M-f C-b C-f Go to the beginning of the command. Delete the character under the cursor.
The final RET terminates the command line editing sequence. Change the rest of the word to upper case.
Editing and re-submission
text C-f text DEL C-d M-d C-k C-y C-t M-l M-c RET Insert text at the cursor.
On most terminals. Insert (yank) the last “saved” text here. Go back one word. Delete the previous character (left of the cursor). Go to the end of the line. Transpose the character under the cursor with the next. Append text after the cursor. Re-submit the command to R. Go back one character. Go forward one character. you can also use the left and right arrow keys instead of C-b and C-f.
. Delete the rest of the word under the cursor. Go forward one word. Delete from cursor to end of command. and “save” it. Change the rest of the word to lower case.