possum n. 1 Any of many chiefly herbivorous, long-tailed, tree-dwelling, mainly Australian marsupials, some of which are gliding animals (e.g. brush-tailed possum, flying possum). 2 a mildly scornful term for a person. 3 an affectionate mode of address. From the Australian Oxford Paperback Dictionary, 2nd ed, 1996.

45
4 0

50

55

34
3 2

36

38

40

42

ii

Introduction...........................................................................................................................................................1 The R System..................................................................................................................................................1 The Look and Feel of R ..................................................................................................................................1 The Use of these Notes ...................................................................................................................................2 The R Project ..................................................................................................................................................2 Web Pages and Email Lists.............................................................................................................................2 Datasets that relate to these notes....................................................................................................................2 _________________________________________________________________________ .......................2 1. Starting Up ........................................................................................................................................................3 1.1 1.2 1.3 Getting started under Windows.................................................................................................................3 Use of an Editor Script Window................................................................................................................4 A Short R Session .....................................................................................................................................5 Entry of Data at the Command Line...................................................................................................6

or to know when naïve methods are inadequate.1. there can be recourse to the huge range of more advanced abilities that R offers. allowing the calculations such as 2+3. for Unix and for Macintosh OS X. and exposed to higher standards of scrutiny than most other systems. visualize and manipulate data. Adaptation of available abilities allows even greater flexibility. Chapters 1 and 2 indicate the range of abilities that are immediately available to novice users.) It is available through the Comprehensive R Archive Network (CRAN). likely to be an educational experience. Simple calculations and analyses can be handled straightforwardly.
The R System
R implements a dialect of the S language that was developed at AT&T Bell Laboratories by Rick Becker. there are no pointers.5 has an example. The skills needed for the computing are not on their own enough. to queries. Users who want a point and click interface should investigate the R Commander (Rcmdr package) interface. that are tightly linked with its analytic abilities. Books that provide a more extended commentary on the methods illustrated in these examples include Maindonald and Braun (2003). most computation is handled using functions. It is a community that is sensitive to the potential for misuse of statistical techniques and suspicious of what might appear to be mindless use. Anyone with a contrary view may care to consider whether a butcher’s meat-cleaving skills are likely to be adequate for effective animal (or maybe human!) surgery. The citation for John Chambers’ 1998 Association for Computing Machinery Software award stated that S has “forever altered how people analyze. Important differences are that R has no header files. It is often possible and desirable to operate on objects – vectors. This largely avoids the need for explicit loops. John Chambers and Allan Wilks. or 3^11. Be grateful for whatever help is given.1
Introduction
These notes are designed to allow individuals who have a basic grounding in statistical methodology to work through examples that demonstrate the use of R for a range of types of data manipulation. arrays. While R is as reliable as any statistical software that is available. Expect scepticism of the use of models that are not susceptible to some minimal form of data-based validation. there are traps that call for special care. with a limited tradition of experience of the limitations and pitfalls. Beyond this. (There are are older versions of R that support 8. Section 2. Versions of R are available. Some of the model fitting routines are leading edge. graphical presentation and statistical analysis.
. users have no right to expect attention. The R system is developing rapidly.
1
The structure of an R program has similarities with programs that are written in C or in its successors C++ and Java. If simple methods prove inadequate. and vectors of text strings can be defined and manipulated directly. for 32-bit versions of Microsoft Windows for Linux. The R community is widely drawn. check each step with care. Hurrah for the R development team!
The Look and Feel of R
R is a functional language. The action of quitting from an R session uses the function call q(). Neither R nor any other statistical system will give the statistical expertise needed to use sophisticated abilities.” The R project enlarges on the ideas and insights that generated the S language. New features and abilities appear every few months. and especially when there is some element of complication. from application area specialists as well as statistical specialists. Experience with the use of R is however. at no cost. most declarations are implicit. on the R-help list or elsewhere.6 and 9. Because R is free.1 There is a language core that uses standard forms of algebraic notation. leading to clearer code. Whatever the statistical system. more than with most systems. lists and so on – as a whole. Web addresses are given below. The implementation of R uses a computing model that is based on the Scheme dialect of the LISP language. Here are points that potential users might note: R has extensive and powerful graphics abilities.

edu. or for adaptation to other systems. go to http://cran. and systems that are specifically designed for such manipulation may be preferable. CSIRO). and provided several of the data sets and a number of the exercises. and for other information. typically. Note that data structure is. a regression with 500. go to http://wwwmaths.6 explains how to access the datasets. Development of R is now overseen by a `core team’ of about a dozen people.R.ac. etc.aarnet. Andreas Ruckstuhl (Technikum Winterthur Ingenieurschule. Switzerland) and John Braun (University of Western Ontario) gave me exemplary help in getting the earlier S-PLUS version of this document somewhere near shipshape form.stat. written by the R Development Core Team. which will usually appear within a matter of weeks. for copies of various versions of R. Reported bugs are commonly fixed in the next minor-minor release.g. to various scientists named in the notes who have allowed me to use their data. The r-help email list gives access to an informal support network that can be highly effective..) using the notes as commentary.
_________________________________________________________________________
Jeff Wood (CMIS.anu. an even more important issue for large data sets than for small data sets. The R language environment is designed to facilitate the development of new scientific computational tools. R’s routines can readily process data sets that by historical standards seem large – e.R-project. to open a display file window and transfer the commands across from the that window. are available from the web site for the R project at http://www. Source-code is available for inspection.
Datasets that relate to these notes
Copy down the R image file usingR. if computational efficiency demands. Exposing code to the critical scrutiny of highly expert users has proved an extremely effective way to identify bugs and other inadequacies.edu.RData from http://wwwmaths. Computerintensive components can.edu. Like Linux. I take full responsibility for the errors that remain. Australian users may wish to go directly to http://mirror.r-project. ch3-4.au/~johnm/r/dsets/ Section 1.anu. Datasets are also available individually.au/pub/CRAN There is no official support for R. as demonstrated in Section 1. be handled by a call to a function that is written in the C or Fortran language. and as a result of continuing improvements in the efficiency of R’s coding and memory management. and of other lists that serve the R community.auckland. Readers of these notes may find it helpful to have available for reference the document: “An Introduction to R”. widely drawn from different institutions worldwide.org and find the nearest CRAN (Comprehensive R Archive Network) mirror site. The packages give access to up-to-date methodology from leading statistical and other researchers.
. John Braun gave valuable help with proofreading. repeated smaller analyses with subsets of the total data may give insight that is not available from a single global analysis. on a Unix machine with 2GB of memory. Additionally. the main issue is often manipulation of data. Novice users will notice small but occasionally important differences between the S dialect that R implements and the commercial S-PLUS implementation of S. Those who write substantial functions and (more importantly) packages will find large differences. With very large datasets. R is an “open source” system. Details of the R-help list.R. Under Windows an alternative to typing the commands at the console is. and to elicit ideas for enhancement. hosted at http://www.nz/r-downunder.org/ Note also the Australian and New Zealand list. I am grateful.
The R Project
The initial version of R was developed by Ross Ihaka and Robert Gentleman. supplied with R distributions and available from CRAN sites.2.
Web Pages and Email Lists
For a variety of official and contributed documentation. With the very large address spaces now possible. also.au/~johnm/r/dsets/individual-dsets/ A number of the datasets are now available from the DAAG or DAAGxtras packages.000 cases and 100 variables is feasible.2 The Use of these Notes
The notes are designed so that users can run the examples in the script files (ch1-2. both from the University of Auckland. Email archives can be searched for questions that may have been previously answered.

i. type 2+2 and press the Enter key. Here is what appears on the screen:
> 2+2 [1] 4 >
Here the result is 4.
The command line prompt. see the README file that is included with the R distribution. follow the installation instructions appropriate to the operating system. a little strangely. right click|rename on the copy to rename it3. In interactive use under Microsoft Windows there are several ways to input commands to R. Figure 1 shows the command window as it appears when version 2. perhaps the name itself. and follow instructions. Users of Microsoft Windows may wish to create a separate icon for each such working directory.1
Getting started under Windows
Click on the R icon. click on its icon to start installation. right click|paste to place a copy on the desktop.0. then left click on the copy menu item”.3
1. 1: The upper left portion of the R console (command line) window. Figures 1 and 2 demonstrate two of the possibilities.
Fig. For example. Here. It pays to have a separate working directory for each major project. and then finally go to right click|properties to set the Start in directory to be the working directory that was set up earlier.
.0 of R. Either or both of the following may be used at the user’s discretion.
1.
2 3
This is a shortcut for “right click.e. the >. it.exe from the relevant base directory on the nearest CRAN site.
Enter the name of your choice into the name field. For ease of remembering. Installation is now especially straightforward for Windows users. Packages that do not come with the base distribution must be downloaded and installed separately. there is just one element. Starting Up
R must be installed on your system! If it is not. Then right click|copy2 to copy an existing R icon. For more details. This document mostly assumes that users will type commands into the command window. Copy down the latest SetupR. The > indicates that R is ready for another command. choose the icon that corresponds to the project that is in hand. is an invitation to start typing in your commands. for version 2. immediately after starting up.0.0 of R has just been started. choose a name that closely matches the name of the workspace directory. First create the directory. For this demonstration I will click on my r-notes icon. The[1] says. “first requested element will follow”. at the command line prompt. Or if there is more than one icon.

or new text added. click on Open script…. The icons are. 3: This shows the five icons that appear when the focus is on a script file window. or to click on the in the top right hand corner of the R window.R. Click on the `Run line or selection’ icon.
Fig. If a new blank window is required. 2: The focus is on an R display file window. The text in a script file window can be edited. There will be a message asking whether to save the workspace image. the standard form of input is the command line interface. In Figure 2 the file was firstSteps. which have a somewhat similar set of icons but do not allow editing. starting from the left: Open script. 2 and 3.
4
This requires emacs.4
For later reference.
Fig. Run line or selection. To get a script file window. Under Microsoft Windows. go to the File menu.2
Use of an Editor Script Window
The screen snapshot in Figure2 shows a script file window. with the console window in the background. Highlight the commands that are intended for input to R. you will be asked for the name of a file whose contents are then displayed in the window. To load an existing file. are another possibility. note that the exit or quit command is
> q()
Alternatives to the use of q() are to click on the File menu and then on Exit. or that have been typed or copied into the window. Look under Software|Other on the CRAN page. click on New script. Display file windows. Both are free. Save script. which is the middle icon of the script file editor toolbar in Figs. Under both Microsoft Windows and Linux (or Unix). Return focus to console. and ESS which runs under emacs. and Print. Clicking Yes (the safe option) will save all the objects that remain in the workspace – any that were there at the start of the session and any that have been added since. to send commands to R.
1.
Under Unix. a further possibility is to run R from within the emacs editor4.
. This allows input to R of statements from a file that has been set up in advance.

5 4 6 11 21 3 8 11 17 38 4941 6182 6836 7579 9640
We will read into R a file that holds population figures for Australian states and territories.
.txt”. and total population. is a “plugin” for WinEdt. It means “is assigned to”. Use of header=TRUE causes R to use= the first line to get header information for the columns.5
attractive options are to use either the R-WinEdt utility that is designed for use with the shareware WinEdt editor.
Figure 4: Population of Australian Capital Territory. header=TRUE)
The <. for WinEdt R-WinEdt and various other editors that work with R. data=austpop. at various times since 1917. a data frame. in R parlance. then using this file to create a graph." "Qld" "SA" "WA" "Tas. 306 392 457 502 688 879 193 211 233 257 326 375 NT ACT Aust.read. which is free. Qld 683 873 SA 440 565 WA Tas. If column headings are not included in the file. . Here is a plot (Figure 4) of the Australian Capital Territory (ACT) population between 1917 and 1997 (Figure 4).
5
The R-WinEdt utility.table(“d:/austpop. look under Software|Other on the CRAN web page." "Vic. year has a lower case y. Now type in austpop at the command line prompt.txt on a disk in drive d:
> austpop <. or to use the free tinn-R editor5. For links to the relevant web pages. .
The object austpop is. at various times between 1917 and 1997.is a left diamond bracket (<) followed by a minus sign (-). 5 4 3 8 4941 6182 1 1917 1904 1409 2 1927 2402 1727 . Data frames that consist entirely of numeric data have the same form of rectangular layout as numeric matrices. Qld 683 873 993 SA 440 565 589 646 873 WA Tas. The data in the file are:
1917 1904 1409 1927 2402 1727 1937 2693 1853
1947 2985 2055 1106 1957 3625 2656 1413
1967 4295 3274 1700 1110
62 103 11799
1977 5002 3837 2130 1286 1204 1987 5617 4210 2675 1393 1496 1997 6274 4605 3401 1480 1798
415 104 214 14192 449 158 265 16264 474 187 310 18532
The following reads in the data from the file austpop." "NT"
Note that in the version of the dataset that is in the DAAG package. Code to plot the graph is:
plot(ACT ~ Year.
1. 306 392 193 211 NT ACT Aust. the argument can be omitted. pch=16)
We first of all remind ourselves of the column names:
> names(austpop) [1] "Year" [9] "ACT" "NSW" "Aust.3
Year
A Short R Session
NSW Vic. displaying the object on the screen:
> austpop Year NSW Vic.

182.edit(elasticband)
Figure 5: Editor window. for each amount by which an elastic band is stretched over the end of a ruler. type
elasticband <.”). Command alternatives to the default use of a space are sep=". users can control over the choice of missing value character or characters.173.141. The default is header=FALSE.54.data.
. This last choice makes tabs separators.strings=". distance=c(148. data input will fail. the distance that the band moved when released: stretch distance 46 148 54 182 48 173 50 166 44 109 42 141 52 166
The function data. optionally various parameters additional to the file name that holds the data.3.3 Options for read. There are several variants of read. specify na. Similarly.table() takes. replace ‘Year’ by ‘year’
The option pch=16 sets the plotting character to a solid black dot.3.166.2 Entry and/or editing of data in an editor window
To edit the data frame elasticband in a spreadsheet-like format.1
Entry of Data at the Command Line
A data frame is a rectangular array of columns of data.fields() to report the number of fields that were identified on each separate line of the file. We can specify more informative axis labels. In addition users can specify the separator character or characters. showing the data frame elasticband.
1.table() detects that lines in the input file have different numbers of fields.table() that differ only in having different default parameter settings.csv().table()
The function read. and so on." and sep="\t".48.42.frame() can be used to input these (or other) data directly at the command line.109. Note in particular read. pch=16) # For the DAAG version. which by default is NA.52). If read.frame(stretch=c(46. This plot can be improved greatly. It is then often useful to use the function count.44. Specify header=TRUE if there is an initial row of header names. with an error message that draws attention to the discrepancy. change size of the text and of the plotting symbol. and both columns will be numeric. We will give the data frame the name elasticband:
elasticband <.166))
1.
1.". If the missing value character is a period (“.50. Here we will have two columns.6
A simple way to get the plot is:
> plot(ACT ~ Year. The following data gives.3. data=austpop. which has settings that are suitable for comma delimited (csv) files that have been generated from Excel spreadsheets.

search() and apropos() can be a huge help in finding what one wants. Anything that follows a # on the command line is taken as comment and ignored by R.
. By default.6
The Loading or Attaching of Datasets
The recommended way to access datasets that are supplied for use with these notes is to attach the file usingR.) > apropos(“matrix”)
(This lists all function names that include the text “matrix”.) as the separator. and then use key words to do a search..RData. from the datasets or MASS packages) should then be available without need for any further action. To get help on a specific R function. e.g. the command is still not complete. For file names however. For the names of R objects or commands.start() opens a browser window that gives access to the full range of documentation for syntax.g.
1. with the semicolon (.4
>
Further Notational Details
As noted earlier. Typing q on its own. but omit the + that will appear on the commands window if the command is typed in as we show it. the continuation prompt is
+
In these notes. the Microsoft Windows conventions apply. type in
> help(plot)
The two search functions help.3.
1.
6
Multiple commands may appear on the one line. Try it!
1. the command line prompt is R commands (expressions) are typed following this prompt6. and case does not distinguish file names. On Unix systems letters that have a different case are treated as different. type:
> help()
In R for Windows. Experimentation often helps clarify the precise action of an R function.. used when.search("matrix") (This lists all functions whose help pages have a title or alias in which the text string “matrix” appears.5
On-line Help
To get a help window (under R for Windows) with a list of help topics. Place this file in the working directory and.7 1. There is also a continuation prompt.3.RData")
Files that are mentioned in these notes. case is significant. and that are not supplied with R (e. and the size and colour of the plotting symbol. Note: Recall that. available from the author's web page. following a carriage return. Details will be given in Section 3.4 Options for plot() and allied functions
The function plot() and related functions accept parameters that control the plotting symbol. Examples of their use are:
> help. displays the text of the function on the screen. without the parentheses. type:
> attach("usingR. an alternative is to click on the help menu item.) The function help. This is because q is a function. we often continue commands over more than one line. in order to quit from the R session we had to type q(). packages and functions. from within the R session. Thus Austpop is different from austpop. plot().

To save keystrokes. 2. These are the data. iv. The following ten observations.5 1978 14. iii Use the hist() command to plot a histogram of the snow cover values. taken during the years 1970-79.3.7
Exercises
1.)
Temperature Erosion (F) 53 57 63 70 70 75 3 1 1 1 1 0 Blowby 2 0 0 0 0 2 Total incidents 5 1 1 1 1 1 incidents incidents
Enter these data into a data frame.5 1971 12. Plot total incidents against temperature. Distinguish between the attaching of image files and the attaching of data frames.2
i. (Refer back to Section 1. the difference being that with attach() datasets are loaded into memory only when required for use. Input the following data.0 1972 14.8
Users can also load (use load()) or attach (use attach()) specific files. from the DAAG package. In the data frame elasticband from section 1. Enter the data into R.
1. These have a similar effect. (Snow cover is in millions of square kilometers):
year snow. 3. erosion.1.1).9 1973 10. are on October snow cover for Eurasia. [Section 1.cover 1970 6.0 1974 10. with (for example) column names temperature.cover versus year.5 1979 9. The attaching of data frames will be discussed later in these notes. blowby and total.9 1976 21. on damage that had occurred in space shuttle launches prior to the disastrous launch of Jan 28 1986. for 6 launches out of 24.3. that were included in the pre-launch charts that were used in deciding whether to proceed with the launch. Repeat ii and iii after taking logarithms of snow cover. enter the successive years as 1970:1979] ii. plot distance against stretch. (Data for the 23 launches where information is available is in the data set orings.
.3.1 showed one way to do this. Plot snow.9 1977 12.7 1975 7.

then take sin() [1] 0.60 # Assumes hills. compounded annually [1] 435. We will discuss graphical summaries in the next section.95 1st Qu. As a first example.500 Median: 6.6 minutes. for five years
2.2 R will provide numerical or graphical summaries of data
A special class of object. There is one column for each variable.5000000 0.Rdata is in the working directory
We may for example require information on ranges of variables.: 4. Thus the range of distances (first column) is from 2 miles to 28 miles.0000000 # at 7. An Overview of R 2. where the rows are observations and the columns are variables.60. For now. climb.1 R may be used as a calculator.
7
There are also versions in the DAAG package and in the Venables and Ripley MASS package.1 The Uses of R
2.529 3rd Qu.000 climb Min.: 725 Median:1000 Mean:1815 3rd Qu. appears on subsequent lines
> 2+2 [1] 4 > sqrt(10) [1] 3. think of data frames as matrices. radius is 6378 km [1] 40074. consider the data frame hills that accompanies these notes7.6293 > > pi # R knows about pi [1] 3. Typing in summary(hills)gives summary information on these variables.000 1st Qu. Data frames are central to the way that all the more recent R routines process data.1000 # Interest on $1000.a.162278 > 2*3*4*5 [1] 120 > 1000*(1+0.: 2. while the range of times (third column) is from 15.00 Median: 39.000 Max.: 300 1st Qu. stores rectangular arrays in which the columns may be vectors of numbers or factors or text strings.
R evaluates and prints out the result of any expression that one types in at the command line in the console window.141593 > 2*pi*6378 #Circumference of Earth at Equator.: 8. Expressions are typed following the prompt (>) on the screen. if any. The result.95 (minutes) to 204. and time.16 > sin(c(30. in km.:7500 time Min. thus:
> load("hills.
.:28.1.:2200 Max.: 28.75 Mean: 57.:204.1.62 Max.075)^5 . with the names distance.90)*pi/180) # Convert angles to radians.Rdata") > summary(hills) distance Min.: 68. called a data frame.8660254 1.88 3rd Qu.5% p.9
2.000 Mean: 7.: 15. This has three columns (variables).

Why has this happened? Straight Line Regression: Here is a straight line regression calculation. for specifying tick marks and tick labels.
.805 0. The data are stored in the data frame elasticband (DAAG package).652 0. Correlation: We calculate the correlation matrix for the hills data:
> options(digits=3) > cor(hills) distance climb distance climb time time 1.000
Suppose we wish to calculate logarithms.724 0. and so on.890 0.10 2. and time as log(time). shown in Figure 5.700 0. For this. and between time and climb.920 0.652 1. A helpful graphical summary for the hills data frame is the scatterplot matrix.3 R has extensive graphical abilities
The main R graphics function is plot(). have reduced. In addition to plot() there are functions for adding points and lines to existing graphs. and then calculate correlations. We can do all this in one step.000
Unfortunately R was not clever enough to relabel distance as log(distance). climb as log(climb). There are various other alternative helpful forms of graphical summary. for placing text at specified positions. for labelling axes.4 R will handle a variety of specific analyses
The examples that will be given are correlation and regression.920 0.89 0. thus:
> cor(log(hills)) distance climb distance climb time time 1.805 1. type:
> pairs(hills)
Figure 5: Scatterplot matrix for the Scottish hill race data
2.1. Notice that the correlations between time and distance.00 0.1. The formula that is supplied to the lm() command asks for the regression of distance travelled by the elastic band (distance) on the amount by which it is stretched (stretch). The variable names are the names of columns in that data frame.000 0.70 1.000 0.000 0.724 1.

RData in the working directory. Careful housekeeping may be needed to distinguish between objects that are
8
Type in help(ls) and help(grep) to get details. for use in the next session in the same workspace. etc. An alternative to ls() is objects().data=elasticband) Call: lm(formula = distance ~ stretch. data = elasticband) Coefficients: (Intercept) -63.g.25:30 > fahrenheit <.RData")
Image files. file="tempscales.image() # Save contents of workspace. by default in the file . it makes sense to save the contents of the working directory from time to time.9.554
More complete information is available by typing
> summary(lm(distance~stretch. Typing the name of an object causes the printing of its contents.11
> plot(distance ~ stretch. They can all be operated on as data.RData") ls(pos=2) # Check the contents of the file that has been attached
The parameter pos gives the position on the search list.571 stretch 4. from the working directory or (with the path specified) from another directory.0 78. Try typing q.image(file="archive. or collections of objects into a named image file.0
2. e.2 R Objects
All R entities.5 R is an Interactive Programming Language
We calculate the Fahrenheit temperatures that correspond to Celsius temperatures 25. The pattern matching conventions are those used for grep(). pch=16) > elastic.RData # Save into the file archive.6 82. which is modelled on the Unix grep command. including functions and data structures. in Section 2.8 80.2 86. exist as objects.data=elasticband. In both cases there is provision to specify a particular pattern.data=elasticband))
Try it!
2. ….RData save. Fahrenheit=fahrenheit) > print(conversion) Celsius Fahrenheit 1 2 3 4 5 6 25 26 27 28 29 30 77.
. can be attached.lm(distance~stretch. Some possibilities are:
save. 30:
> celsius <. fahrenheit. R offers the option of saving the workspace image. This allows the retention. In a long session. It is also possible to save individual objects. any objects that were created in the current session. For example
attach("tempscales.9/5*celsius+32 > conversion <.RData")
save(celsius.data=elasticband) > lm(distance ~stretch.1. starting with the letter `p’8. 26. (The search list is discussed later in this chapter. mean.data.lm <. into the file .4 84.) Important: On quitting.frame(Celsius=celsius. Type in ls() to see the names of all objects in your workspace. thus making objects in the file available on request.

FALSE.
*92.j+answer} > answer [1] 173
The calculation iteratively builds up the object answer. Then j=51.4 Vectors
Examples of vectors are
c(2.8 [1] 27..12
to be kept and objects that will not be used again.0 82. Initially.
2.3. i.e.1) 3:10 # The numbers 3.0 84.51.”Newcastle”. Before typing q() to quit.TRUE.51.0 > for (j in c(31.91)) [1] 173
Skilled R users have limited recourse to loops. Saving the workspace image will then save everything remains. numeric or character11.0 78. and answer is assigned the value 91 + 81 = 173. Then the procedure ends. 10 c(TRUE. There are often. better alternatives. and answer is assigned the value 51 + 31 = 82.51.”Sydney”.FALSE. or a sequence of statements that are enclosed within braces
..91).1 More on looping
Here is a long-winded way to sum the three numbers 31. a vector with elements of mode character). Finally.3 Looping
A simple example of a for loop is10
for (i in 1:10) print(i)
Here is another example of a for loop.2.1.
9
Asterisks (*) identify sections that are more technical and might be omitted at a first reading Other looping constructs are: repeat <expression> ## break must appear somewhere inside the loop
10
while (x>0) <expression> Here <expression> is an R statement. using the successive values of j listed in the vector (31.”Darwin”)
Vectors may have mode logical.91)){answer <. j=31. as in this and earlier examples.e. and answer is assigned the value 31 + 0 = 31.TRUE. to do in a complicated way what we did very simply in section 2.FALSE.e. . The first two vectors above are numeric.0 80. 51 and 91:
> answer <.5:
> # Celsius to Fahrenheit > for (celsius in 25:30) + print(c(celsius.4 [1] 29. the third is logical (i. a vector with elements of mode logical).FALSE) c(”Canberra”. The workspace image will be automatically loaded upon starting another session in that directory.5. 4. and the fourth is a string vector (i.2 [1] 30 86
2.3. 9/5*celsius + 32)) [1] 25 77 [1] 26. and the contents of answer can be examined by typing in answer and pressing the Enter key.6 [1] 28. use rm() to remove objects that are no longer required. j=91.7. There is a more straightforward way to do this calculation:
> sum(c(31.

c(2. >=. Thus suppose we want to extract values of x that are greater than 10.5.c(10. Jeff=183)[c("John". and != tests for inequality. can be included as an element of a vector.7. 12 # Extract elements (rows) 2 and 4
One can use negative numbers to omit elements:
> x <.15.
12
A third more subtle method is available when vectors have named elements. Use the function class() to determine the class of an object. == tests for equality.
2.4)] [1] 11 15 # Assign to x the values 3. 8. NA. 1. <=."Jeff")] John Jeff 185 183
. Specify the numbers of the elements that are to be extracted. 0. 3. the meaning is: “Join these numbers together in to a vector.
2.
> x>10 # This generates a vector of logical (T or F) [1] F T F T T > x[x>10] [1] 11 15 12
Arithmetic relations that may be used in the extraction of subsets of vectors are <.1 Joining (concatenating) vectors
The c in c(2. ==.11. NA)
11
It will.2 Subsets of Vectors
There are two common ways to extract subsets of vectors12.12) > y [1] 10 15 12 > z <.c(3.c(x.
> x <.2.e.15.8.13
The missing value symbol. One can then use a vector of names to extract the elements.g. 7. >. thus:
> c(Andreas=178.12) > x[-c(2.4. 3. which is NA. which we then concatenate to form a vector z:
> x <.11. 1) above is an acronym for “concatenate”. plot() and summary().
2.8. Existing vectors may be included among the elements that are to be concatenated. John=185.c(1. Specify a vector of logical values.15.4. be important to know the “class” of such objects.12) > x[c(2.3. i. The elements that are extracted are those for which the logical value is T. Set
y <. 5.3 The Use of NA in Vector Subscripts
Note that any arithmetic operation or relation that involves NA generates an NA. This determines how the method used by such generic functions as print().c(3.4. and !=.3)] [1] 3 15 12
2. y) > z [1] 2 3 5 2 7 1 10 15 12
The concatenate function c() may also be used to join lists.1) > x [1] 2 3 5 2 7 1 > y <. 11. e. 15. The first four compare magnitudes. later in these notes. In the following we form vectors x and y.

Incorrect spelling of the level names will generate an error message.c(rep(“female”. we can create a vector of strings that that holds the values thus:
gender <.0
2. and there is no assignment. use
y[is. They are crucial in the representation of categorical effects in model and graphics formulae. The class attribute of a factor has. The reason is that all elements of y==NA evaluate to NA. In most cases where the context seems to demand a character string.factor(gender.
.691). 3. created as above [1] "female" "male"
The order of the levels in a factor determines the order in which the levels appear in graphs that use this information.692))
(The usage is that rep(“female”.relevel(gender.691). use
gender <. levels=c(“Male”. the space required for storage is reduced. Hence:
> levels(gender) # Assumes gender is a factor. Consider a survey that has data on 691 females and 692 males. The values “female” and “male” are the levels of the factor."male" rows now hold missing values table(gender) rm(gender) # Remove gender
13
The attributes() function makes it possible to inspect attributes. the levels are in alphanumeric order. not surprisingly. rep(“male”. the value “factor”. If the first 691 are females and the next 692 males. “female”)) # Erroneous .) We can change the vector to a factor.4 Factors
A factor is stored internally as a numeric vector with values 1.factor(gender. levels=c(“male”. levels=c(“male”.4. ref=“male”)
An alternative is
gender <.na(y)] <. 2. or later when one wishes to change the order of levels in an existing factor. For example
attributes(factor(1:3))
The function levels() gives a better way to inspect factor levels. To cause “male” to come before “female”. and similarly for the creation of 692 copies of “male”. By default. “female”))
This last syntax is available both when the factor is first created. Factors provide a compact way to store character strings. the 1 is translated into “female” and the 2 into “male”. 691) creates 691 copies of the character string “female”. An attributes table gives the ‘level’ for each integer value13.factor(gender. rep(“male”.692))) table(gender) gender <. Try
gender <. “female”)) table(gender) gender <.factor(gender)
Internally the factor gender is stored as 691 1’s. so that “female” precedes “male”.0 leaves y unchanged. and in tables. by entering:
gender <. where k is the number of levels.factor(c(rep(“female”. This does not select an element of y. k. followed by 692 2’s.14
Be warned that y[y==NA] <. It has stored with it the table:
1
2
female male
Once stored as a factor. To replace all NAs by 0.

summary)) Compact.g.frame(). which gives girth. in which different columns may have different modes.1 Data frames as lists
A data frame is a list15 of column vectors.5.. all of equal length.summary[. # Take the object that is stored
2. vectors. They may be any mixture of scalars. Here it is:
> Cars93.names(Cars93.
Here is summary information on this data frame
14 15
Also legal is Cars93. etc.4] type <.summary)) are Min.summary[.cars. Just as with any other list. In general forms of list. 4].summary[. As noted above.summary[.summary[[4]] or Cars93. e.passengers Max.summary[1. All elements of any column must however have the same mode. This gives a data frame with the single column Type.Cars93. use Cars93.Cars93. available both with read. an alternative is
as. in the datasets package. i. then storing it in the vector type.3 Built-in data sets
We will often use data sets that accompany one of the R packages. The parameter setting stringsAsFactors=TRUE.2 Inclusion of character string vectors in data frames
When data are input using read.summary. A data frame is a generalisation of a matrix. Max.frame() function is used to create data frames. functions.
2. The use of matrix-like subscripting. and abbrev.4] or Cars93. .”abbrev”] type <.4] to extract the column vector. all numeric or all factor. One such data frame.5. or when the data. Columns can be vectors of any mode. Among the data sets in the DAAG package is Cars93. Cars93. or all character.Cars93. The column abbrev could equally well be stored as a factor.summary. will if needed ensure that character strings are input without such conversion.5 Data Frames
Data frames are fundamental to the use of the R modelling and graphics functions.
2.passengers No.passengers (i.of. The first three columns have mode numeric.5. For read. vectors of character strings are by default turned into factors. . is trees.
type <.15
2.table() and with data. No. elements that are of arbitrary type.is=TRUE. subscripting extracts a list. usually stored as data frames.summary Min. the minimum number of passengers for cars in this category).summary.of.summary[2]. Large.e.table().
. which is the fourth column vector of Cars93. created from information in the Cars93 data set in the Venables and Ripley MASS package.cars abbrev Compact Large Midsize Small Sporty Van 4 6 4 4 2 7 6 6 6 5 4 8 16 11 22 21 14 9 C L M Sm Sp V
The data frame has row labels (access with row.passengers. takes advantage of the rectangular structure of data frames.Cars93. Any of the following14 will pick out the fourth column of the data frame Cars93. The column names (access with names(Cars93.summary$abbrev type <. Thus Cars93. height and volume for 31 Black Cherry Trees.summary[4] is a data frame with a single column. . and the fourth has mode character.e.table().summary[[4]] # in the fourth list element.

There is an example that shows how to use it to count the number of missing values in each column of data.rm.17
The functions mean() and range(). fraseri Acmena smithii B. called databases. Specify x <.] 4 3 8 105 2 135 0. This can be changed in the course of a session.NA. then determine the range
One can specify na.rm=TRUE as a third argument to the function sapply. For example
> range(rainforest$branch.GlobalEnv" "package:methods" "package:stats"
.1 Numbers of NAs in subgroups of the data
The following gives information on the number of NAs in subgroups of the data:
> table(rainforest$species.e. For example:
> sapply(rainforest[. !is.7 Making Tables
table() makes a table of counts.8) > x <.
2.factor(x) > x [1] 1 Levels: [1] 1 Levels: 5 5 NA 8 1 5 8 NA 8 1 5 8 NA
> factor(x. To get a full list of these directories.5. na.factor(x. range. annoyingly. The action needed to get NAs tabulated under a separate NA category depends.8 The Search List
R has a search list where it looks for objects. and a number of other functions.exclude=NULL)
2. This argument is then automatically passed to the function that is specified in the second argument position.-7].exclude=NULL)
> x <.rm=TRUE) dbh wood bark root rootsk branch [1. barley$site) Grand Rapids Duluth University Farm Morris Crookston Waseca 1932 10 1931 10 10 10 10 10 10 10 10 10 10 10
WARNING: NAs are by default ignored.7. type:
> search() [1] ".3 24.na(rainforest$branch)) FALSE TRUE Acacia mabellae C.c(1. myrtifolia 6 0 15 1 10 12 11 10
Thus for Acacia mabellae there are 6 NAs for the variable branch (i.
2. If the vector is a factor then it is necessary to generate a new factor that includes “NA” as a level. Specify one vector of values (often a factor) for each table margin that is
required. out of a total of 16 data values. on whether or not the vector is a factor. take the parameter na. specify exclude=NULL. For example:
> library(lattice) # The data frame barley accompanies lattice > table(barley$year. number of branches over 2cm in diameter).rm=TRUE) [1] 4 120 # Omit NAs. na.0 4 120 56 1530
Chapter 8 has further details on the use of sapply(). If the vector is not a factor.] [2.

to. its columns can be accessed by giving their names. BUCHANAN. mean(Bodywt))
2.0 62. Note also the function with().2 A Plotting function
The data set florida has the votes in the 2000 election for the various US Presidential candidates. the data frame becomes a database. Remember to detach an attached data frame. in miles
The function will do the conversion for several distances all at once. detach(primates).1 An Approximate Miles to Kilometers Conversion
miles. In technical terms. which attaches the data frame that is given as its first argument for the duration of the calculation that is specified by its second argument.GlobalEnv" [4] "package:stats" [7] "package:utils" [10] "package:base" "package:MASS" "package:graphics" "package:datasets" "package:methods" "package:grDevices" "Autoloads"
Use of attach() likewise extends the search list. This is however an uncommon construction
.9. xlab="Bush". specify:
> miles. which is searched as required for objects that the user may specify.g.9 Functions in R
We give two simple examples of R functions.300)) [1] 160 320 480
2. without quotes) or image (.RData) files (the file name is placed in quotes).km <. ylab="Buchanan") detach(florida) # In S-PLUS.
attach(florida) plot(BUSH. For example:
> av <. The following demonstrates the attaching of the data frame primates:
> names(primates) [1] "Bodywt" > Bodywt Error: Object "Bodywt" not found > attach(primates) > Bodywt [1] 10.
2. specify detach("florida")
18
Alternatively a return value may be given using an explicit return() statement. This function can be used to attach data frames or lists (use the name.km(175) [1] 280 # Approximate distance to Sydney.. Use the function thus
> miles.9.to.200.km(c(100.2 # R will now know where to find Bodywt "Brainwt"
Once the data frame primates has been attached.function(miles)miles*8/5
The return value is the value of the final (and in this instance only) expression that appears in the function body18.to.8 52. e. county by county in the state of Florida. without further reference to the name of the data frame.0 6.with(primates.
> library(MASS) > search() [1] ".0 207. The following plots the vote for Buchanan against the vote for Bush.18
[4] "package:graphics" [7] "package:datasets" "package:grDevices" "package:utils" "Autoloads" "package:base"
Notice that the loading of a new package extends the search list. 200 and 300 miles to distances in kilometers. when it is no longer required. To convert a vector of the three distances 100.

in the US 2000 Presidential election. additional to those covered above. predict the result.
2.2) Unexpected consequences (e.g.11 Exercises
1. how would you expect prod() to work? Try it! 3. Look up the help for the function prod(). conversion of columns of numeric data into factors) from errors in data (7. i. Figure 6 shows the graph produced by plot.
Figure 6: Election night count of votes received. line=1.j+answer }
c) answer <.75. parameter settings are left at their defaults.florida(xvar=”GORE”.e. xlab=xvar.19
Here is a function that makes it possible to plot the figures for any pair of candidates. “Votes in Florida.yvar] plot(x. y. Alternatively. the function allows.10 More Detailed Information
Chapters 7 and 8 have a more detailed coverage of the topics in this chapter. at this point.0
for (j in 3:5){ answer <. to glance through chapters 7 and 8. Then do the computation: a) answer <. Add up all the numbers from 1 to 100 in two different ways: using for and using sum. It may pay. As well as plot.ylab=yvar) mtext(side=3.xvar] y <.florida[. that may be important for relatively elementary uses of R include: The entry of patterned data (7.3) The handling of missing values in subscripts when vectors are assigned (7.10
for (j in 3:5){ answer <. by county.function(xvar=”BUSH”.
plot. and use prod() to do the calculation in 1(c) above.
plot.florida(). in \nthe 2000 US Presidential election”) }
Note that the function body is enclosed in braces ({ }).10
for (j in 3:5){ answer <. Remember also to use R’s help pages and functions.g. yvar=”BUCHANAN”){ x <.florida().j*answer }
2. Now apply the function to the sequence 1:100. by county. For each of the following code sequences. yvar=”NADER”)
2. Topics from chapter 7. e.florida(yvar=”NADER”) # yvar=”NADER” over-rides the default plot.1.j+answer }
b) answer<.florida[. What is its action?
.4.1).florida <.

Which columns are ordered factors? [Use is. 20 find the corresponding volumes and print the results out in a table. Multiply all the numbers from 1 to 50 in two different ways: using for and using prod.
. 6.1. 5. Use the technique of section 2.factor to each column of the supplied data frame tinting. Use sapply() to apply the function is. determine the levels. The volume of a sphere of radius r is given by 4 r3/3. …. For spheres having radii 3.5 to construct a data frame with columns radius and volume.20
4. For each of the columns that are identified as factors. 5. 4.ordered()].

enter
demo(graphics)
Press the Enter key to move to each new graph. Try
plot((0:20)*pi/10. Explicitly setting type = "p" causes either function to plot points.92. Is it obvious that these points lie on a sine curve? How can one make it obvious? (Place the cursor over the lower border of the graph sheet. text(). data=austpop. and for lines() is type = "l".
par(cex=1. points(). The mtext() function places text in one of the margins. type="l") plot(ACT ~ Year. specify detach(“austpop”)
3. ACT). sin((1:30)*0. Drag the border in towards the top border. type = "l" gives lines.
3. data=austpop. type="l") detach(austpop) # Fit smooth curve through points # In S-PLUS. The axis() function gives fine control over axis ticks and labels. until it becomes a doublesided arror. The default setting for points() is type = "p". lines().1. axis(). When it is necessary to change parameter settings for a subsequent plot. For example. for each numeric variable. form a suite that plots points.21
3. The text() function adds text at specified locations. the par() function does this. type="b")
The points() function adds points to a plot. making the graph sheet short and wide.
19
Actually these functions differ only in the default setting for the parameter type.
.2 Fine control – Parameter settings
The default settings of parameters. a normal probability plot. Plotting
The functions plot(). Here is a further possibility
attach(austpop) plot(spline(Year.1 Plot methods for other classes of object
The plot function is a generic function that has special methods for “plotting” various different classes of object.) Here are two further examples.92))
Comment on the appearance that these graphs present. To see some of the possibilities that R offers. For example.
attach(elasticband) # R now knows where to find distance & stretch plot(distance ~ stretch) plot(ACT ~ Year. Try
plot(hills) # Has the same effect as pairs(hills)
3. identify() etc. y) # Use a formula to specify the graph #
Obviously x and y must be the same length. sin((0:20)*pi/10)) plot((1:30)*0. The lines() function adds lines to a plot19. plotting a data frame gives. such as character size. are often adequate. lines and text. mtext(). Plotting the lm object that is created by the use of the lm() modelling function gives diagnostic and other information that is intended to help in the interpretation of regression results.25) # character expansion
increases the text and plot symbol size 25% above the default.1 plot () and allied functions
The following both plot y against x:
plot(y ~ x) plot(x.

2 The shape of the graph sheet
Often it is desirable to exercise control over the shape of the graph page.2.25) # mex=1.exit(par(oldpar))
Type in help(par) to get details of all the parameter settings that are available with par().names(primates) <.2 115 406 1320 179 440
Observe that the row names store labels for each row20.25.0 62.8 52. and pulling. They can be assigned directly.graph() or x11() that set up the Windows screen take the parameters width (in inches).log(brain)) detach(Animals) par(mfrow=c(1. enter par(oldpar). To restore the original parameter settings at some later time.1 Multiple plots on the one page
The parameter mfrow can be used to configure the graphics sheet so that subsequent plots appear row by row. it is often useful to store existing settings. on the one page.par(cex=1. For this.25 expands the margin by 25%
This stores the existing settings in oldpar. e.g. which must be attached
3.g."Human". holding down the left mouse button. in a two by two layout:
par(mfrow=c(2. lines and text
Here is a simple example that shows how to use the function text() to add text labels to the points on a plot. Here is an example:
attach(elasticband) oldpar <."Chimp")
. mex=1.
row. specify
oldpar <.g. use mfcol instead. For a column by column layout. (brain)^0.c("Potar monkey". pch=1) # Restore to 1 figure per page # This dataset is in the MASS package. It is the relative sizes of these parameters that matter for screen display or for incorporation into Word and similar programs.22
On the first use of par() to make changes to the current device. The setting of pointsize (default =12) determines character heights.0 207.par(cex=1.25) on.par(cex=1. e.5) plot(distance ~ stretch) par(oldpar) detach(elasticband) # Restores the earlier settings
Inside a function specify. so that the individual plots are rectangular rather than square.3 Adding points. "Rhesus monkey".
20
Row names can be created in several different ways.
3. The R for Windows functions win.0 6.1) plot(log(body).
3. Graphs can be enlarged or shrunk by pointing at one corner. brain) plot(sqrt(body).
oldpar <. then changes parameters (here cex and mex) as requested.2). pch=16) attach(Animals) plot(body. height (in inches) and pointsize (in 1/72 of an inch). so that they can be restored later. one after the other in a rectangular layout.1). In the example below we present four different transformations of the primates data.1. e.
> primates Bodywt Brainwt Potar monkey Gorilla Human Rhesus monkey Chimp 10. sqrt(brain)) plot((body)^0."Gorilla".2.

The parameter pch controls the choice of plotting symbol. pos=4) # Labels with symbol number
When using read. colour and choice of plotting symbol
For plot() and points() the parameter cex (“character expansion”) controls the size. 7).
Figure 7: Plot of the primates data. 7). Figure 7A would be adequate for identifying points.
.1350)) # Specify xlim so that there is room for the labels text(x=Bodywt. added to the plot that results from the above three statements. xlim=c(0. a column that holds the row names.names(primates). We now show how to improve it. Try
plot(1. col=1:7. 1.5. pos=2)
3. rep(2.names(primates). ylim=c(1. pch=0:6) text(1:7.5. the parameter row. pch=(0:6)+14) # Plot symbols 7 to 13 # Label with symbol number # Plot symbols 14 to 20
text((1:7). specify
text(x=Bodywt.5.rep(3. col=1:7) # Do not plot points
The following. pch=16. paste((0:6)+7). paste((0:6)+14). y=Brainwt. ylab="") box() points(1:7.rep(2. with labels on points. plot(Bodywt. cex=1:7. This helps make the points stand out against the labelling. cex=1:7. while col (“colour”) controls the colour of the plotting symbol. by number or name. labels=row. xlim=c(1.names is available to specify. rep(2. but is not a presentation quality graph. labels=row. 7. It uses the parameter setting pos=4 to move the labelling to the right of the points. The parameters cex and col may be used in a similar way with text(). pch=(0:6)+7) text((1:7).table() to input data. labels=row. Figure 7B uses the xlab (x-axis) and ylab (y-axis) parameters to specify meaningful axis titles. y=Brainwt. pos=4) points(1:7. xlab="Body weight (kg)". demonstrates other choices of pch.280).7).rep(2.7). type="n".7).5). ylab="Brain weight (g)".23
attach(primates) # Needed if primates is not already attached. labels=paste(0:6).5. rep(4.names(primates). It sets pch=16 to make the plot character a heavy black dot. xlab="".3.1 Size. Figure 8B improves on Figure 8A. pos=4) detach(primates)
To place the text to the left of the points.5).
points(1:7. pos=4) # pos=4 positions text to the right of the point
Figure 7A shows the result. ylim=c(0. Brainwt) text(x=Bodywt.7). y=Brainwt. Here is the R code for Figure 7B:
plot(x=Bodywt. axes=F. y=Brainwt.75.

4=right). cex=3) text(10. rainchars. unless the setting of the parameter n is reached first."Y". Draw the graph first. 3(top) and 4."O".5. col=palette()[1:6].. colours and sizes A variety of color palettes are available. 0.colors".function(){ plot(1.1:9) # Split “cm. cex=2.1 identify()
This function requires specification of a vector x. paste(1:6).colors” into its 9 characters text(1:9.4 Identification and Location on the Figure Region
Two functions are available for this purpose. A click with the right mouse button signifies that the identification or location task is complete. axes=F. rep(0. rep(2. The sides are numbered 1(x-
axis).5. "cm.14)."B". "rainbow(7)". 2. 1:9. 1.5.5) text(10. .
3."I". then call one or other of these functions. rep(1. 1. and clicks the left mouse button. adj=0) cmtxt <.5. xlim=c(0.24
Note the use of pos=4 to position the text to the right of the point (1=below. enter
view. "Default palette". a vector y.
3. text."V") text(1:7.substring("cm. For identify() the default setting of n is the number of data points. line.colors(9).5.5) text(10.3. identify() labels points. xlab="".ylab="") text(1:6. 2(y-axis). adj=0) }
To run the function.9). and clicks the left mouse button.colours <. cmtxt. while for locator() the default setting is n = 500.c("R". adj=0) rainchars <.3). county by county in
. and a vector of text strings that are available for use a labels. Here is a function that displays some of the possibilities:
view.colours()
3. type="n". ylim=c(0. locator() prints out the co-ordinates of points.4.colors(9)". col=cm.2 Adding Text in the Margin
mtext(side. One positions the cursor near the point that is to be identified. cex=2. col=rainbow(7).5.7).) adds text in the margin of the current plot. The data set florida has the votes for the various Presidential candidates. 3=top. One positions the cursor at the location for which coordinates are required."G". 2=left. Here (Figure 8) is the plot:
Figure 8: Different plot symbols.6).

boxplots and normal probability plots.5.2 locator()
Left click at the locations whose coordinates are required
attach(florida) locator() detach(florida) # if not already attached plot(BUSH.") detach(possum)
Density plots.5. BUCHANAN.5 Plots that show the distribution of data values
We discuss histograms. in fact. .
3. density plots. When labelling is terminated (click with the right mouse button). A histogarm is. County) detach(florida)
Click to the left or right. 77. so that a density curve can be readily superimposed.
attach(florida) plot(BUSH.1 Histograms and density plots
The shapes of histograms depend on the placement of the breaks. BUCHANAN. ylab=”Buchanan”)
The function can be used to mark new points (specify type=”p”) or lines (specify type=”l”) or both points and lines (specify type=”b”). Here is the code used to plot the histogram in Figure 10A. again. as Figure 9 illustrates:
Figure 9: Panel A shows a histogram with a frequency scale. now that they are available. ylab=”Buchanan”) identify(BUSH.. 22). breaks = 72.sex == "f" hist(totlngth[here].
3. superimposed.25
the state of Florida. xlab=”Bush”. BUCHANAN. in order.
. then invoking identify() so that we can label selected points on the plot. are often a preferred alternative to a histogram. xlab="Total length".
attach(possum) here <.5. main ="A: Breaks at 72. so that the histogram gives a rather different impression of the distribution of the data. xlab=”Bush”. Panel C has a different choice of breakpoints. We plot the vote for Buchanan against the vote for Bush. The density curve is. depending on the preferred positioning of the label. and slightly above or below a point. Panel B is drawn with a density scale.4. a very crude form of density plot.
3..5 + (0:5) * 5. the row numbers of the observations that have been labelled are printed on the screen. ylim = c(0.

xlab="Total length"..density(totlngth[here]) xlim <. main="") lines(dens) detach(possum)
3. does affect the appearance of the plot.e.
3. In order to calibrate the eye to recognise plots that indicate nonnormal variation. The code for Figure 9B is:
attach(possum) here <. The main effect is to make it more or less smooth. The code for the boxplot is
attach(possum) boxplot(totlngth[here]. so that the histogram gives a rather different impression of the distribution of the data. i.range(dens$x) ylim <. The choice of width and type of window.3 Boxplots
We now (Figure 10) make a boxplot of the above data:
Figure 10: Boxplot of female possum lengths.26
Density plots do not depend on a choice of breakpoints. probability = TRUE. The only difference between Figure 9C and Figure 9B is that a different choice of breakpoints is used for the histogram.range(dens$y) hist(totlngth[here].5. it is helpful to do several normal probability plots for random samples of the relevant size from a normal distribution. ylim = ylim. The following will give a density plot:
attach(possum) plot(density(totlngth[here]).5 + (0:5) * 5.sex == "f" dens <. side=1) detach(possum)
Figure 12 adds information that is intended to assist in the interpretation of boxplots. horizontal=TRUE) rug(totlngth[here]. with additional labelling information. height=6) # This is a better shape for this plot
. The points of this plot will lie approximately
on a straight line if the distribution is Normal.4 Normal probability plots
qqnorm(y) gives a normal probability plot of the elements of y.type="l") detach(possum)
In Figure 9B the y-axis for the histogram is labelled so that the area of a rectangle is the density for that rectangle.
x11(width=8. This gives a scale that is appropriate for superimposing a density curve estimate. controlling the nature and amount of smoothing. xlim = xlim.5. the frequency for the rectangle is divided by the width of the rectangle. breaks = 72.

sex == "f" par(mfrow=c(2. sport and body-size dependency of hematology in highly trained athletes.off() # A 2 by 4 layout of plots
Figure 11 shows the plots. R.B. Otherwise the points for the female possum lengths are as close to a straight line as in many of the plots for random normal data. R.4)) y <. If data are from a normal distribution then points should fall. The idea is an important one. ylab = “% Body fat”) panel. main="Possums") for(i in 1:7)qqnorm(rnorm(43). ylab="Simulated lengths". It is a way to train the eye. xlab = “Height”.smooth() plots points.xlab="". and Cunningham. The plot in the top left hand corner shows the 43 lengths of female possums.smooth(ht[here]. type dev. 1991: Sex. By default. Medicine and Science in Sports and Exercise 23: 788-794.6 Other Useful Plotting Functions
For the functions demonstrated here. examine a number of randomly generated samples of the same size from a normal distribution. For example:
attach(ais) here<. we use data on the heights of 100 female athletes21.
3.sex=="f" plot(pcBfat[here]~ht[here].1 Scatterplot smoothing
panel.D. then adds a smooth curve through the points.
3. rnorm() generates random samples from a distribution with mean 0 and standard deviation 1.pcBfat[here]) detach(ais)
21
Data relate to the paper: Telford.6. main="Simulated") detach(possum) # Before continuing.
Figure 11: Normal probability plots.
.xlab="". along a line. The other plots are for independent normal random samples of size 43. ylab="Length". In order to judge whether data are normally distributed. There is one unusually small value.27
attach(possum) here <. approximately.totlngth[here] qqnorm(y.

25. xlab = “Height”. The following is better.1.6. The parameters may be an intercept and slope. or a vector that holds the intercept and slope.5 Dotcharts
These can be a good alternative to barcharts.0)) boxplot(ht[here].2 Adding lines to plots
Use the function abline() for this. ylab = "Height". cex=0. 0.
3. vertical bars showing the distribution of values of x. but shrinks the sizes of the points so that they almost disappear:
dotchart(islands.1. side = 1) par(oldpar) detach(ais)
The parameter boxwex controls the width of the boxplot.1.6.1.6. and there is substantial overlap.
attach(ais) here<.3 Rugplots
By default rug(x) adds. Alternatively it is possible to draw a horizontal line (h = <height>). ylab = “% Body fat”) abline(lm(pcBfat[here] ~ ht[here])) detach(ais)
3.35.5)
. boxwex = 0. It can however be particularly useful for showing the actual values along the side of a boxplot.1). Here is the code
attach(ais) here <.
3.sex == "f" oldpar <. They have a much higher information to ink ratio! Try
dotchart(islands) # vector of named numeric values. along the x-axis of the current plot.75.par(mar=c(4. with a rugplot added on the y-axis.6.1. or an lm object.3 demonstrated the use of the pairs() function. or a vertical line (v = <ordinate>). horizontal=TRUE) rug(ht[here]. mgp=c(2.sex=="f" plot(pcBfat[here] ~ ht[here].
Figure 12: Distribution of heights of female athletes.28 3. in the datasets base package
Unfortunately there are many names. 1.1. Figure 12 shows a boxplot of the distribution of height of female athletes.4 Scatterplot matrices
Section 2. The bars on the left plot show actual data values.

or whatever. in feet. Use colour or different plotting symbols to distinguish different groups. although what will appear on the plot is Area = pi*r^2. Notice that in expression(Area == pi*r^2). Take care to use colours that contrast. It is pointless to present information on a source of error that is of little or no interest. with a minimum waste of ink. either or both of xlab and ylab can be an expression. Explain what source of `error(s)’ is represented. Figure 13 was produced with
x <. when there is little or no overlap.
. large enough to stand out relative to other features. Thus in scientific papers contour plots are much preferable to surface plots or twodimensional bar graphs. for example analytical error when the relevant source of `error’ for comparison of treatments is between fruit.29
3. there is a double equals sign (“==”). Explain clearly how error bars should be interpreted — SE limits. In plot(). use open plotting symbols. See help(plotmath) for detailed information on the plotting of mathematical expressions. with a single equals sign.8 Guidelines for Graphs
Design graphs to make their point tersely and clearly. which is exactly what one wants.
3.9 Exercises
1. 95% confidence interval.7 Plotting Mathematical Symbols
Both text() and mtext() will take an expression rather than a text string. Where points are dense. Draw graphs so that reduction and reproduction will not interfere with visual clarity. There is a further example in chapter 12. xlab="Radius". The final plot from
demo(graphics)
shows some of the possibilities for plotting mathematical symbols. and to important grouping in the data. Use graphs from which information can be read directly and easily in preference to those that rely on visual impression and perspective.
3. When there is extensive overlap of plotting symbols. In scatterplots the graph should attract the eye’s attention to the points that are plotted. overlapping points will give a high ink density. Use solid points. Use scatterplots in preference to bar or related graphs whenever the horizontal axis represents a quantitative effect.1:15 plot(x. Label as necessary to identify important features. ylab=expression(Area == pi*r^2))
Figure 13: The y-axis label is a mathematical expression. The data set huron that accompanies these notes has mean July average water surface elevations. SD limits. y.

2. 5.. Dept.rt(10.height)
This plots the mean level at year i against the mean level at year i-1. That is. c) As in the case of many time series. use
lag. b) Use the identify function to label points for the years that correspond to the lowest and highest mean levels. R. Cunningham. Look up the help for rnorm.html#DataVis has links to many different collections of information on statistical graphics. 1995.30
IGLD (1955) for Harbor Beach. Label the axes in untransformed units. on Lake Huron.yorku. type
identify(huron$year. Now generate a sample of size 10 from a normal distribution with mean 170 and standard deviation 4. and Donnelly. (b) and (c). that has mean heights for 1875-1772 only. In the bottom four rows. Station 5014. specify library(MASS)] 3.ca/SCS/StatResource. Check the distributions of head lengths (hdlngth) in the possum23 data set that accompanies these notes. Print out the numbers that you get. U. For example x <rchisq(10. c) a normal probability plot (qqnorm(possum$hdlngth)).1) will generate 10 random values from a t distribution with degrees of freedom one. (b) density plots. but taking samples from a uniform distribution rather than from a normal distribution.1) will generate 10 random values from a chi-squared distribution with degrees of freedom 1.rnorm(10). Michigan. Compare the following forms of display: a) a histogram (hist(possum$hdlngth)). Alongside on the same graphics page. Plot the graph of brain weight (brain) versus body weight (body) for the data set Animals from the
MASS package. D. National Oceanic and AtmosphericAdministration. F. log(brain weight) versus
log(body weight). now working with the logarithms of the data values. 8. B. and d) a density plot (plot(density(possum$hdlngth)). of Commerce. Repeat (a).
22
Source: Great Lakes Water Levels. all from a normal distribution. you might like to try it for some further distributions. Viggers. C. L... The statement x <. (Alternatively work with the vector LakeHuron from the datasets package. Use the row labels to label the points with the three largest body weight values. b) a stem and leaf plot (stem(qqnorm(possum$hdlngth)). display plots for samples of size 1000. show normal probability plots (section 3.) a) Plot mean.
3. Use mfrow() to set up the layout for a 3 by 4 array of plots. The function runif() generates a sample from a uniform distribution. In the middle 4 rows.height against year. Then repeat exercise 6 above. the mean levels are correlated from year to year. display plots for samples of size 100. In the top 4 rows. for 1860-198622. Try x <. Australian Journal of Zoology 43: 449-458. B. National Ocean Survey. [To access this data frame. 1860-1986.height. Trichosurus caninus Ogilby (Phalangeridae: Marsupialia). What shape do the points follow? *7. and print out the numbers you get. For the first two columns of the data frame hills. To see how each year's mean level is related to the previous year's mean level.plot(huron$mean.runif(10). Morphological variation among populations of the mountain brush tail possum. What are the advantages and disadvantages of these different forms of display? 4. K. Make normal probability plots for samples of various sizes from these distributions.huron$mean.4. (c) normal probability plots. Try x <. Label the axes appropriately. 6.math. If you find exercise 6 interesting. by default on the interval 0 to 1. examine the distribution using: (a) histogram.10 References
The web page http://www. press both mouse buttons simultaneously.labels=huron$year)
and use the left mouse button to click on the lowest point and highest point on the plot.
.2) for four separate `random’ samples of size 10.S. To quit. Comment on how the appearance of the plots changes as the sample size changes.
23
Data relate to the paper: Lindenmayer.

the time required for a simple discrimination task) and age are variables.
. but with a limitation to two conditioning factors or variables only. The two different symbols distinguish between low contrast and high contrast targets. The authors were mainly interested in visual recognition tasks that would be performed when looking through side windows. In this data frame. Lattice graphics
Lattice plots allow the use of the layout on the page to reflect meaningful aspects of data structure. hicon) are ordered factors. 2 = female) and agegp (1 = young. lo. T. R.
Figure 14: csoa versus it. Effects of car window tinting on visual performance: a comparison of elderly and young drivers. sex (1 = male. tint (level of tinting: no. To use lattice graphics. They offer abilities similar to those in the S-PLUS trellis library. it (inspection time. and Willson.. The older coplot() function that is in the base package has some of same abilities as xyplot( ). The different symbols are different contrasts. + = high contrast) are shown with different symbols.31
4.
4. The two targets (low. we might plot csoa against it for each combination of sex and agegp. In a simplified version of Figure 16 above. hi) and target (contrast: locon. Providing it is installed. both these packages must be installed. the time in milliseconds required to recognise an alphanumeric target).1 Examples that Present Panels of Scatterplots – Using xyplot()
The basic function for drawing panels of scatterplots is xyplot(). for each combination of females/males and elderly/young. N. csoa (critical stimulus onset asynchrony. For this simplified version. Nettlebeck. 1999. the grid package will be loaded automatically when lattice is loaded. 2 = an older participant. We will use the data frame tinting (supplied) to demonstrate the use of xyplot().e. J. The lattice package sits on top of the grid package. These data are from an experiment that investigated the effects of tinting of car windows on visual performance24. it would be enough to type:
xyplot(csoa ~ it | sex * agegp. data=tinting) # Simple use of xyplot()
Here is the statement used to get Figure 16. Ergonomics 42: 428-443. in the early 20s. i. in the early 70s) are factors. i. Figure 14 shows the style of graph that one can get from xyplot(). e.. White.
24
Data relate to the paper: Burns. M.

outer=TRUE. we do a plot (Figure 15) that uses different symbols (in black and white) for different levels of tinting. it is desirable to use auto. If on the same panel. separately for the two target types:
xyplot(csoa~it|sex*agegp. occur only for elderly males. as we might have expected. par.2 Some further examples of lattice plots
These are given with a minimum of explanation. +=low. Finally. The longest times are for the high level of tinting.settings=simpleTheme(pch=16.key=list(columns=2))
If colour is available. data=tinting. auto. A striking feature is that the very high values. >=high) are shown with different symbols.key to give a simple key.1 Plotting columns in parallel
Use the parameter outer to control whether the columns appear on the same or separate panels.key=list(columns=3))
Figure 15: csoa versus it. outer=FALSE. data=tinting."smooth")
The relationship between csoa and it seems much the same for both levels of contrast. The following puts smooth curves through the data. cex=2) )
. for each combination of females/males and elderly/young. data=grog) xyplot(Beer+Spirit+Wine ~ Year.32
xyplot(csoa~it|sex*agegp. groups=target. It is apparent that the long response times for some of the elderly males occur. outer=TRUE. groups=tint.2.
4. The following use the dataset grog from the DAAGxtras package:
library(DAAGxtras) xyplot(Beer+Spirit+Wine ~ Year | Country. groups=Country. data=tinting. panel=panel. groups=target. for both csoa and it.key=list(columns=3).
xyplot(csoa~it|sex*agegp.
4. different colours will be used for the different groups. data=grog. data=grog) xyplot(Beer+Spirit+Wine ~ Year | Country. with the low contrast target. type=c("p". auto. The different levels of tinting (no. auto.superpose.

34
(b) Construct a contour plot of chest versus belly and totlngth – use levelplot() or contourplot(). (c) Construct box and whisker plots for hdlngth, using site as a factor. (d) Use qqmath() to construct normal probability plots for hdlgth, for each separate level of sex and Pop. Is there evidence that the distribution of hdlgth varies with the level of these other factors. 13. The frame airquality that is in the datasets package has columns Ozone, Solar.R, Wind, Temp, Month and Day. Plot Ozone against Solar.R for each of three temperature ranges, and each of three wind ranges.

35

5. Linear (Multiple Regression) Models and Analysis of Variance 5.1 The Model Formula in Straight Line Regression
We begin with the straight line regression example that appeared earlier, in section 2.1.4. First, plot the data:
plot(distance ~ stretch, data=elasticband)

5.2 Regression Objects
An lm object is a list of named elements. Above, we created the object elastic.lm . Here are the names of its elements:
> names(elastic.lm) [1] "coefficients" "residuals" "effects" "rank"

By default the first, second and fourth plot use the row names to identify the three most extreme residuals. [If explicit row names are not given for the data frame, then the row numbers are used.]

5.3 Model Formulae, and the X Matrix
The model formula for the elastic band example was distance ~ stretch . The model formula is a recipe for setting up the calculations. All the calculations described in this chapter require the use of an model matrix or X matrix, and a vector y of values of the dependent variable. For some of the examples we discuss later, it helps to know what the X matrix looks like. Details for the elastic band example follow. The first 4 rows of the X matrix, with the y-vector alongside, is: X y Stretch (mm) Distance (cm) 1 46 148 1 54 182 1 48 173 1 50 166 … … ... The model matrix relates to the part of the model that appears to the right of the equals sign. The straight line model is y = a + b x + residual which we write as y=1 a+x b + residual

)
Model formulae are widely used to set up most of the model calculations in R.1 173-154.6 + 4.2 166-163. Least squares regression.1. data=elasticband) (Intercept) stretch 1 2 3 4 .6) and b ( = 4.lm as in section 5.8 = 18.1 = -0.55
ˆ y
(Fitted)  -63. .matrix(distance ~ stretch. The aim is to reproduce. glm. chooses a and b so that the sum of squares of the residuals is as small as possible.55 -63. as closely as possible. glm.9
y (Observed) Distance (mm) ! 148 182 173 166 …
ˆ y"y
(Residual) Observed – Fitted 148-145.55 -63. Thus:
> model. the first column by a and the second column by b.55 … Stretch 46 = 145.1 …
!
46 54 48 50 …
Note the use of the symbol
ˆ y [pronounced y-hat] for predicted values.55) that result from least squares regression X Stretch (mm) -63. etc.3 182-182. and adding.55 -63. Thus. Another name is predicted values. the x’s.6 + 4.
5. The X matrix then consists of a single column. y~x + fac + fac:x : lm. The function model..55 -63.3.1 Model Formulae in General
Model formulae take a form such as:
y~x+z : lm. the values in the y-column. which is the form of regression that we describe in this course.8 50 = 163.matrix(elastic.9 = 2. The residuals are the differences between the values in the y-column and the fitted values.6 + 4.
We might alternatively fit the simpler (no intercept) model.6 1 1 1 1 … 4. Fitted values are given by multiplying each column of the model matrix by its corresponding parameter.e. recall that the graph formula for a coplot that gives a plot of y against x for each different combination of levels of fac1 (across the page) and fac2 (up the page) is:
y ~ x | fac1+fac2
.1 48 = 154.matrix() prints out the model matrix.7 = 2.6 + 4. (If fac is a factor and x is a variable. .6 + 4.lm)
The following are the fitted values and residuals that we get with the estimates of a (= -63. For this we have y=x b+e
! where e is a random variable with mean 0.37
The parameters that are to be estimated are a and b. i. etc. is:
model. attr(.”assign”) [1] 0 1 1 1 1 1 46 54 48 50
Another possibility. aov.7 54 = 182. with elastic. fac:x allows a different slope for each different level of fac. Notice the similarity between model formulae and the formulae that are used for specifying coplots.

hard (hardness in `Shore’ units).L.formula(y~x+z)
or
formyxz <.g.nam[2])) > lm(formds.
26
The original source is O.3. This makes it straightforward to paste the argument together from components that are stored in text strings.
formyxz <. We first obtain a scatterplot matrix (Figure 18) :
Figure 18: Scatterplot matrix for the Rubber data frame from the MASS package.names(elasticband) > formds <. and tens (tensile strength in kg/sq m). as just demonstrated.4. data=elasticband) Call: lm(formula = formds. be a text string.38 *5.2 Manipulating Model Formulae
Model formulae can be assigned.
5."~".1 The data frame Rubber
The data set Rubber from the MASS package is from the accelerated testing of tyre rubber26. 119.
.3780 distance 0.1395
Note that graphics formulae can be manipulated in exactly the same way as model formulae. e. For example
> names(elasticband) [1] "stretch" "distance" > nam <. Table 6. Davies (1947) Statistical Methods in Research and Production.1 p.formula(“y~x+z”)
The argument to formula() can. Oliver and Boyd.4 Multiple Linear Regression Models
5. data = elasticband) Coefficients: (Intercept) 26. The variables are loss (the abrasion loss in gm/hr).formula(paste(nam[1].

obtained with termplot().lm(loss~hard+tens.lm().77e-011 71 on 2 and 27 degrees of freedom. data = Rubber) Residuals: Min 1Q Median 3. been fitted through the partial residuals.smooth) par(mfrow=c(1.194 14.84. partial=TRUE.lm) Call: lm(formula = loss ~ hard + tens. We proceed to regress loss on hard and tens.75 Max 65.161 -6.752 0.828 p-value: 1. data=Rubber) > options(digits=3) > summary(Rubber. note the use of termplot().1))
This plot raises interesting questions. F-statistic: Adjusted R-squared: 0.82 3Q 19.8e-14 1.2)) termplot(Rubber. the response is not linear with tensile strength.07 3.38 -14. at the mean of the contributions for the other term.0e-11 1. Figure 19) used the following code:
par(mfrow=c(1.98 -79.
.27 -7. Error t value Pr(>|t|) (Intercept) hard tens 885.3e-07
Residual standard error: 36.61 Coefficients: Estimate Std. smooth=panel.lm <.583 0.lm.571 -1.
Figure 19: Plot.5 on 27 degrees of freedom Multiple R-Squared: 0.374 61. at the upper end of the range. A smooth curve has.33 -11.39
The code is
library(MASS) pairs(Rubber) # if needed (the dataset Rubber is in the MASS package)
There is a negative correlation between loss and hardness.
Rubber. There is a clear suggestion that. in each panel.
In addition to the use of plot. showing the contribution of each of the two terms in the model.

The regressions on height and width give plausible results. as the analyses in Lalonde (1986) demonstrate.2 14.5243 0.lm3<-lm(weight~thick+height+width.7 8.
. height and width are so strong that if one tries to use more than one of them as a explanatory variables.5 19.2 Weights of Books
The books to which the data in the data set oddbooks (accompanying these notes) refer were chosen to cover a wide range of weight to height ratios.0 20.
> logbooks <.8 17.224 1. They contain very similar information. This is more than an academic issue.5 15.0 9.5 23.lm3)$coef Estimate Std.4.2 1075 940 625 400 550 600 450 450 300 690 400 250
Notice that as thickness increases.0124
> logbooks.5 29.121 1.552 0.678 -0.662 3.154 1.log(oddbooks) # We might expect weight to be > > summary(logbooks.data=logbooks) > summary(logbooks.907 0.114 3.data=logbooks) > summary(logbooks.708 0.lm1<-lm(weight~thick.6 10.07 0.7 19. the individual coefficients may be misleading.216 0.1 27.2 21.877 3.3 22.070 0.data=logbooks)
> logbooks. Even though regression equations from observational data may work quite well for predictive purposes.5 18.5 15. as is evident from the scatterplot matrix.434 1.5 12.719 0.356 0.313 2. Error t value Pr(>|t|) (Intercept) thick 9.lm1)$coef Estimate Std.lm2<-lm(weight~thick+height.219 13.8 13. Error t value Pr(>|t|) (Intercept) thick height -1.8 17.4 11.26e-04 # proportional to thick * height * width > logbooks.316 0. weight reduces. Here are the data:
> oddbooks thick height width weight 1 2 3 4 5 6 7 8 9 10 11 12 14 15 18 23 24 25 28 28 29 30 36 44 30. was statistically significant but with the wrong sign! The regression should be fitted only to that part of the data where values of the covariates overlap substantially. Error t value Pr(>|t|) (Intercept) thick height width -0.755 0.0 15.69 -1.6 23. and also from two comparable non-experimental “controls”. the coefficients are ill-determined. The design of the data collection really is important for the interpretation of coefficients from a regression equation.829 0.273 1.117
So is weight proportional to thick * height * width? The correlations between thick.6 12.465 0. They had data from experimental “treatment” and “control” groups. when comparison was with one of the non-experimental controls.7303 0.070 -0. The regression estimate of the treatment effect.117 0.lm2)$coef Estimate Std.472 0.40 5. while the coefficient of the regression on thick is entirely an artefact of the way that the books were selected.5 23.9 6.263 0.35e-08 -4.

00286 -0.04571 -0. and a column of values of rate2.1 Polynomial Terms in Linear Models
The data frame seedrates27 (DAAG package) gvies. Note that polynomial curves of high degree are in general unsatisfactory. C.41
Dehejia and Wahba demonstrate the use of scores (“propensities”) to identify subsets that are defensibly comparable.009911 0.
> seedrates.lm2 <. with fitted quadratic curve:
Figure 20: Number of grain per head versus seeding rate. for the barley seeding rate data. Curves are constructed from linear combinations of transformed values. data=seedrates) > summary(seedrates.066686 0.lm2) Call: lm(formula = grain ~ rate + I(rate^2). (1982) Effect of rates of seeding on barley grown for grain.5.5 Polynomial and Spline Regression
Calculations that have the same structure as multiple linear regression are able to model a curvilinear response. data = seedrates) Residuals: 1 2 3 4 5 0. a column of values of rate.000049 52. the number of barley grain per head.
plot(grain ~ rate.455694 0.09429 -0. constructed by joining low order polynomial curves (typically cubics) in such a way that the slope changes smoothly.00036 0. We will need an X-matrix with a column of ones. for each of a number of different seeding rates.50 0. New Zealand Journal of Agriculture 10: 133-136. Error t value Pr(>|t|) (Intercept) 24. data=seedrates) # Plot the data
Figure 20 shows the data.07294 0.12286 Coefficients: Estimate Std.lm(grain ~ rate+I(rate^2).
5. For this. are in general preferable. H. Summary details are in Maindonald. The propensity is then the only covariate in the equation that estimates the treatment effect.000171 0.01429
27
Data are from McLeod. C. both rate and I(rate^2) must be included in the model formula.
5.73 3. It is impossible to be sure that any method is giving the right answer. (1992).02138 0. with fitted quadratic curve.060000 rate I(rate^2) -0. J. Spline curves.80 -6.
.

5. e. R2 will in general increase as the range of the values of the dependent variable increases. 1. > # For this to work the values of seedrates$rate must be ordered. and compare them thus:
> seedrates. Type:
> model.) A predictive model is adequate when the standard errors of predicted values are acceptably small.972 for the straight line model.996. looked about right for the above data.160714 12.lm2.992 p-value: 0. an x3 column might be added."assign")
This was a (small) extension of linear models.187000 -1 -0. a quadratic curve. to handle a specific form of non-linear relationship.17^2.lm1<-lm(grain~rate.) based on 8 replicate results that were averaged to give each value for number of grains per head.f. > 1-pf(0. Any transformation can be used to form columns of the model matrix.17 on 35 d.predict(seedrates. but given the accuracy of these data it was not good enough! The statistic is an inadequate guide to whether a model is adequate.115 on 2 degrees of freedom Multiple R-Squared: 0. the F-value is now 5. The increase in the number of degrees of freedom more than compensates for the reduction in the F-statistic.17^2.data=seedrates) (Intercept) rate I(rate^2) 1 2 3 4 5 [1] 0 1 2 1 1 1 1 1 50 75 100 125 150 2500 5625 10000 15625 22500
attr(.4 on 1 and 35 degrees of freedom. Once the model matrix has been formed.16/0.07294
The F-value is large. of each of degrees 1 and 2. If we compare the change in the sum of squares (0. (R2 is larger when there is more variation to be explained.g. not when R2 achieves some magic threshold.
Again.0244
Finally note that R2 was 0. check the form of the model matrix.matrix(grain~rate+I(rate^2). i.2 What order of polynomial?
A polynomial of degree 2. How does one check? One way is to fit polynomials.42
Residual standard error: 0. F-statistic: > hat <.seedrates. hat)) > # Placing the spline fit through the fitted points allows a smooth curve.1607.Df Res.
> # However we have an independent estimate of the error mean > # square.Sum Sq Df Sum Sq F value Pr(>F) 1 2 0. we are limited to taking linear combinations of columns. on 1 df) with a mean square of 0. The estimate is 0.0039 256 on 2 and 2 degrees of freedom.lm2<-lm(grain~rate+I(rate^2). Thus. Even for any one context.228 0.026286 2 3 0.data=seedrates)
> anova(seedrates.172 (35 df). on 35 df. Adjusted R-squared: 0.e. and we have p=0.
5.024 .lm2) > lines(spline(seedrates$rate. but on this evidence there are too few degrees of freedom to make a totally convincing case for preferring a quadratic to a line. This may seem good.
.data=seedrates) > seedrates.lm1) Analysis of Variance Table Model 1: grain ~ rate + I(rate^2) Model 2: grain ~ rate Res. However the paper from which these data come gives an independent estimate of the error mean square (0. 35) [1] 0.

6 Using Factors in R Models
Factors are crucial for specifying R models that include categorical or “factor” variables.5. A.df <.ylab="Grains per head") new. ylim = c(15. D.5. xlim = c(50. While one would not usually handle these calculations using an lm model."lwr"]. Confidence bounds for a fitted line spread out more slowly. F.df) lines(rate.lm5."fit"]) > lines(cars$speed.data=cars)
> ci5<-predict(cars.scale"
> lines(cars$speed.lm) > lines(cars$speed. (1994). A Handbook of Small Data Sets. pch = 16. D.df)
The extrapolation has deliberately been taken beyond the range of the data. The places where the cubics join are known as `knots’.5).5. It turns out that. Lunn.se.lty=2) > lines(cars$speed. but are even less believable!
5.data=cars) > hat<-predict(cars. 175).. The fitting of polynomial functions was a simple example of this.lm2. it makes a simple example to illustrate the choice of a baseline level.data.predict(seedrates. Spline functions variables extend this idea further. once the knots are fixed. and depending on the class of spline curves that are used."upr"].3 Pointwise confidence bounds for the fitted curve
Here is code that gives pointwise 95% confidence bounds.hat2$lower. Ostrowski.df.frame(rate = c((4:14) * 12.interval="confidence". data = seedrates. The splines that I demonstrate are constructed by joining together cubic curves. interval="confidence") hat2 <. newdata = new.
28
Data are from Hand.hat.lower=pred2[.ci5$fit[. 22). give output in which some of the numbers are different and must be interpreted differently. hat2$upper.lty=2)
5. data = seedrates) pred2 <.lty=3) # NB assumes values of speed are sorted # try for a closer fit (1 knot) # By default. where each basis function is a transformation of the variable. in such a way the joins are smooth. Consider data from an experiment that compared houses with and without cavity insulation28.frame(fit=pred2[. Daly. spline functions of a variable can be constructed as a linear combination of basis functions.ci5$fit[.lm5 <. J."lwr"].data. Different choices... The data frame cars is in the datasets package.ci5$fit[.xlab="Seeding rate".fit=TRUE) > names(ci5) [1] "fit" "se..lm(grain ~ rate + I(rate^2). upper=pred2[. although they fit equivalent models. readers of this document will be used to the idea that it is possible to use linear models to fit terms that may be highly nonlinear functions of one or more of the variables.4 Spline Terms in Linear Models
By now.fit" "df" "residual.lm2 <. eds.lty=2) lines(rate. Chapman and Hall."fit"]. there are no knots
> cars. in order to show how the confidence bounds spread out.data=cars) > library(splines) > cars.. E.lm(dist~bs(speed.5)) seedrates. hat2$fit) lines(rate."upr"]) attach(new.43 5. Note that these do not combine to give a confidence region for the total curve! The construction of such a region is a much more complicated task!
plot(grain ~ rate.
.lm<-lm(dist~bs(speed). and a set of contrasts.lty=2) detach(new.
> plot(dist~speed.

c(10225. Here are the X and the y that are used for the calculations. 7980+0 7980+0 .. 12669. 14683. 8541. Hence: Average for Insulated Houses = 7980 To get the estimate for uninsulated houses take 7980 + 3103 = 10993.factor(c(rep("without". corr=F) Call: lm(formula = kWh ~ insulation) Residuals: Min -4409 1Q Median -979 132 3Q 1575 Max 3690
Coefficients: Estimate Std.022. 8).... then 7 with > kWh <..8e-07 0.05) that we can distinguish between the two types of houses. which may be taken to indicate (p < 0.
> insulation <.022
Residual standard error: 2310 on 13 degrees of freedom Multiple R-Squared: 0. 8022)
To formulate this as a regression model.. 9708-7980 6700-7980 . 10315. 8022) > insulation. we take kWh as the dependent variable.
The p-value is 0. 9708.. Note that the first eight data values were all withouts: Contrast 7980 1 1 . 10315. 12086. 8162.44
We begin by entering the data from the command line:
insulation <.. then 7 with # `with’ precedes `without’ in alphanumeric order.factor(c(rep("without". & is the baseline kWh <..59 5. 6584...0223 F-statistic: 6. 6584.. 8017. 9708 6700 Residual 10225-10993 10689-10993 ..1 The Model Matrix
It often helps to keep in mind the model matrix or X matrix. we need to know that the factor levels are.. 4307. 8541..341. But what does the “intercept” of 7890 mean.lm(kWh ~ insulation) > summary(insulation. taken in alphabetical order.6.03 2. and the factor insulation as the explanatory variable.. The standard error of the difference is 1196.. 7))) # 8 without.. So with comes before without.. 12467. 7))) > # 8 without.lm <. and that the initial level is taken as the baseline.c(10225. 8). rep("with".
5. Error t value Pr(>|t|) (Intercept) insulation 7890 3103 874 1196 9. Add to get 7980+3103=10993 7980+3103=10993 kWh Compare with 10225 10689 . Adjusted R-squared: 0. 4307. 6700.73 on 1 and 13 degrees of freedom. 14683. 10689. 8162. and with is the baseline. 0 0 . 12086. + 12669... 1 1 . 9708. 12467. 8017. and what does the value for “insulation” of 3103 mean? To interpret this.. by default. 3103 1 1 . 10689..
. rep("with".lm.. 6700.29 p-value: 0.

"contr. With the sum contrasts the baseline is the overall mean. Here is the output from use of the `sum’ contrasts29:
> options(contrasts = c("contr. Adjusted R-squared: 0.4e-10 0. so that it becomes the baseline.lm(kWh ~ insulation) > summary(insulation.45
Type
model. “mean for with”) = 9442 To get the estimate for uninsulated houses (the first level). levels=c("without".0223 F-statistic: 6. there are other choices of contrasts. In technical jargon.C(insulation. i. With the “sum” contrasts the baseline is the mean over all factor levels. One obvious alternative is to make without the first factor level.lm.]
30
In general.73 on 1 and 13 degrees of freedom.poly".6.
*5. the user has to calculate it as minus the sum of the remaining effects. with a statement such as: insulation <. [Use the function ordered() to create ordered factors. The effect for the first level is omitted.sum". digits = 2) # Try the `sum’ contrasts > insulation <.341. For this.factor(insulation.poly").2 Other Choices of Contrasts
There are other ways to set up the X matrix. Error t value Pr(>|t|) (Intercept) insulation 9442 1551 598 598 15.
29
The second string element. "with")) # Make `without' the baseline > insulation.022
Residual standard error: 2310 on 13 degrees of freedom Multiple R-Squared: 0. The sum contrasts are sometimes called “analysis of variance” contrasts.78 2.
.29 p-value: 0. corr=F) Call: lm(formula = kWh ~ insulation) Residuals: Min -4409 1Q Median -979 132 3Q 1575 Max 3690
Coefficients: Estimate Std.
It is possible to set the choice of contrasts for each factor separately. So the effect for the second level (`with’) is -1551.e. Novices should avoid them30.matrix(kWh~insulation)
and check that it gives the above model matrix. baseline="without") # Make `without’ the baseline
Another possibility is to use what are called the “sum” contrasts. These are not at all intuitive and rarely helpful. Thus to get the estimate for insulated houses (the first level). take 9442 + 1551 = 10993 The `effects’ sum to one. even though S-PLUS uses them as the default.59 7.relevel(insulation.1551 = 7980. use either the treatment contrasts or the sum contrasts. "contr. is the default setting for factors with ordered levels. contr=treatment)
Also available are the helmert contrasts. take 9442 .
Here is the interpretation: average of (mean for “without”. specify:
insulation <.lm <.

lm2 <. . or interactions between factors. For the second group (Stellena styx.21858 -0.lm2) (Intercept) factor(species) logweight 1 2 . x1 = 0) the constant term is a1 and the slope is b1.111 on 14 degrees of freedom Multiple R-Squared: 0. data=dolphins)
Check the model matrix:
> model. while for the second group (Delphinus delphis. x1 = 1) the constant term is a1 + a2. .133 0. In the example that follows. 8 .treatment" 1 0 3.7 Multiple Lines – Different Regression Lines for Different Species
The terms that appear on the right of the model formula may be variables or factors.e. .989 1 1 1 1 3.5e-07 # Plot the data
Residual standard error: 0.lm1 <.522 0.46
5. 16 [1] 0 1 2 attr(. i. A and B:
> plot(logheart ~ logweight. Then possibilities we may want to consider are: A: A single line: y = a + b x2 B: Two parallel lines: y = a1 + a2 x1 + b x2 [For the first group (Stellena styx. x1 = 1) the constant term is a1 + a2. . We take x2 to be body weight. we had weights for a porpoise species (Stellena styx) and for a dolphin species (Delphinus delphis).838.lm(logheart ~ logweight. Error t value Pr(>|t|) (Intercept) logweight 1.827 p-value: 6.555 3.325 1. We take x1 to be a variable that has the value 0 for Delphinus delphis. data = dolphins) > summary(cet.
For model B (parallel lines) we have
> cet.133 2. Here we take advantage of this to fit different lines to different subsets of the data. data = dolphins) Residuals: Min 1Q Median 0. and 1 for Stellena styx. or interactions between variables and factors.024 6. Adjusted R-squared: 0. .54 8. .04981 Max 0.51e-007 F-statistic: 72.] C: Two separate lines: y = a1 + a2 x1 + b1 x2 + b2 x1 x2 [For the first group (Delphinus delphis.lm(logheart ~ factor(species) + logweight."contrasts") [1] "contr.00274 3Q 0. data=dolphins) > options(digits=4) > cet.08249 Coefficients: Estimate Std.738
.951 attr(.15874 -0. corr=F) Call: lm(formula = logheart ~ logweight.6 on 1 and 14 degrees of freedom.lm1. x1 = 0) the constant term is a1.] We show results from fitting the first two of these models.matrix(cet."assign") 1 0 3.52 0. and the slope is b1 + b2.

For the remaining two blocks shelter effects. P.001194
31
Data relate to the paper: Snelgar. Adjusted R-squared: 0.3710 0. Further details. (If this is not specified. But the F-statistics and p-values will be wrong.1176 0.527 -1. P.0411 1. including a diagram showing the layout of plots and vines and details of shelter.0877
Residual standard error: 0.2 Shading of Kiwifruit Vines
These data (yields in kilograms) are in the data frame kiwishade that accompanies these notes.2788 0.)
> kiwishade$shade <.2627 Max 1.331 1.06 Pr(>F) 1 2496. Influence of time of shading on flowering and yield of kiwifruit vines.17 20.aov) Df Cult Date Residuals 2 Sum Sq Mean Sq F value 909.179e-09 9..846 on 2 and 27 degrees of freedom.8. shading from December to February.relevel(kiwishade$shade. Each treatment appeared once in each of the three blocks.J. in one case from the east and in the other case from the west. Journal of Horticultural Science 67: 481-487. and shading from February to May.51
464. They are from an experiment31 where there were four treatments . Manson. the four vines within a plot constitute subplots.data=kiwishade) > summary(kiwishade.2788 25.65 47. Results are given for each of the four vines in each plot.35 125. The two papers have different shorthands (e.
> VitC. In experimental design parlance.4180 -0. shading from August to December. Variation in the vitamin C levels seems relatively consistent between cultivars.57 86..data=cabbages) > summary(VitC.3690 -1. are in Maindonald (1992).40
*5.48
Residuals: Min 1Q Median 3Q 0.6234 on 27 degrees of freedom Multiple R-Squared: 0. > help(cabbages) > names(cabbages) [1] "Cult" "Date" "HeadWt" "VitC" > coplot(HeadWt~VitC|Cult+Date.1944 0.g. The northernmost plots were grouped in one block because they were similarly affected by shading from the sun.15 2496.0320 -0.074879
3 1394. 1992.0002486 56 2635. W.1971 0.J.772 <2e-16 0.aov) Error: block:shade Df block shade Residuals 2 6 Sum Sq Mean Sq F value 172.aov<-aov(yield~block+shade+Error(block:shade). ref="none") > ## Make sure that the level “none” (no shade) is used as reference > kiwishade.01591 F-statistic: 4. Martin.30 454.P.15 53.no shading.84 22.aov<-aov(VitC~Cult+Date. were thought more important. Error t value Pr(>|t|) (Intercept) grouptrt1 grouptrt2 5.
.data=cabbages)
# cabbages is from the MASS package
Examination of the plot suggests that cultivars differ greatly in the variability in head weight.93 Pr(>F) 4.2641.0710 -0. The block:shade mean square (sum of squares divided by degrees of freedom) provides the error term. Sept-Nov versus Aug-Dec) for describing the time periods for which the shading was applied.6609 0.2112 0. one still gets a correct analysis of variance breakdown.2096 p-value: 0.0060 Coefficients: Estimate Std.4940 0.

189. Does the quadratic curve give a useful improvement to the fit? If you have studied the dynamics of particles. and sum contrasts. carry out a regression of strength on SpecificGravity and Moisture. regress time on distance and climb. 54. What is the key difference between the two sets of data? Using the data frame beams (in the data sets accompanying these notes). Mater.18
5. 114. the values are: stretch (mm): 46. Using an outcome variable
y <. Hunter. 249. 228. the values are: stretch (mm): 25. and plot the line. For the data set elastic2. 50. Type
hosp<-rep(c(”RNC”.430000 3. formally. Fit a line to this relationship. obtained by supplying the name of the lm object as the first parameter to plot(). Using the data frame hills (in package MASS). Then try fitting and plotting a quadratic curve. 217.10. Repeat the exercise for the helmert contrasts. Do the two sets of results appear consistent. 217.49
Error: Within Df Sum Sq Mean Sq F value Pr(>F) Residuals 36 438.”Mater”). whether one needs different regression lines for the two data frames elastic1 and elastic2. Compare the two sets of results.428333 12. 6.5. 50 30.281667 -7. Which of these regression equations would you prefer? Use the method of section 5. 50. 55.58 > coef(kiwishade.7 check the form of the model matrix (i) for fitting two parallel lines and (ii) for fitting two arbitrary lines when one uses the sum contrasts. 208.9)
. 127. can you find a theory that would tell you how stopping distance might change with speed?5. 42.030833 -10. What can you learn from the diagnostic plots that you get when you plot the lm object? Try also regressing log(time) on log(distance) and log(climb). 45. treatment contrasts. including the same rubber band. Here are two sets of data that were obtained the same apparatus. as the data frame elasticband. 35. 249. For each of the data sets elastic1 and elastic2.2) hosp fhosp<-factor(hosp) levels(fhosp)
Now repeat the steps involved in forming the factor fhosp.aov) (Intercept) : (Intercept) 96. determine the regression of stretch on distance. 178. 40.7 to determine. stopping distance) versus speed. What does this indicate? Using the data frame cars (in the datasets package). Do this using helmert contrasts. Using a different symbol and/or a different colour.”Hunter”.8.5327 block:shade : blocknorth 0. 48.9 Exercises
1. 291.e. In section 5.3. 187. this time keeping the factor levels in the order RNC. 44. 52 distance (cm): 183. plot the data from the two data frames elastic1 and elastic2 on the same graph. Use contrasts(fhosp) to form and print out the matrix of contrasts.993125 Within : numeric(0) blockwest shadeAug2Dec shadeDec2Feb shadeFeb2May -3. 150. Carefully examine the regression diagnostic plot. plot distance (i. 60 distance (cm): 71. In each case determine (i) fitted values and standard errors of fitted values and (ii) the R2 statistic.c(2. For the data set elastic1. 196.

Statistical design. examine the dependence of y (amount of heat produced) on x1. B. 1999. Weisberg. A. 397–402. New Zealand Journal of Agricultural Research 35: 121-141.. Lee. A. To what extent is it potentially misleading? Also do the analysis where the block:shade term is omitted altogether. Comment: Aspects of diagnostic regression analysis. Introductory Statistics with R. B.. Comment on what you get. D. Fox. but replacing Error(block:shade) with block:shade. perhaps by taking log(x/(100-x))? What alternative strategies one might use to find an effective prediction equation? In the data set pressure (datasets package). G. Modern Applied Statistics with S. E. and Mark. New York. Geology 11: 195-1976. F. For which choice(s) of contrasts do the parameter estimates change when you re-order the factor levels? In the data set cement (MASS package). L. quadratic. W. analysis and presentation issues. the direction of the curvature changes. Maindonald J H and Braun W J 2007. do they require transformation. Consider transformation of pressure. repeating the fit for each of the three different choices of contrasts. 1986. N. 1996. A logarithmic transformation is too extreme.. Causal effects in non-experimental studies: re-evaluating the evaluation of training programs. Lalonde. S.
5. Wiley. American Economic Review 76: 604-620. Multivariate and Tree-Based Methods
. examine the dependence of pressure on temperature. C. R. Venables.. 1999. Wiley. 1986. and (b) a cubic curve. Statistics in Medicine 15: 361-387. R. As the explanatory variables are proportions. Evaluating Assumptions and Adequacy. 1988. and Weisberg. Cook.5. Transformations Unmasked. R. x2. P. and Measuring and Reducing Errors.
Maindonald J H 1992.2). An R and S-PLUS Companion to Applied Regression. Data Analysis and Graphics Using R – An Example-Based Approach. Comment on that analysis. 4th edn 2002. 2nd edn. Cambridge University Press. Statistical Science 1. [Transformation of temperature makes sense only if one first converts to degrees Kelvin. Dehejia. Sage Books. Applied Linear Regression. Springer. Improper use of regression equations in the earth sciences. Applied Regression including Computing and Graphics. Tutorial in Biostatistics. P 2002. Which of the three possibilities (line. and Wahba. x3 and x4 (which are proportions of four constituents). with accompanying 95% confidence bounds. Williams. S. J 2002. Journal of the American Statistical Association 94: 1053-1062. 1985. Begin by examining the scatterplot matrix. Comment on the output that you get from summary(). Harrell. curve) is most plausible? Can any of them be trusted? *Repeat the analysis of the kiwishade data (section 5. D. Technometrics 30: 311-318.10 References
Atkinson. Dalgaard. 1983. S.50
fit the model lm(y~fhosp). Evaluating the economic evaluations of training programs.
Atkinson. C. Multivariable Prognostic Models: Issues in Developing Models. with accompanying 95% pointwise confidence bounds. 2nd edn. What family of transformations might one try? Modify the code in section 5. D.H.8. and Ripley. New York. Springer.3 to fit: (a) a line. K.

2] ~ possum. A plot in which the sites and/or the sexes are identified will draw attention to any very strong structure in the data. col=palette()[as. With biological measurement data.6. The first principal component is the component (linear combination of the initial variables) that explains the greatest part of the variation.. Cunningham. groups=possum$site. 1995. The measure of variation used is the sum of the variances of variables.
. relative variability).princomp(log(possum[here. and hence with the variance-covariance matrix. R. it is usually desirable to begin by taking logarithms.key=list(columns=3)) # Principal components # Print scores on second pc versus scores on first pc. The idea is to replace the initial set of variables by a small number of “principal components” that together may explain most of the variation in the data. gives a greater weight to variables that have a large variance. scaling to give equal variances is unnecessary. D.6:14]. This may draw attention to gross errors in the data. L. C. This is because the ratio of largest to smallest value is relatively small.51
6. Here are some of the scatterplot matrix possibilities:
pairs(possum[. never more than 1. and so on.na(possum$footlgth) print(sum(!here)) # We need to exclude missing values # Check how many values are missing
We now look (Figure 21) at particular views of the data that we get from a principal components analysis:
possum.e. The second principal component is the component that. B.6:14]. The common alternative –scaling variables so that they each have variance equal to one – is equivalent to working with the correlation matrix. auto. and Donnelly. Viggers. Morphological variation among columns of the mountain brushtail possum.6:14])) # by populations and sex. for some or all of the variables. and because the standard deviations are typically fairly comparable. perhaps after scaling so that they each have variance one.prc$scores[. An analysis that works with the unscaled variables. F. It is good practice to begin by examining relevant scatterplot matrices.prc$scores[..1 Multivariate EDA.. trapped at seven sites from southern Victoria to central Queensland32.prc <. for all variables. The data set possum has nine morphometric measurements on each of 102 mountain brushtail possums.
32
Data relate to the paper: Lindenmayer. The standard deviations then measure the logarithm of relative change. Australian Journal of Zoology 43: 449-458. Taking logarithms of these data does not make much difference to the appearance of the plots. K. Because all variables measure much the same quantity (i. among linear combinations of the variables that are uncorrelated with the first principal component.1]|possum$Pop[here]+possum$sex[here]. and Principal Components Analysis
Principal components analysis is often a useful exploratory tool for multivariate data. Multivariate and Tree-based Methods 6. Trichosurus caninus Ogilby (Phalangeridae: Marsupiala). col=palette()[as. B. identified by site xyplot(possum. explains the greatest part of the remaining variation.integer(possum$site)]) here<-!is.integer(possum$sex)]) pairs(possum[. For example one site may be quite different from the others.

for the possum data. In hierarchical agglomeration each observation starts as a separate group. to distinguish animals from different sites. It is “unsupervised” because the clusters are not known in advance.lda(site ~ hdlngth+skullw+totlngth+ # Only if not already attached. by population and by sex. Groups that are “close” to one another are then successively merged. How does one get the initial classification? Typically.e. we may wish to predict.
6. i. and algorithms based on iterative relocation.3 Discriminant Analysis
We start with data that are classified into several groups. Our interest is in whether it is possible. There are two types of algorithms – algorithms based on hierachical agglomeration.na(possum$footlgth) > possum. which future patients will remain free of disease symptoms for twelve months or more. In the language of Ripley (1996).
> library(MASS) > here<. A cruder distinction is between populations. Because it has little on the distribution of variable values. based on prognostic measurements and outcome information for previous patients. using the lda() function from the Venables & Ripley MASS package. The output yields a hierarchical clustering tree that shows the relationships between observations and between the clusters into which they are successively merged. In iterative relocation.
6.lda <. sites in Victoria (an Australian state) as opposed to sites in other states (New South Wales or Queensland).2 Cluster Analysis
In the language of Ripley (1996)33. with a choice of methods available. cluster analysis is a form of unsupervised classification. on the basis of morphometric measurements. and want a rule that will allow us to predict the group to which a new data value will belong.52
Figure 21: Second principal component versus first principal component. The mva package has the cluster analysis routines. I have not thought it necessary to take logarithms. the algorithm starts with an initial classification. Here are calculations for the possum data frame. by a prior use of a hierarchical agglomeration algorithm. that it then tries to improve. I discuss this further below. For example. The function kmeans() (k-means clustering) implements iterative relocation.
. The function dist() calculates distances. A judgement is then needed on the point at which further merging is unwarranted. The function hclust() does hierarchical agglomerative clustering.!is. our interest is in supervised classification.
33
References are at the end of the chapter.

1860 1. as in Figure22 3. Tree-based methods seem more suited to binary regression and classification than to regression with an ordinal or continuous dependent variable.lda. A tree-based regression approach is available for use for regression problems.
6.7578 > > plot(possum. from MASS package glass. Output includes a “decision tree” that is immediately useful for prediction. method.9372 # Examine the singular values 3. One advantage of such methods is that they automatically handle non-linearity and interactions.lda$svd [1] 15.tree)
. Clearly canonical variates after the third will have little if any discriminatory power.
library(rpart) # Use fgl: Forensic glass fragment data. for the canonical variates in turn.5078 1. Note that there may be interpretative advantages in taking logarithms of biological measurement data. Tree-based models. or discrimination. dimen=3) > # Scatterplot matrix for scores on 1st 3 canonical variates. text(glass. The standard against which patterns of measurement are commonly compared is that of allometric growth. subset=here) > options(digits=4) > possum. may be suitable for regression and classification problems when there are extensive data. The reader may wish to repeat the above analysis. Differences between different sites are then indicative of different patterns of allometric growth.7772
Figure 22: Scatterplot matrix of first three canonical variates.lda() to get (among other information) scores on the first few canonical variates. but working with the logarithms of measurements.1420 0. which implies a linear relationship between the logarithms of the measurements.data=possum. data=fgl) plot(glass. also known as “Classification and Regression Trees” (CART). One can use predict.tree). The singular values are the ratio of between to within group sums of squares.tree <.53
+ taillgth+footlgth+earconch+eye+chest+belly.4 Decision Tree models (Tree-based models)
We include tree-based classification here because it is a multivariate supervised classification.rpart(type ~ RI+Na+Mg+Al+Si+K+Ca+Ba+Fe.

uk/pub/SWin/rpartdoc. An Introduction to Recursive Partitioning Using the RPART Routines. J.. Cambridge UK. and (ii) the six different types (Type) of car. Cambridge University Press. 1999. B.5 Exercises
1. Investigate discrimination between plagiotropic and orthotropic species in the data set leafshape34. Therneau.ox. and Maindonald. containing de-identified data on the survival status of patients diagnosed with AIDS before July 1 1991. D. Tree architecture in relation to leaf dimensions and tree stature in temperate and tropical rain forests. Journal of Ecology 87: 1012-1024.
6. 3.3. 1992. It integrates cross-validation with the algorithm for forming trees.zip Venables.tree)
To use these models effectively. 1997.ac. 2. Friedman. T.H.e. Now create a new data set in which binary factors become columns of 0/1 data. and Ripley. This is one of two documents included in: http://www. Plot a scatterplot matrix of the first three principal components. R. Cambridge University Press. 5. T. Using the columns of continuous or ordinal data. D. and Atkinson. Drawing. and about crossvalidation.A. D. determine scores on the first and second principal components.
6. apply principal components analysis to the scores for Composition. New York. 4. Methods for reduction of tree complexity that are based on significance tests at each individual node (i. Springer. Use tree-based classification (rpart()) to identify major influences on survival. Statistical Models in S. Data Analysis and Graphics Using R – An Example-Based Approach.6 References
Chambers. The Atkinson and Therneau rpart (recursive partitioning) package is closer to CART than is the S-PLUS tree library. Using the data set painters (MASS package). Investigate the comparison between (i) USA and non-USA cars. J. N. W. Pacific Grove CA. 2nd edn 1997. Pattern Recognition and Neural Networks. Examine the loadings on the first three principal components. 2nd edn. Repeat the calculations of exercises 1 and 2. The MASS package has the Aids2 data set. Ripley. B.
34
Data relate to the paper: King. Wadsworth and Brooks Cole Advanced Books and Software. Additive logistic regression: A statistical view of boosting. as in section 6. and include these in the principal components analysis. branching point) typically choose trees that over-predict. J. Available from the internet. but this time using the function lda() from the MASS package to derive canonical discriminant scores.54
summary(glass. M.stats. and Tibshirani. Colour. and Hastie.
. J. and Expression. Modern Applied Statistics with S-Plus. The data set Cars93 is in the MASS package. 1996. T. Maindonald J H and Braun W J 2008. J. M. you also need to know about approaches to pruning trees. using different colours or symbols to identify the different schools.. (1998). Hastie. E.

1 Subsets of Vectors
Recall (section 2. <=. meaning “keep on repeating the sequence until the length is length. 6.5). 15. Entering 15:5 will generate the sequence in the reverse order. !=. c(4. …. John=185.3.4. >. Users who do not carefully consider implications for expressions that include Nas may be puzzled by the results.out.c(4.4)) we could enter rep(c(2. for vectors that have named elements:
> c(Andreas=178.5).NA) > is. note that the function rep() has an argument length.4. The first four compare magnitudes. 4) thus:
> rep(c(2. Be sure to use is. Thus suppose we want to extract values of x that are greater than 10. note that x==NA generates NA.3)).out. == tests for equality.3.6.5).3.3.3.1. and != tests for inequality. Jeff=183)[c("John".3. As x==NA gives a vector of NAs.”
7.6.
. the missing value symbol is NA.
> rep(c(2.4.2 Patterned Data
Use 5:15 to generate the numbers 5. The elements that are extracted are those for which the logical value is T. rep(4. Specify a vector of logical values.5).c(1. enter rep(c(2.2.5).4. The following demonstrates a third possibility. and otherwise FALSE TRUE # All elements are set to NA [1] FALSE FALSE FALSE [1] NA NA NA NA
Missing values in subscripts: In R
y[x>2] <. 3.3)."Jeff")] John Jeff 185 183
A vector of names has been used to extract the elements.4)). numeric or character. ==. To repeat the sequence (2.4) [1] 2 3 5 2 3 5 2 3 5 2 3 5 >
If instead one wants four 2s.55
*7.2) two common ways to extract subsets of vectors: Specify the numbers of the elements that are to be extracted. >=.1.na(x) > x==NA > NA==NA [1] NA # TRUE for when NA appears. then four 5s.
7. In addition to the above.1 Vectors
Recall that vectors may have mode logical. For example
> x <. Specifically. c(4. you get no information at all about x.2 Missing Values
In R.3.5). in place of c(4. each=4)
Note further that. This applies also to the relations <.
7. Any arithmetic operation or relation that involves NA generates an NA.4)) [1] 2 2 2 2 3 3 3 3 5 5 5 5 # An alternative is rep(c(2. then four 3s.x[x>2]
generates the error message “NAs are not allowed in subscripted assignments".na(x) to test which values of x are NA. enter rep(c(2. So a further possibility is that in place of rep(c(2. 5) four times over. One can use negative numbers to omit elements. R Data Structures 7.5).4) we could write rep(4.

use of e.3. To get information on the data sets that are included in the datasets package.36667 No. of the properties of matrices.3 Data frames
The concept of a data frame is fundamental to the use of most of the R modelling and graphics functions."yield") we could have written.
7.2 Data Sets that Accompany R Packages
Type in data() to get a list of data sets (mostly data frames) associated with all packages that are in the current search path. To attach any other installed package.4).33333
The first column holds the row labels. In most packages.1 Extraction of Component Parts of Data frames
Consider the data frame barley that accompanies the lattice package:
> names(barley) [1] "yield" "variety" "year" "Duluth" "Waseca" "site" "University Farm" "Morris" > levels(barley$site) [1] "Grand Rapids" [5] "Crookston"
We will extract the data for 1932.36667
120 Wisconsin No.6. Other packages must be explicitly installed.56
Users are advised to use !is. specify
data(package="datasets")
and similarly for any other package. at the Duluth site. 38 29. more simply. All elements of any column must however have the same mode. 457 22. 462 22. i. or all character.60000 No. Data frames where all columns hold numeric data have some.86667 Svansota 22.matrix(). and several other packages.50000 Peatland 31. 475 27. data(airquality) to load the data set airquality (datasets package) is unnecessary The out-of-the-box Windows and other binary distributions include many commonly required packages.na(x) to limit the selection. to those elements of x that are not NAs.3.46667 Trebi 30. c(2.e.56667 Glabron 25. which comes with the default distribution. use the library() command.
7. We will have more to say on missing values in the section on data frames that now follows. are automatically attached at the beginning of the session. + c("variety".
> Duluth1932 <. A data frame is a generalisation of a matrix. Lists are discussed below. in which different columns may have different modes. on one or both sides as necessary. will be used from time to time. in section 7.barley[barley$year=="1932" & barley$site=="Duluth"."yield")] variety 66 72 78 84 90 96 102 108 114 yield Manchuria 22. the MASS package.
.g.. use as.23333 Velvet 22. In place of c("variety". but not all. For remaining sections of these notes. The base package. data from an attached package are automatically available. all numeric or all factor. There are important differences that arise because data frames are implemented as lists. which in this case are the numbers of the rows that have been extracted.
7.70000 No. To turn a data frame of numeric data into a matrix of numeric data.

If you use any missing value symbols other than the default (NA). The default separator is white space.csv".table() is straightforward for reading in rectangular arrays of data that are entirely numeric.table("a:/myfile. The row names are the city names. it is sometimes necessary to use tab (“\t”) or comma as the separator. Problems may arise when small mistakes in the data cause R to interpret a column of supposedly numeric data as character strings. specify sep="\t". They have an attribute levels that holds the level names.table() expects missing values to be coded as NA.strings=c(“NA”. the column is by default stored as a factor with as many different levels as there are unique text strings35.2 Missing values when using
read. as in the above example. use read.”*”. the first column (country) has the name of the country. The data frame islandcities that accompanies these notes holds the populations of the 19 island nation cities with a 1995 urban centre population of 1. you need to make this explicit see section 7.4."). unless you set na.
36
Specifying as.3.
. 36 If the column is later required for use as a factor in a model or graphics formula. They are stored as integer vectors. When.table()
With data from spreadsheets37. They provide an economical way to store vectors of character strings in which there are many multiple occurrences of the same strings.".57
7. e. Factors have a dual identity.4.is = TRUE. and the second column (population) has the urban centre population.csv and is on drive a:.is = T prevents columns of (intended or unintended) character strings from being converted into factors. Here is a table that gives the number of times each country occurs
Australia Cuba Indonesia Japan Philippines Taiwan United Kingdom 3 1 4 6 2 1 2
35
Storage of columns of character strings as factors is efficient when a small number of distinct strings that are of modest length are each repeated a large number of times.) and blank (in a case where the separator is something other than a space) will cause to whole column to be treated as character data.
7. which are automatically turned into factors.1 Idiosyncrasies
The function read.". you will probably want to set na. it may be necessary to make it into a factor at that time. sep=". More crucially.table()
The function read. with each of the values interpreted according to the information that is in the table of levels38. Users who find this default behaviour of read.g.csv text file with comma as the separator.strings=c(". in millions. Otherwise any appearance of such symbols as *. For example there may be an O (oh) somewhere where there should be a 0 (zero).3 Separators when using read.4 Data Entry Issues
7. There may be multiple missing value indicators.table() confusing may wish to use the parameter setting as. one of the columns contains text strings.6. If you have a text file that has been output from SAS.") to read the data into R. To set tab as the separator.
37
One way to get mixed text and numeric data across from Excel is to save the worksheet in a .4. [But watch that none of the cell entries include commas.""). period(.
7.strings to recognise other characters as missing value indicators. Some functions do this conversion automatically. they have a central role in the incorporation of qualitative effects into model and graphics formulae.4.4 million or more.]
38
Factors are column objects that have mode numeric and class “factor”. This copes with any spaces which may appear in text strings. If for example file name is myfile. The "" will ensure that empty cells are entered as NAs. or an el (l) where there should be a one (1).5 Factors and Ordered Factors
We discussed factors in section 2.
7. na.2 below.

stress<"medium" TRUE FALSE FALSE TRUE > ordf.6 Ordered Factors
Actually. and of any classes from which it inherits. Try
> stress. R sorts the level names in alphabetical order. 2.3.5. levels=c("low". Ordered factors inherit the attributes of factors.stress [1] low Levels: [1] medium high low medium high low < medium < high TRUE FALSE FALSE TRUE TRUE
> ordf.character(country). To be sure of getting the country names. it is their levels that are ordered.4. as."medium".] Printing the contents of the column with the name country gives the names.stress<-ordered(stress. and have a further ordering attribute.3.levels(islandcities$country) > lev[c(7.factor(islandcities$country.5.stress) [1] "ordered" "factor"
.58
[There are 19 cities in all. Enclose the vector of character strings in the wrapper function I() if it is to remain character."high")) > ordf. or to turn a factor into an ordered factor.stress>="medium" [1] FALSE TRUE FALSE
Later we will meet the notion of inheritance.2) > ordf. from North to South. We might prefer countries to appear in order of latitude. We can change the order of the level names to reflect this desired order:
> lev <. To extract the numeric levels 1. As in most operations with factors. you get details both of the class of the object. When you ask for the class of an object. use the function ordered(). There are though annoying exceptions that can make the use of factors tricky. 3. not the integer values.6. The levels of an ordered factor are assumed to specify positions on an ordinal scale.2. specify e. R by default turns it into a factor.e."medium". specify
as.
7. factors with ordered levels.6. Factors have the potential to cause a few surprises.level. specify
unclass(islandcities$country)
By default. this is the order that is used:
> table(islandcities$country) Australia Cuba Indonesia Japan Philippines Taiwan United Kingdom 3 1 4 6 2 1 2
This order of the level names is purely a convenience. so be careful! Here are two points to note: When a vector of character strings becomes a column of a data frame. …. R does the translation invisibly. To be sure of getting the vector of text strings. levels=lev[c(7.1)] [1] "United Kingdom" "Japan" [5] "Philippines" > table(country) United Kingdom Japan Taiwan Cuba Philippines Indonesia Australia 2 6 1 1 2 4 3 "Indonesia" "Taiwan" "Australia" "Cuba"
> country <. There are some contexts in which factors become numeric vectors.g.character(islandcities$country)
To get the integer values.1)])
In ordered factors."high"). Thus:
> class(ordf.numeric(country).level<-rep(c("low". If we form a table that has the number of times that each country appears. To create an ordered factor.2. i. specify as. there are inequalities that relate factor levels.4.

residual" "model" [5] "fitted. matrices or more general arrays.lm[1] $coefficients (Intercept) -63.214 1. or all character.786
The tenth list element documents the function call:
> elastic.matrix(1:6. enter matrix(1:6. Thus a matrix is a more restricted structure than a data frame. data = elasticband) > mode(elastic.lm$coefficients.4) created thus:
elastic. You can get the names of these objects by typing in
> names(elastic.553571
The second list element is a vector of length 7
> options(digits=3) > elastic.lm <. For this.lm$residuals 1 2.lm) [1] "coefficients" [9] "xlevels" "residuals" "call" "effects" "qr" "terms" "rank" "df.ncol=3) > xx [.1. data=elasticband)
It is readily verified that elastic. Thus consider
> xx <.1] [. scalars.lm[1]. For example. The information is preceded by $coefficients.] 1 2 3 4 5 6 # Equivalently.107 2 -0.321 3 18. and often are.] [2.values" "assign"
The first list element is:
> elastic. functions.321 7 -7.000 4 5 6 13.59
7.553571
Alternative ways to extract this first list element are:
elastic. etc.nrow=2)
.lm[“coefficients”] or elastic.lm (c.e.lm(distance~stretch. meaning “list element with name coefficients”. Note that matrices are stored columnwise. We will use for illustration the list object that R creates as output from an lm calculation. Matrices are likely to be important for those users who wish to implement new regression and multivariate methods.571429 stretch 4.571429 stretch 4.893 -27.3] [1. Lists can be.8 Matrices and Arrays
All elements of a matrix have the same mode.lm[["coefficients"]] elastic. consider the linear model (lm) object elastic. f.2] [.lm consists of a variety of different kinds of objects. sections 1. i.7 Lists
Lists make it possible to collect an arbitrary set of R objects together under a single name. One reason for numeric matrices is that they allow a variety of mathematical operations that are not available for data frames. There is a subtle difference in the result that is printed out. specify elastic. stored as a list.lm$call) [1] "call"
*7.lm[[1]]
We can alternatively ask for the sublist whose only element is the vector elastic. all numeric.1. which may have more than two dimensions. You might for example collect together vectors of several different modes and lengths.lm$call lm(formula = distance ~ stretch.4 and 2.lm$coefficients (Intercept) -63. The matrix construct generalises to array. a rag-tag of different objects.
> elastic.

Store the numbers obtained . 4) for (j in 2:length(answer)){ answer[j] <. for each of the columns of the data frame airquality (datasets package). 3.. (c) extract the vector of values of Wind for values of Ozone that are above the upper quartile. and finally nine 3s. Determine which columns of the data frame Cars93 (MASS package) are factors. but not for data frames. print out the levels vector. Generate the sequence consisting of eight 4s.] 15 18 21 24
7.answer[j-1])} answer <. Refer to the Eurasian snow data that is given in Exercise 1. and range. For each of the following calculations. 7. upper and lower quartiles. 5. Generate four repeats of the sequence of numbers (4. then seven 6s. …. for each of these data sets. 1. Write brief notes. 4) for (j in 2:length(answer)){ answer[j] <.max(answer[j]. Which of these are ordered factors? Use summary() to get information about data in the data frames attitude (both in the datasets package). extract a data frame which holds only information for small and sporty cars. 7. 102.answer[j-1])}
In the built-in data frame airquality (datasets package): (a) Determine. 12. 3).matrix() to handle any conversion that may be necessary. Use as.
. 112. 1.c(2.61
[3.2 Conversion of Numeric Data frames into Matrices
There are various manipulations that are available for matrices.sum(answer[j]. Find the mean of the snow cover (a) for the oddnumbered years and (b) for the even-numbered years. Create a vector consisting of one 1. then two 2’s. For each of these factor columns.9 Exercises
Generate the numbers 101. mean. in the columns of a 6 by 4 matrix. 5. on what this reveals. 12. From the data frame mtcars (MASS package) extract a data frame mtcars6 that holds only the information for cars with 6 cylinders. in order.
7.6 . etc. 3. the median. and store the result in the vector x. and cpus (MASS package). (b) Extract the row or rows for which Ozone has its maximum value. what you would expect? Check to see if you were right! a) b)
answer <. three 3’s. 6.c(2.8. and ending with nine 9’s. From the data frame Cars93 (MASS package).

rep(3. the test is invalid if there is clustering in the data.rep(1:5. Alternatively.
library(MASS) # if needed
To find the position at which the first space appears in the information on make of car. I note two of the simpler functions.
*8. For example. fixed=TRUE).brandnames <.x %in% c(2. This allows both a one-sample and a two-sample test.5)) > x [1] 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 > two4 <. <last position>) nchar(<vector of text strings>) ## Returns vector of number of characters in each element.2 Chi-Square tests for two-way tables
Use chisq.2 Matching and Ordering
> match(<vec1>.3.
8.62
8. <vec2>) > order(<vector>) ## For each element of <vec1>. The operator %in% can be highly useful in picking out subsets of data.test() for a test for no association between rows and columns in the output from table().
8.character(Cars93$Make). Character vectors will be sorted in alphanumeric order.1 Operations with Vectors of Text Strings – A Further Example
We will work with the column Make in the dataset Cars93 from the MASS package. we might do the following:
> car.4) > two4 [1] FALSE FALSE FALSE [11] TRUE TRUE TRUE TRUE FALSE FALSE FALSE TRUE TRUE FALSE FALSE FALSE
> # Now pick out the 2s and the 4s > x[two4] [1] 2 2 2 4 4 4
8. For example:
> x <. This test assumes that counts enter independently into the cells of a table. returns the ## position of the first occurrence in <vec2> ## Returns the vector of subscripts giving ## the order in which elements must be taken ## so that <vector> will be sorted.brandnames[1:5] [1] "Acura" "Acura" "Audi" "Audi" "BMW"
.1.sapply(strsplit(as.1 Functions for Confidence Intervals and Tests
Use the help to get complete information. the argument may be a matrix.3 String Functions
substring(<vector of text strings>. > rank(<vector>) ## Returns the ranks of the successive elements.1 The t-test and associated confidence interval
Use t.1. <first position>. Functions 8.
8. " ".test(). + function(x)x[1])
> car. Below.
Numeric vectors will be sorted in numerical order.

summary Min.Date("1/1/1960".Date() to convert text strings into dates.64
5 6 c39 c52 d21 2. by.Cars93 <.6 Merging Data Frames
The data frame Cars93 (MASS package) holds extensive information on data from 93 cars on sale in the USA in 1993. help(as. I have created a data frame Cars93.Date(c("2003/08/24".passengers Max. This can be avoided by making sure that Type has NA as one of its levels.of.y=Cars93. in both data frames.4.summary. dates are stored using January 1 1970 as origin. is Type. By default. See help(Dates).summary[.
> Cars93.47
The syntax for tapply() is similar. The default is that the year comes first. and we can proceed thus:
new.Date) and help(format.74 d21 1. If there had been rows with missing values of Type. and convert it into a dates object. the date package has been superseded by functions for working with dates that are in R base.integer() is used to convert a date into an integer value. then the month. This becomes apparent when as.y="row. Use as.
*8.passengers No. Here are examples:
> as. thus:
> # Electricity Billing Dates > dd <. The function as."2003/11/23". NA is returned. Where there are no data values for a particular combination of factor levels.9."2004/02/22".Date() will take a vector of character strings that has an appropriate format. 91 days
Use format() to set or change the way that a date is formatted. stored as a factor. One of the variables. these would have been omitted from the new data frame. and then the day of the month.as. as number %a: abbreviated weekday name (%A: unabbreviated) %m: month (00-12) %b: month abbreviated name (%B: unabbreviated) %y: final two digits of year (%Y: all four digits)
The default format is "%Y-%m-%d".0. Here however our demands are simple.
8.names")
This creates a data frame that has the abbreviations in the additional column with name “abbrev”. format="%d/%m/%Y")
.Date) for detailed information.7 Dates
Since version 1.x="Type".cars abbrev Compact Large Midsize Small Sporty Van 4 6 4 4 2 7 6 6 6 5 4 8 16 11 22 21 14 9 C L M Sm Sp V
We proceed thus to add a column that has the abbreviations to the data frame. The following are a selection of the symbols used:
%d: day.by. suitable for use in plotting.merge(x=Cars93. while a later column holds two character abbreviations of each of the car types. except that the name of the second argument is INDEX rather than by.drop=F]."2004/05/23")) > diff(dd) Time differences of 91. The output is an array with as many dimensions as there are factors. 91. in which the row names are the distinct values of Type.

Writing Functions and other Code
We have already met several functions.function(fahrenheit=32:40)(fahrenheit-32)*5/9 > # Now invoke the function > fahrenheit2celsius(c(40. just to see what it does. On the right hand side.and.
> dec1 <.format="%d:%m:%Y") [1] "1960-12-01" > as. Following the closing “)” the function body appears. The return value usually appears on the final line of the function body. You can.Date("1960-12-1")-as. Unless the result from the function is assigned to a name. SD=sd) + } > > # Now invoke the function > mean.314 1619. See help(format.444444 10.
. so that users can run the function without specifying a parameter.8.Date("1960-1-1") Time difference of 335 days > as. format="%b %d %Y") [1] "Dec 01 2004" > format(dec1.1 Syntax and Semantics
A function is created using an assignment.Date("1/1/2000".65
[1] "1960-01-01" > as. this is enclosed between curly braces ({ }). the result is printed. this was the vector consisting of the two named elements mean and sd.function(x=1:10){ + av <.Date("1/1/1970". Except where the function body consists of just one statement. In the example above.500000 3.sqrt(var(x)) + c(mean=av.50.and.Date("1:12:1960". In the example above the default was x = 1:10. the parameters appear within round brackets.integer(as.sd() mean SD 5.000000 15."%d/%m/%Y") [1] "1960-12-31" > as.151
8.sd <.integer(as. give a default.8."%d/%m/%Y") [1] 0 > as.Date("2004-12-1") > format(dec1. More generally.sd(hills$climb) mean SD 1815."%d/%m/%Y")) [1] 10957
The function format() allows control of the formatting of dates.and. format="%a %b %d %Y") [1] "Wed Dec 01 2004"
8.Date("31/12/1960". if you wish.027650 > mean.as.60)) [1] 4.mean(x) + sd <. a function returns the value of the last statement of the function. Here is a function that prints out the mean and standard deviation of a set of numbers:
> mean.555556
The function returns the value (fahrenheit-32)*5/9. Here is a function to convert Fahrenheit to Celsius:
> fahrenheit2celsius <.Date).

so that they are meaningful. where this seems appropriate.logical(match(newnames. Choose meaningful names for arguments. Settings that may need to change in later use of the function should appear as default settings for parameters. to group together parameters that belong together conceptually.function(objnames = dsetnames) { newnames <. The code needed to implement this feature has the side-effect of showing by example what the function does. even if this means that they are longer than one would like.df() that encapsulates the calculations. function(x)if(!is. nomatch = 0)) newnames[!existing] }
At some later point in the session. The following function faclev() uses is. Remember that they can be abbreviated in actual use.
> faclev <. faclev)
We can alternatively give the definition of faclev directly as the second argument of sapply.8.objects()
Now suppose that we have a function additions(). Thus dead. we may want to do similar calculations on a number of different data frames.2]. we encountered the function sapply() that can be used to repeat a calculation on all columns of a data frame. Use meaningful names for R objects. It prints out 0 if x is not a factor. It will allow us to determine whether x is a factor.4 Issues for the Writing and Use of Functions
There can be many functions. defined thus:
additions <.3 Compare Working Directory Data Sets with a Reference Set
At the beginning of a new session. Choose their names carefully. with appropriate labelling. Consider the use of names rather than numbers when you pull out individual elements. thus:
dsetnames <.8. data structures and code. Ensure that the names used reflect the hierarchies of files. Here is the definition of check.
check. Where appropriate.objects(pos=1) existing <.df <.] To apply faclev() to all columns of the data frame moths we can specify
> sapply(moths. As far as possible. make code self-documenting.as. we can enter
additions(dsetnames)
to get the names of objects that have been added since the start of the session. how many levels it has. we might store the names of the objects in the working directory in the vector dsetnames. So we create a function check."dead"] is more meaningful and safer than dead.
. objnames. function(x)if(!is.tot[. and may be useful for debugging.8. the first argument of sapply() may be a list. The built-in function is. columns etc. and otherwise the number of levels of x. and if a factor. R allows the use of names for elements of vectors and lists.factor() will return T if x is a factor.factor(x))return(0) else length(levels(x)))
Finally. provide a demonstration mode for functions. [More generally.66 8.tot[.factor() to test whether x is a factor. thus
> sapply(moths.factor(x))return(0) else length(levels(x))
Earlier.function(x)if(!is.factor(x))return(0) else length(levels(x)))
8.function(df=moths) sapply(df. and otherwise F.2 A Function that gives Data Frame Details
First we will define a function that accepts a vector x as its only argument.df(). and for rows and columns of arrays and dataframes.
8. Use lists. Such a mode will print out summary information on the data and/or on the results of manipulations prior to analysis.

The test if(!is. or whether identification by name is preferable. For this reason substantial chunks of code should be incorporated into functions sooner rather than later. given the data frame elasticband giving the amount of stretch and resulting distance of movement of a rubber band. choose file names that reflect it. and especially where you may want to retrace your steps some months later. If this number is less than . This structures the code. Use user-defined data frame attributes to document data. but where the parameter if it appears may be any one of several types of data structure. to be coerced to 1 (TRUE) or 0 (FALSE). For example. Use as the name of the list a unique and meaningful identification code. Thus if there is a factorial structure to the data files. One of R’s big pluses is its tight integration of computation and graphics. the use of switch().runif(100) correct. Each question corresponds to an independent Bernoulli trial with probability of success equal to 0. we say that the student guessed correctly. which is the same as the probability that the student guesses incorrectly. Often a good strategy is to use as defaults parameters that will serve for a demonstration run of the function. This helps ensure that functions contain welltested and well-understood components. while a 0 is recorded each time the student is wrong.7 A Simulation Example
We would like to know how well such a student might do by random guessing. This will work. Write any new “primitives” so that they can be re-used. on a multiple choice test consisting of 100 questions each with five alternatives.answers thus contains the results of the student's guesses.3) for useful functions for routine tasks.2 is exactly . Structure computations so that it is easy to retrace them. using paste() to glue the separate portions of the name together.1*(guesses < . Thus.2. Watch the r-help electronic mail list (section 13. otherwise. while the probability that a uniform random variable exceeds . because the probability that a uniform random variable is less than .2 is exactly .8. give parameters sensible defaults. Structure code to avoid multiple entry of information. A 1 is recorded each time the student correctly guesses the answer.8. We can simulate the correctness of the student for each question by generating an independent uniform random number. also. with the identification code used to determine what switch() should pick out.5 Functions as aids to Data Management
Where data.null()) then determines whether one needs to investigate that parameter further. and makes the function a source of documentation for the data.
. Concentrate in one function the task of pulling together data and labelling information. to pull out specific information and data that is required for a particular run.8.2. Lists are a useful mechanism for grouping together all data and labelling information that one may wish to bring together in a single set of computations.
8. we say that the student guessed incorrectly. NULL is a useful default where the parameter mostly is not required. labelling etc must be pulled together from a number of sources. Bear in mind. from a number of separate files. The vector correct. Wherever possible.6 Graphs
Use graphs freely to shed light both on computations and on data.2. take the same care over structuring data as over structuring code. Consider whether you should include objects as list items.
8. R can do this as follows:
guesses <.2). and Resulting Distance”
8.answers
The multiplication by 1 causes (guesses<.2) correct. which is calculated as TRUE or FALSE. You can then generate the file names automatically. perhaps with some subsequent manipulation. the uniform random number generator is simulating the student.8.67
Break functions up into a small number of sub-functions or “primitives”.answers <. We can get an idea by using simulation. Re-use existing functions wherever possible. one might specify
attributes(elasticband)$title <“Extent of stretch of band.

Apply the function to the data (in the data set huron that accompanies these notes) for the levels of Lake Huron. 9.3. with their row and column labels. We pursue the Poisson distribution in an exercise below. Write a function that generates 100 independent observations on a uniformly distributed random variable on the interval [3. 12. 3. data=milk.range(milk) par(pin=c(6. xlim=xyrange. by 6. Write a function that prints. 2. Suppose for example traffic accidents occur at an intersection with a Poisson distribution that has a mean rate of 3. Repeat for a moving average of order 3. Use it to generate 50 random integers between 0 and 99. Write a function that will take as its arguments a list of response variables.
.
plot.8. a data frame. (This means that we do not allow any number to be sampled a second time. Now look up the help for sample(). according to R. given the name of a data frame and of any two of its columns. Look up the help for the sample() function.8 Poisson Random Numbers
One can think of the Poisson distribution as the distribution of the total for occurrences of rare events. 5. Now modify the function so that you can specify an arbitrary interval. but insisting on the same range for the x and y axes. (Compare your results with the theoretical results: the expected value of a uniform random variable on [0. We leave this as an exercise. variance and standard deviation of such a uniform random variable. only those elements of a correlation matrix for which abs(correlation) >= 0. 10.75 in.
8.one <.75. Create a function to compute the average. it will plot the second named column against the first named column. Find a way of computing the moving averages in exercise 3 that does not involve the use of a for loop. showing also the line y=x.) 11. we can specify rpois(10. 7. for each patient. If no response variable is specified. 5. an accident at an intersection on any one day should be a rare event. The total number of accidents over the course of a year may well follow a distribution that is close to Poisson. Here is a function that plots.75)) abline(0.7 per year. Find the mean. Why?] The function using rpois() generates Poisson random numbers. and January 1 in the year 2000 4. For example.
8. Determine the number of days. 6.1) } # Calculates the range of all values in the data frame # Set plotting area = 6. pch=16)
Rewrite this function so that. and January 1 in the year 1800 January 1 in the year 1998. the four result against the one result.9 Exercises
1. Use the round function together with runif() to generate 100 random integers between 0 and 99. for example the mean or the median. with replacement. one sample with four units of additive and the other with one unit of additive.0833.1]. Write your own wrapper function for one-way analysis of variance that provides a side by side boxplot of the distribution of values by groups. a list of factors.7).8].7. and the variance of such a random variable is 0.) Now. one.75 in. between the following dates: January 1 in the year 1700. Seventeen people rated the sweetness of each of two samples of a milk product on a continuous scale from 1 to 7.1] is 0. # Line where four = one
plot(four. generate 50 random integers between 0 and 9. 6. sampled without replacement. variance and standard deviation of 1000 randomly generated uniform random numbers. Write a function that computes a moving average of order 2 of the values in a given vector. It will return a data frame in which each value for each combination of factor levels is summarised in a single statistic. ylim=xyrange.68
One can thus write an R function that simulates a student guessing at a True-False test consisting of some arbitrary number of questions. and use it for the same purpose. To simulate the annual number of accidents for a 10-year period.9.5. and a function such as mean or median. on [0.function(){ xyrange <. the function will generate random normal data (no difference between groups) and provide the analysis of variance and boxplot information for that. The supplied data frame milk has columns four and one. [However the total number of people injured is unlikely to follow a Poisson distribution.

Using simulation. Generate random samples from normal. Write a function that simulates the repeated calculation of the coefficient of variation (= the ratio of the mean to the standard deviation). . . Given data on trees at times 1. and so on. t (2 d. Compare with the theoretical values of . 4. f. What are the theoretical values? 16.8 accidents per year. here
1 1 1 1 1 2 2 2 3 . 18. Find the mean and variance of the student's answers. .99 of working. Write a function that calculates the minimum of a quadratic. calculates the median of the absolute values of the deviations from the sample median. 3.69
13.] t 2 " t1
[The relative growth rate may be defined as
!
!
.] *21.. . df=1)
Apply the function from exercise 17 to each sample. 15. estimate the expected value and variance of the random variable X. *20. and t (1 d. 2. . for any sample. the second is the number of occurrences of level 2. and is otherwise zero. and outputs the vector of factor levels. Write a function that.). .5 and . where each bulb has a probability . *22. Compare with the theoretical values of . plotting the result in a suitable way. and the value of the function at the minimum. Write an R function that simulates a student guessing at a True-False test consisting of 40 questions.). 14.2 and . *23. 3. where there is chance of 1 in 5 of getting the right answer to each question. 4.16. The vector x consists of the frequencies
5. df=2) b) xe<-rexp(100) d) xt2<-rt(100. *19. which is 1 if the light bulb works and 0 if the light bulb does not work. . Form a vector in which 1’s appear whenever the factor level is incremented. Assume that the number of accidents in a year follows a Poisson distribution. . A “between times” correlation matrix. for independent random samples from a normal distribution. Write a function that calculates the average of the correlations for any given lag. . 2. exponential. Compare with the standard deviation in each case. Write a function that takes any such vector x as its input.25. has been calculated from data on heights of trees at times 1. Write an R function that simulates a student guessing at a multiple choice test consisting of 40 questions. . Hence it is reasonable to calculate the average w dt dt log w 2 " log w1 over the interval from t1 to t2 as . .
[You’ll need the information that is provided by cumsum(x). 1. write a function that calculates the matrix of “average” relative growth rates over the several intervals. thus: a) xn<-rnorm(100) c) xt2<-rt(100. 6
The first element is the number of occurrences of level 1. 4. Write an R function that simulates the number of working light bulbs out of 500. 3. Write a function that does an arbitrary number n of repeated simulations of the number of accidents in a year.
1 dw d log w = . Find the mean and variance of the student's answers. f. Run the function assuming an average rate of 2. 17.

and the other extensions we describe. In the discussion of this section.. and write
!
y = g(a + b1 x1 ) + " =
exp(a + b1 x1 ) +" 1+ exp(a + b1 x1 )
Here g(.1 A Taxonomy of Extensions to the Linear Model
R allows a variety of extensions to the multiple linear regression model. logit model) y = g(a + b1x1) + Here g(.70
*9.. 1# "
We can turn this model around.. The GAM (Generalized Additive Model) model is a further extension.z p = " p (x p ) may be smoothing functions. We can add more explanatory variables: a + b1x1 + .+ " p (x p ) + #
Additive models are a generalization of lm models. Generalized Linear Models. and General Non-linear Models
GLM models are Generalized Linear Models.. in essence. The basic model formulation39 is:
Observed value = Model Prediction + Statistical Error
Often it is assumed that the statistical error values (values of in the discussion below) are independently and identically distributed as Normal. Models which have this form may be nested within other models which have this basic form. The constant term gets absorbed into one or more of the " s. The various other models we describe are... Generalized Linear Model (e. They extend the multiple regression model. !
Some of Generalized Additive Model
!
39
There are various generalizations.. and " log( ) = logit(" ) is log(odds). allow a variety of non-normal distributions.) is selected from one of a small number of options..
.
! Use glm() to fit generalized linear models.) undoes the logit transformation. In 1 dimension y = "1 (x1 ) + #
! !
z1 = "1 (x1 ).. For logit models. . and we leave until later sections the discussion of different possibilities for the “error” distribution. . where " log( ) = a + b1 x1 1# "
Here
!
!is an expected proportion. while others may be the usual linear model terms. our focus is on the form of the model prediction.z2 = " 2 (x 2 ). generalizations of this model.+
pxp +
Use lm() to fit multiple regression models. GLM.
9. In this chapter we describe the alternative functional forms.g. Multiple regression model y= +
1x1 + 2x2 +
.
Additive Model
y = "1 (x1 ) + " 2 (x 2 ) + . Thus there may be `predictions’ and `errors’ at different levels within the total model.
y = " + # . + bpxp.

1" p
The linear model predicts... We can transform to get the model y
= " p (x p ) may be smoothing functions. but log(
! !
Figure 23: The logit or log(odds) transformation. The fitting of spline (bs() or ns()) terms in a linear model or a generalized linear model can be a good ! ! alternative to the use of a full generalized additive model. !
The reason is that g(.to + .. It is not satisfactory to use a linear model to predict proportions. Some such transformation (`link’ function) as the logit is required. it is not reasonable to expect that the expected proportion will be a linear function of x.+ " p (x p )) + #
Generalized Additive Models are a generalisation of Generalized Linear Models. The function g(.) and g(. i. and log(odds) = log(
p ) = log(p) -log(1-p) 1" p p )...71
y = g("1 (x1 ) + " 2 (x 2 ) + .
9. that horse A will win the race). while others may be the
= g(z1 + z2 + . One of the commonest is the complementary log-log function. as in a logistic regression model. It is however in order to use a linear model to predict logit(proportion). we may still want to retain both "1(.99.) doing anything more that seems necessary.e...g.
The logit or log(odds) function turns expected proportions into values that may range from . Shown here is a plot of log(odds) versus proportion.+ z p ) + "
! Notice that even if p = 1. y = g("1 (x1 )) + # .z p usual linear model terms. then the corresponding odds are p/(1-p)..) does as much as it can of the task of transformation.z 2 = " 2 (x 2 ). For example. g(.. such as the inverse of the logit function.
.
!
Some of z1 = "1 (x1 ).) may be the function that undoes the logit transformation. Figure 23 shows the logit transformation. Notice how the range is stretched out at both ends.). If p is a probability (e. not p. The values from the linear model may well lie outside the range from 0 to 1. With proportions that range from less than 0.1 to 0. with "1 (.) is a specific function.2 Logistic Regression !
We will use a logistic regression model as a starting point for discussing Generalized Linear Models. A good way to think about logit models is that they work on a log(odds) scale. The logit function is an example of a link function. There are various other link functions that we can use with proportions.

so we begin by tabulating proportion that have the nomove outcome against concentration. The horizontal line estimates the expected proportion of nomoves if the concentration had no effect. the logit model and the complementary log-log model. To understand the output. Australian National University) for allowing me use of these data.1 Anesthetic Depth Example
Thirty patients were given an anesthetic agent that was maintained at a pre-determined [alveolar] concentration for 15 minutes before making an incision40.5 0 2 2
_____________________________________________ Table 1: Patients moving (0) and not moving (1).
. Alveolar Concentration Nomove 0 1 Total 0.f. We can fit the models either directly to the 0/1 data. for reasons that will emerge later.8 6 1 7 1 4 1 5 1.72 9. jerked or twisted.e.) We prefer models with a small mean residual sum of squares. Figure 24 then displays a plot of these proportions.
40
I am grateful to John Erickson (Anesthesia and Critical Care. of proportion of patients not moving. or to the proportions in Table 1.) We prefer models with a small mean deviance. A deviance has a role very similar to a sum of squares in regression.2. University of Chicago) and to Alan Welsh (Centre for Mathematics & its Applications. The interest is in estimating how the probability of jerking or twisting varies with increasing concentration of the anesthetic agent. versus concentration. Thus we have: Regression degrees of freedom sum of squares mean sum of squares (divide by d.6 0 4 4 2. i. The response is best taken as nomove.
We fit two models.4 2 4 6 1. There is a small number of concentrations. for each of six different alveolar concentrations.
Figure 24: Plot. It was then noted whether the patient moved.f. you need to know about “deviances”. Logistic regression degrees of freedom deviance mean deviance (divide by d.2 2 4 6 1.

47 conc 5.77 -0.72
(Dispersion Parameter for Binomial family taken to be 1 ) Null Deviance: 41. Here is the R code:
> anaes. they are drawn at random from some larger population. of log(odds) [= logit(proportion)] of patients not moving.981
Fig.42 2. with the same probability. With such a small sample size it is impossible to do much that is useful to check the adequacy of the model.logit) Call: glm(formula = nomove ~ conc.73
If individuals respond independently. family = binomial(link = logit).07 Coefficients: Value Std. Error t value (Intercept) -6.glm(nomove ~ conc. data = anesthetic) Deviance Residuals: Min 1Q Median 3Q Max -1. Justification for assuming the same probability will arise from the way in which individuals are sampled. versus concentration. Try also plot(anaes.68 2.0341 0.logit <. The line is the estimate of the proportion of moves.744 0. each time a new individual is taken. While individuals will certainly be different in their response the notion is that.04 -2.8 on 28 degrees of freedom Number of Fisher Scoring Iterations: 5 Correlation of Coefficients: (Intercept) conc -0. then we have Bernoulli trials.687 2. based on the fitted logit model. family = binomial(link = logit).5 on 29 degrees of freedom Residual Deviance: 27.57 2. + data = anesthetic)
The output summary is:
> summary(anaes. 25 is a graphical summary of the results:
Figure 25: Plot.logit).
.

glm)
9.1 Dewpoint Data
The data set dewpoint41 has columns mintemp. data=anesthetic)
The family parameter specifies the distribution for the dependent variable.3. i. for each combination of mintemp and maxtemp.
9. family = binomial(link = logit). maxtemp and dewpoint. but often the default works quite well. R does not at present have a facility for plots that show the contribution of each term to the model. The default link is then the identity.lm <. as for an lm model. Geography Department.logit <. data = dewpoint) options(digits=3) summary(dewpoint. We fit the model:
dewpoint = mean of dewpoint + smooth(mintemp) + smooth(maxtemp)
Taking out the mean is a computational convenience.
41
I am grateful to Dr Edward Linacre. of monthly data from a number of different times and locations. Below we give further examples.
9. You need to install the splines package.
9. One can request a `smooth’ b-spline or n-spline transformation of a column of the X matrix. The link that is commonly used here is log. Visiting Fellow. Australian National University.3. The dewpoint values are averages. data = airquality) # Assumes gaussian family. Engineering or business failures can be modelled using this same methodology.74
9.lm(dewpoint ~ bs(mintemp) + bs(maxtemp).2 Data in the form of counts
Data that are in the form of counts can often be analysed quite effectively assuming the poisson family. The log link transforms from positive numbers to numbers in the range .5 Survival Analysis
For example times at which subjects were either lost to the study or died (“failed”) may be recorded for individuals in each of several treatment groups. Also it provides a more helpful form of output.e.
. Here are details of the calculations:
dewpoint.4.to + that a linear model may predict. This way of formulating an lm type model does however have the advantage that one is not restricted to the identity link. for making these data available.3 glm models (Generalized Linear Regression Modelling)
In the above we had
anaes.3 The gaussian family
If no family is specified.
# Dataset airquality.glm<-glm(Ozone^(1/3) ~ Solar. normal errors model summary(air.4 Models that Include Smooth Spline Terms
These make it possible to fit spline and other smooth transformations of explanatory variables. One can control the smoothness of the curve. from datasets package air. The R survival package has state of the art abilities for survival analysis. In place of x one specifies bs(x)or ns(x).glm(nomove ~ conc.R + Wind + Temp. then the family is taken to be gaussian. There is an optional argument that allows us to specify the link function.lm)
9.

J.. W.8 Further Elaborations
Generalised Linear Models were developed in the 1970s. both for linear models and for a wide class of nonlinear models. 1983.. Their practical implementation built on the powerful computational abilities that. They unified a huge range of diverse methodology.
9. J. An Introduction to Statistical Modelling. N. 1990. Allow different intercepts for different habitats. A.lm() on an aov or glm object. A. Modern Applied Statistics with S-Plus. summary. and Nelder. Use log(meters) as a covariate. J. Hastie.
9. 2nd edn 1997. D. Chapman and Hall. Cambridge University Press. Generalized Additive Models.
.7 Model Summaries
Type in
methods(summary)
to get a list of the summary methods that are available. J. Each such new development builds on the theoretical and computational tools that have arisen from earlier developments. B. London. Chapman and Hall. Springer.75
9. 1989. P. You may want to mix and match. New York. by the 1970s. and Ripley. An important elaboration is to the incorporation of more than one term in the error structure. Chapman and Hall. They have now become a stock-in-trade of statistical analysts. London. McCullagh. Exciting new analysis tools will continue to appear for a long time yet. The R nlme package implements such extensions. Fit a Poisson regression model to the data in the data frame moths that Accompanies these notes. Data Analysis and Graphics Using R – An Example-Based Approach.. e. 2nd edn. Generalized Linear Models.6 Non-linear Models
You can use nls() (non-linear least squares) to obtain a least squares fit to a non-linear function.
Venables. This is fortunate.10 References
Dobson. T.
9. The output may not be what you might expect. Most professional users of R will regularly encounter data where the methodology that the data ideally demands is not yet available. Maindonald J H and Braun W J 2003.g.9 Exercises
1. and Tibshirani. Practical data analysis demands further elaborations. R. had been developed for handling linear model calculations. So be careful!
9.

2 for details of these data.")) > kiwishade.867629 6 -5.867629 6 -3. plot within block. data=kiwishade) > summary(kiwishade.0073 Correlation: (Intr) shdA2D shdD2F shadeAug2Dec -0.761621 36 56.53 0. Including Repeated Measures Models
Models have both a fixed effects structure and an error structure.1 The Kiwifruit Shading Data. The function lme has associated with it highly useful abilities for diagnostic checking and for various insightful plots. The fixed effects are block and treatment (shade). In the time series context there is usually just one realisation of the series. Multi-level Models. and units within each block/plot combination.490378 Fixed effects: yield ~ shade Value Std.Error DF t-value p-value (Intercept) 100.478639 3.4556 -125.53 shadeDec2Feb -0. The random effects are block (though making block a random effect is optional.50 shadeFeb2May -0. for purposes of comparing treatments). Version 3 of lme is broadly comparable to Proc Mixed in the widely used SAS statistical package. There is a strong link between a wide class of repeated measures models and time series models. in an inter-laboratory comparison there may be variation between laboratories.1558 shadeDec2Feb -10.0001 shadeAug2Dec 3.9831 Random effects: Formula: ~1 | block (Intercept) StdDev: 2.53 0. between observers within laboratories. from the Pinheiro and Bates nlme package.lme) Linear mixed-effects model fit by REML Data: kiwishade AIC BIC logLik 265.8.1 Multi-Level Models. Repeated Measures and Time Series 10.0015 shadeFeb2May -7.97741 0.88086 <. If we treat laboratories and observers as random.20250 1.9663 278.50520 0. which may however be observed at a large number of time points.03083 1.50 0. kiwishade$shade.76
*10. For example. handle models in which a repeated measures error structure is superimposed on a linear (lme) or non-linear (nlme) model.42833 1.867629 6 1.random=~1|block/plot. The functions lme() and nlme(). Again
Refer back to section 5.28167 1.50 Standardized Within-Group Residuals:
.62282 0.
10. and between multiple determinations made by the same observer on different samples. the only fixed effect is the mean. In the repeated measures context there may be a large number of realisations of a series that is typically quite short.1. sep=".lme<-lme(yield~shade. Here is the analysis:
> library(nlme) Loading required package: nls > kiwishade$plot<-factor(paste(kiwishade$block.019373 Formula: ~1 | plot %in% block (Intercept) Residual StdDev: 1.

it (inspection time. J. data=tinting. Data are in the data frame tinting.e. random=~1|id.871441 1 vs 2
42
Data relate to the paper: Burns.58 12. We start with a model that is likely to be more complex than we need (it has all possible interactions):
itstar. random=~1|id.1 we encountered data from an experiment that aimed to model the effects of the tinting of car windows on visual performance42.11093 0.it1. random=~1|id. Here is what we get:
> anova(itstar.146187 91.lme it1. R.method="ML")
A reasonable guess is that first order interactions may be all we need.lme it2.7288 0.074879
3 1394.138171 26.430915 2 vs 3 22.data=kiwishade) > summary(kiwishade.2112 0. Here we examine the analysis for log(it).001194
Error: Within Df Sum Sq Mean Sq F value Pr(>F) Residuals 36 438. we need to work with log(csoa) and log(it).lme 1 26 3 8 AIC BIC logLik Test L. and between individuals.method="ML")
Note that we have fitted all these models by maximum likelihood.. N. csoa (critical stimulus onset asynchrony. The authors are mainly interested in effects on side window vision.1 make it clear that. Nettlebeck.
. there is the very simple model. to get variances that are approximately homogeneous.lme<-lme(log(it)~tint*target*agegp*sex. i.926906 1. White.aov<-aov(yield~block+shade+Error(block:shade).it2. M.8.lme. So we need to specify id (identifying the individual) as a random effect.72523 18. the time required for a simple discrimination task) and age are variables.18
10.1.88105
2 17 -3.method="ML")
Finally.45036 21. allowing only for main effects:
it1.lme) Model df itstar.77022 7.e.2 The Tinting of Car Windows
In section 4. Effects of car window tinting on visual performance: a comparison of elderly and young drivers. 1999..lme<-lme(log(it)~(tint+target+agegp+sex).742883 50. Plots such as we examined in section 4. data=tinting. In this data frame.aov) Error: block:shade Df block shade Residuals 2 6 Sum Sq Mean Sq F value 172. and hence in visual recognition tasks that would be performed when looking through side windows.1176 0. The variable sex is coded 1 for males and 2 for females. and Willson. i.2:
> kiwishade.57 86. while the variable agegp is coded 1 for young people (all in their early 20s) and 2 for older participants (all in the early 70s).0065 8.84 22.e. T.35 125.lme. Ergonomics 42: 428-443. the time in milliseconds required to recognise an alphanumeric target).
it2. We have two levels of variation – within individuals (who were each tested on each combination of tint and target).78
Now fsee where these same pieces of information appeared in the analysis of variance table of section 5. data=tinting. This is so that we can do the equivalent of an analysis of variance comparison. while tint (3 levels) and target (2 levels) are ordered factors.93 Pr(>F) 4.Ratio p-value 6.17 20.51
464. i.lme<-lme(log(it)~(tint+target+agegp+sex)^2.

10.2 Time Series Models
The R stats package has a number of functions for manipulating and plotting time series, and for calculating the autocorrelation function. There are (at least) two types of method – time domain methods and frequency domain methods. In the time domain models may be conventional “short memory” models where the autocorrelation function decays quite rapidly to zero, or the relatively recently developed “long memory” time series models where the autocorrelation function decays very slowly as observations move apart in time. A characteristic of “long memory” models is that there is variation at all temporal scales. Thus in a study of wind speeds it may be possible to characterise windy days, windy weeks, windy months, windy years, windy decades, and perhaps even windy centuries. R does not yet have functions for fitting the more recently developed long memory models. The function stl() decomposes a times series into a trend and seasonal components, etc. The functions ar() (for “autoregressive” models) and associated functions, and arima0() ( “autoregressive integrated moving average models”) fit standard types of time domain short memory models. Note also the function gls() in the nlme package, which can fit relatively complex models that may have autoregressive, arima and various other types of dependence structure. The function spectrum() and related functions is designed for frequency domain or “spectral” analysis.

10.3 Exercises
1. Use the function acf() to plot the autocorrelation function of lake levels in successive years in the data set huron. Do the plots both with type=”correlation” and with type=”partial”.

*11. Advanced Programming Topics 11.1. Methods
R is an object-oriented language. Objects may have a “class”. For functions such as print(), summary(), etc., the class of the object determines what action will be taken. Thus in response to print(x), R determines the class attribute of x, if one exists. If for example the class attribute is “factor”, then the function which finally handles the printing is print.factor(). The function print.default() is used to print objects that have not been assigned a class. More generally, the class attribute of an object may be a vector of strings. If there are “ancestor” classes – parent, grandparent, . . ., these are specified in order in subsequent elements of the class vector. For example, ordered factors have the class “ordered”, which inherits from the class “factor”. Thus:
> fac<-ordered(1:3) > class(fac) [1] "ordered" "factor"

Here fac has the class “ordered”, which inherits from the parent class “factor”. The function print.ordered(), which is the function that is called when you invoke print() with an ordered factor, could be rewritten to use the fact that “ordered” inherits from “factor”, thus:
> print.ordered function (x, quote = FALSE) { if (length(x) <= 0) cat("ordered(0)\n") else NextMethod(“print”) cat("Levels: ", paste(levels(x), collapse = " < "), "\n") invisible(x) }

The system version of print.ordered() does not use print.factor(). The function print.glm() does not call print.lm(), even though glm objects inherit from lm objects. The mechanism is avaialble for use if required.

11.2 Extracting Arguments to Functions
How, inside a function, can one extract the value assigned to a parameter when the function was called? Below there is a function extract.arg(). When it is called as extract.arg(a=xx), we want it to return “xx”. When it is called as extract.arg(a=xy), we want it to return “xy”. Here is how it is done.
extract.arg <function (a) { s <- substitute(a) as.character(s) }

> extract.arg(a=xy) [1] “xy”

If the argument is a function, we may want to get at the arguments to the function. Here is how one can do it
deparse.args <function (a) { s <- substitute (a) if(mode(s) == "call"){

collapse=””))
For example:
> deparse.101:110 > y <.exp . Setting
my. the syntax is checked and it is turned into code that the R computing engine can more immediately evaluate.exp are a little different from what is printed out.exp) [1] 131 > eval(my. and want to turn it into an expression? The answer is to use the function parse().e. The expression is evaluated.exp <. function (x) paste (deparse(x). s[1].expression(mean(x+y))
stores this unevaluated expression in my.args(list(x+y.3 Parsing and Evaluation of Expressions
When R encounters an expression such as mean(x+y) or cbind(x.21:30 > my.expression(mean(x+y)) > my.expression("mean(x+y)") > eval(my. and then evaluate it
> my.parse(text="mean(x+y)") > eval(my. there are two steps: The text string is parsed and turned into an expression. i. just the arguments.exp <. print(paste(“The function is: lapply (s[-1].exp2.exp2) [1] 131
.y).txt) [1] "mean(x+y)"
What if we already have “mean(x+y)” stored in a text string.83
# the first element of a 'call' is the function called # so we don't deparse that. A text string is a text string is a text string. Note that expression(mean(x+y)) is different from expression(“mean(x+y)”). R gives you as much information as it thinks helpful.”()”. collapse = "\n")) } else stop ("argument is not a function call") } “. but indicate that the parameter is text rather than a file name. as is obvious when the expression is evaluated. unless one explicitly changes it into an expression or part of an expression. Thus
> parse(text="mean(x+y)") expression(mean(x + y))
We store the expression in my. The actual contents of my.exp2 <.txt <. Upon typing in
expression(mean(x+y))
the output is the unevaluated expression expression(mean(x+y)). Let’s see how this works in practice
> x <. foo(bar))) [1] "The function is: [[1]] [1] "x + y" [[2]] [1] "foo(bar)" list ()"
11.

parse(text=right)
.call(“data.paste("data. .leftright[2] xname <.leftright[1] right <.. "on the right of the equation")) if(length(list(. i. given without explanation.frame(".exit(detach(old. argtxt. we attach the data frame so that we can leave off the name of the data frame. sep = "")) df <.vars(expr) if(length(xname) > 1)stop(paste("There are multiple variables.") listexpr <.paste(colnames. illustrates some of the possibilities.. For example
make. collapse = ". ")".do...df() NSW ACT 1 1904 . argtxt.call is used it is only necessary to use parse() in generating the argument list.4 Plotting a mathematical expression
The following. colnames = c("NSW".all.names(list(. "ACT")) { attach(old. collapse=" & "). sep = "") expr <.)) if(nam!=xname)stop("Clash of variable names") # The part to the left of the "=" # The part to the right of the "="
expr <.df = austpop. and use only the column names
make.parse(text = exprtxt) df <. eval(listexpr)) names(df) <.paste(xname.strsplit(equation.colnames df }
To verify that the function does what it should. When do.new..df)) argtxt <.df) on. ")".df <function(old.call() makes it possible to supply the function name and the argument list in separate text strings. colnames = c("NSW".){ leftright <.".colnames df }
11. Once in the function. split = "=")[[1]] left <. It needs better error checking than it has at present:
plotcurve <function(equation = "y = sqrt(1/(1+x^2))". . type in
> make..84
Here is a function that creates a new data frame from an arbitrary set of columns of an existing data frame.))==0)assign(xname. 1:10) else { nam <.") exprtxt <. "ACT")) { attach(old.paste(colnames. collapse = ".new.exit(detach(old.df) on. .e. 3
The function do.parse(text=paste("list(".df = austpop.frame”.df <. .new.df)) argtxt <.eval(expr) names(df) <.function(old.

watch for announcements on the electronic mailing lists r-help and r-announce. Also. huron and islandcities.r-project.org/doc/bib/R-books. Cambridge University Press. look in
http://mirror.
[This has become a text book for the use of S-PLUS and R for applied statistical analysis. B. Springer. N. 4th edn 2002. The units of time for hills are however hours (not minutes) in DAAG. go to:
http://cran.3 Data Sets Referred to in these Notes
Data sets included in the image file usingR. Modern Applied Statistics with S.2 Contributed Documents and Published Literature
Dalgaard.au/pub/CRAN/
For Windows 95 etc binaries.] Venables. see
http://www. Introductory Statistics with R.86
12.1 R Packages for Windows
To get information on R packages.edu. So it pays to check the CRAN site from time to time. Data Sets from Package datasets
airquality attitude cars Islands LakeHuron
Data Sets from Package MASS
Aids2 cement Animals cpus Cars93 fgl PlantGrowth michelson Rubber mtcars
. except beams.summary dewpoint huron moths rainforest ais dolphins islandcities oddbooks seedrates anesthetic elasticband kiwishade orings tinting austpop florida leafshape possum beams hills milk primates
All of these datasets.
[This is an R-based introductory text. 2002. Venables.RData
Cars93.org/doc/bib/R-other. New York.aarnet.r-project. Appendix 1 12.]
Maindonald J H and Braun W J 2008. florida.aarnet. and a number of other statisticians.html
12. and Ripley. New packages are being added all the time.edu. Together with the ‘Complements’ it gives brief introductions to extensive packages of functions that have been written or adapted by Ripley.. P. [This is an intermediate level text.html
http://www.org
The Australian link (accessible only to users in Australia) is:
http://mirror. It assumes a fair level of statistical sophistication. Springer. Explanation is careful. but often terse. For other books and literature. W.
12.au/pub/CRAN/windows/windows-9x/
Look in the directory contrib for packages.r-project. 2nd edn. Data Analysis and Graphics Using R – An Example-Based Approach. D. NY. with a biostatistical emphasis. are in the DAAG package.