text:
This action might not be possible to undo. Are you sure you want to continue?
Using R for Data Analysis and Graphics Introduction, Code and Commentary
J H Maindonald Centre for Mathematics and Its Applications, Australian National University.
©J. H. Maindonald 2000, 2004, 2008. A licence is granted for personal study and classroom use. Redistribution in any other form is prohibited. Languages shape the way we think, and determine what we can think about (Benjamin Whorf.). This latest revision has corrected several errors. I plan, in due course, to post a new document that will largely replace this now somewhat dated document, taking more adequate account of recent changes and enhancements to the R system and its associated packages since 2002. 19 January 2008
C a m ba r v ille B e llbir d
W h ia n W h ia n B y r a n ge r y 60 65 70
C o n o n da le A lly n R iv e r 75
B ulbur in
f e m a le
m a le
le n gt h
le n gt h
ear co n ch le n gt h
Lindenmayer, D. B., Viggers, K. L., Cunningham, R. B., and Donnelly, C. F. : Morphological variation among populations of the mountain brushtail possum, trichosurus caninus Ogibly (Phalangeridae:Marsupialia). Australian Journal of Zoology 43: 449-459, 1995.
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.
4 0
3 2
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
1.3.2 Entry and/or editing of data in an editor window...................................................................................6 1.3.3 Options for read.table() ..........................................................................................................................6 1.3.4 Options for plot() and allied functions ...................................................................................................7 1.4 1.5 1.6 1.7 Further Notational Details.......................................................................................................................7 On-line Help.............................................................................................................................................7 The Loading or Attaching of Datasets......................................................................................................7 Exercises ..................................................................................................................................................8
2. An Overview of R..............................................................................................................................................9 2.1 The Uses of R ................................................................................................................................................9 2.1.1 R may be used as a calculator. ...............................................................................................................9 2.1.2 R will provide numerical or graphical summaries of data .....................................................................9 2.1.3 R has extensive graphical abilities .......................................................................................................10 2.1.4 R will handle a variety of specific analyses .........................................................................................10 2.1.5 R is an Interactive Programming Language .........................................................................................11 2.2 R Objects.....................................................................................................................................................11 *2.3 Looping .....................................................................................................................................................12 2.3.1 More on looping...................................................................................................................................12 2.4 Vectors ........................................................................................................................................................12 2.4.1 Joining (concatenating) vectors............................................................................................................13 2.4.2 Subsets of Vectors................................................................................................................................13 2.4.3 The Use of NA in Vector Subscripts....................................................................................................13 2.4.4 Factors..................................................................................................................................................14 2.5 Data Frames ...............................................................................................................................................15 2.5.1 Data frames as lists ..............................................................................................................................15 2.5.2 Inclusion of character string vectors in data frames.............................................................................15 2.5.3 Built-in data sets ..................................................................................................................................15
...........................................................................................................2 Adding lines to plots .............29 3............................................................................................................................32 4.................10 More Detailed Information ................25 3........................................................27 3........................................................3 Rugplots ....................................................................4 Normal probability plots ............9...........................1 Plot methods for other classes of object..6...................27 3...................16 2..................................28 3...............4 Scatterplot matrices.................................................................................................22 3............................................................................................................................................................................................6..................................................................30 4............................................................5.....1 Examples that Present Panels of Scatterplots – Using xyplot() ..............................................................11 Exercises ...................29 3..............1.............................................................................. lines and text ................................................3 Boxplots ...............................................................1 Applying a function to all columns of a data frame ....................................31 4..................5 Dotcharts ...........6 Common Useful Functions.....................................................................3..................................................................................................................................................2 The shape of the graph sheet....................................................................................17 2.......................................................1 Histograms and density plots ........................................................................................................7 Plotting Mathematical Symbols ........................28 3.........................................25 3....................7 Making Tables......................................23 3...................................31 4......................................21 3..........1 An Approximate Miles to Kilometers Conversion....24 3.........................17 2.................................................................18 2...................................................................6..................................2 Fine control – Parameter settings.....................................4 Identification and Location on the Figure Region ..............22 3.........................1 Plotting columns in parallel ....................................... Lattice graphics....................................................................................26 3..................................................... colour and choice of plotting symbol .............................5........................................................16 2.....................................10 References ..................................................8 Guidelines for Graphs..................................21 3...............................................................32 iv ..............................................................................................2..................................................19 2......................................................................21 3...................2 A Plotting function......18 2...................................................................................................9...................5 Plots that show the distribution of data values ..........................................................................................................................................6.....................................................................26 3.............................6 Other Useful Plotting Functions ....4.......................................................2 Adding Text in the Margin.............24 3.........................................................2 locator()..........................................................................................................................................................................2 Some further examples of lattice plots .........................................................................................................................................................................................................7.....................1 Scatterplot smoothing .1 Size.....8 The Search List .....21 3.2.................................................................................................................................................................................. Plotting..........................................4..................................................17 2...6......6.............................................................................................3 Adding points..........................................19 3.....1 identify() ....................................................28 3................................1 plot () and allied functions .................9 Exercises .1 Multiple plots on the one page ......................................................................1 Numbers of NAs in subgroups of the data .......................................5....25 3...................................................9 Functions in R .................29 3.......................................................................iv 2.28 3...................................2...........................................................................................3..24 3.18 2..22 3............................................................................
.............................................35 5..................................................................................5.................1 Extraction of Component Parts of Data frames...................................................8......1 The Model Formula in Straight Line Regression..................................................................................................................................... and Principal Components Analysis........3 An incomplete list of lattice Functions..................................................................................................1 Vectors .................................4 Decision Tree models (Tree-based models) ........................................................................52 6..............2 Manipulating Model Formulae .......................55 7..........2 Patterned Data .................................................................................................................................2...............................................................................................................................................................2 Regression Objects.....................................................................3 Pointwise confidence bounds for the fitted curve ...........................................5 Exercises ......................................................................6.............................................................................. sliced and free scales...........................1 Plant Growth Example ....................1 Multivariate EDA...........................................................................38 5......4 Spline Terms in Linear Models..............33 4..................................................53 6..........................................................................................................................................41 5........................................55 7...8..........43 5.........................................................................................................................1 Model Formulae in General .....................................................................4.............................48 5..................5..............................................................................................................................................................55 7....41 5....................................35 5..................6....................................................................1 Polynomial Terms in Linear Models....................3 Discriminant Analysis ...36 5..................................................................5...............................3.....................49 5....1 The Model Matrix ................................................................55 7......................47 *5..............43 5.........................................................................40 5.....................................................................8 aov models (Analysis of Variance)................v 4.................................................................33 4.................................................................2 Weights of Books......................................................................................................1.............................................................42 5.......................... R Data Structures ..........................................1 Subsets of Vectors...................................1 The data frame Rubber...................................10 References ..................................................3.....................................46 5................................3 Model Formulae..............47 5...........................54 *7......................................38 5........35 5....................................................................................50 6.......................56 7............................................................................................2 Cluster Analysis .................................................................. Multivariate and Tree-based Methods..2 Fixed..................................................................................... and the X Matrix.......................2 What order of polynomial? ...................................................................................................................................37 *5................................................................51 6...........3........................44 *5......................................5.............................................................5 Polynomial and Spline Regression....4 Exercises .........................................6 References ..................7 Multiple Lines – Different Regression Lines for Different Species ......................................................................................................................3 Data frames..................................4..............................................................................43 5......55 7........2 Missing Values .........54 6.....................................................................6 Using Factors in R Models ..........................2 Shading of Kiwifruit Vines ....... Linear (Multiple Regression) Models and Analysis of Variance ......38 5.....................................................................................4 Multiple Linear Regression Models ....................45 5.............................................................33 5.................................................................................................................................................................................................................9 Exercises ............51 6.............56 v .........2 Other Choices of Contrasts ............................................................................52 6...................1...............
....................57 7............................... and General Non-linear Models .....8.....................................................................................57 7...................................................56 7.57 7.....7 A Simulation Example ..................................................................................................1 The t-test and associated confidence interval.........................................................................................................8....................................8 Matrices and Arrays....................................................................................................2 Matching and Ordering ..........................67 8.........................................................2 Logistic Regression ...........66 8...............................9 Exercises ........................................................................62 8.......................... Writing Functions and other Code...63 8....................................62 8.........6 Ordered Factors...........1.................................................................................................................................63 *8...........................................................................................8.....................1 Idiosyncrasies..1 Anesthetic Depth Example..................................4......64 8...........65 8................2 A Function that gives Data Frame Details ..62 8.............................4 Application of a Function to the Columns of an Array or Data Frame ................................... Functions........71 9...........................................59 7...........................57 7.......................................................................................................59 *7..................8....................................3...................................................................................................66 8...............62 *8.................................................................................................62 8.............4 Data Entry Issues..........62 8..................68 8................4.....70 9........................................................................... GLM.74 vi .....................61 7........................................................................................5 Functions as aids to Data Management..2 Chi-Square tests for two-way tables..5 Factors and Ordered Factors ...........1.................................................................2 sapply() .................................................6 Merging Data Frames..................table()............................................................3 Compare Working Directory Data Sets with a Reference Set.................57 7...7 Lists.......4............................................................................63 *8...........................8....................................................................................................................67 8.........1 A Taxonomy of Extensions to the Linear Model ...............table()....................................................1 apply() .....................................................................................................2.........................8............................................................58 7........................................................................1 Functions for Confidence Intervals and Tests..................................................................................................5 aggregate() and tapply() ...7 Dates ...4 Issues for the Writing and Use of Functions .........................61 8......................................................................................4.......................................................................................2 Conversion of Numeric Data frames into Matrices.........................8....3 Separators when using read......................................................................................................................................72 9.................................................................................................66 8..................................................................................................................................2 Missing values when using read..........9 Exercises ................................................................................................................1 Syntax and Semantics ......................................................................................................3 glm models (Generalized Linear Regression Modelling).............................................................................................................................................................8............8 Poisson Random Numbers ..............................................................62 8.......................................64 8..........................4.........70 9............................................................................................................................vi 7.....................68 *9..............3......................65 8...........................................8..1 Arrays........................................................63 8..................................................................8........6 Graphs .....8.............................................................................................................2 Data Sets that Accompany R Packages............................................................................................................67 8............................................3 String Functions...................................................................60 7....................1 Operations with Vectors of Text Strings – A Further Example ...................................
........................86 12............................................................................................76 10..........................................................................................................82 11..................76 10............3.... Including Repeated Measures Models .....3 Parsing and Evaluation of Expressions ...1.................................................................................................................74 9....................6 ...........................................................................4.............84 11..................................86 12.................................................................................9 ................................75 9.......................1 Multi-Level Models...........................................................................................................................................87 Section 1.....................3 The gaussian family .............85 12.......... Methods.................2 Time Series Models ..................79 10..........................1.........vii 9................................................................................................................................................................................................................................................................................................................78 10...... Multi-level Models..................................................................................................................................................................................................... Repeated Measures and Time Series ...........................87 Section 2................................................3 Exercises .......................................80 10................................................1.................................2 Data in the form of counts..................80 10....................................................................76 10..........................................75 9..........6 Non-linear Models .............................................................................................................................................74 9... ............................2 Extracting Arguments to Functions ...........................................................................................................................................5 Survival Analysis.........................................86 12....................................................................................................2 Contributed Documents and Published Literature ...............1 The Kiwifruit Shading Data........................................................ Advanced Programming Topics .................4 References ..4 Plotting a mathematical expression ....7 Model Summaries......82 11...........................................................................83 11..................81 *11........................75 9......................................9 Exercises ........4 Models that Include Smooth Spline Terms......................................................................................................3 Data Sets Referred to in these Notes..................................................1 R Packages for Windows..7 ................................................................................................................74 9......................................................74 9........................................................................................ Again ...............................................................1 Dewpoint Data ................................................................................................................................................82 11................................... Appendix 1.................1..............................................................................................................................................................................3...................................................5 Searching R functions for a specified token.........................................87 vii ............................................................9 .........................................87 Section 3.............86 12............................................................................8 Further Elaborations .....................................................................................................................................................................................75 *10.............................................................................3 The Michelson Speed of Light Data ..............................87 Section 7...........10 References ...............................................................75 9.................................................................4 Answers to Selected Exercises .......2 The Tinting of Car Windows ..........................74 9.
viii viii .
1 There is a language core that uses standard forms of algebraic notation. Here are points that potential users might note: R has extensive and powerful graphics abilities. check each step with care. The action of quitting from an R session uses the function call q(). with a limited tradition of experience of the limitations and pitfalls. Books that provide a more extended commentary on the methods illustrated in these examples include Maindonald and Braun (2003). there are traps that call for special care. Neither R nor any other statistical system will give the statistical expertise needed to use sophisticated abilities. Important differences are that R has no header files." The R project enlarges on the ideas and insights that generated the S language. Adaptation of available abilities allows even greater flexibility. Section 2.5 has an example. The implementation of R uses a computing model that is based on the Scheme dialect of the LISP language. Some of the model fitting routines are leading edge. While R is as reliable as any statistical software that is available. for 32-bit versions of Microsoft Windows for Linux. more than with most systems. The R system is developing rapidly. visualize and manipulate data. The R community is widely drawn. lists and so on – as a whole. Chapters 1 and 2 indicate the range of abilities that are immediately available to novice users. Hurrah for the R development team! The Look and Feel of R R is a functional language. Because R is free. . 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. Experience with the use of R is however. at no cost. 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. If simple methods prove inadequate. It is often possible and desirable to operate on objects – vectors. on the R-help list or elsewhere. Versions of R are available. Whatever the statistical system. The citation for John Chambers' 1998 Association for Computing Machinery Software award stated that S has "forever altered how people analyze. Web addresses are given below. users have no right to expect attention. that are tightly linked with its analytic abilities.6 and 9. or to know when naïve methods are inadequate. Users who want a point and click interface should investigate the R Commander (Rcmdr package) interface. allowing the calculations such as 2+3. and especially when there is some element of complication. 1 The structure of an R program has similarities with programs that are written in C or in its successors C++ and Java. The R System R implements a dialect of the S language that was developed at AT&T Bell Laboratories by Rick Becker. likely to be an educational experience. The skills needed for the computing are not on their own enough. John Chambers and Allan Wilks. or 3^11. (There are are older versions of R that support 8. and vectors of text strings can be defined and manipulated directly. New features and abilities appear every few months.) It is available through the Comprehensive R Archive Network (CRAN). 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. there can be recourse to the huge range of more advanced abilities that R offers. Beyond this. arrays. This largely avoids the need for explicit loops. to queries. most declarations are implicit. and exposed to higher standards of scrutiny than most other systems. Simple calculations and analyses can be handled straightforwardly. from application area specialists as well as statistical specialists. graphical presentation and statistical analysis. leading to clearer code. there are no pointers. for Unix and for Macintosh OS X. Expect scepticism of the use of models that are not susceptible to some minimal form of data-based validation. most computation is handled using functions.1.
if computational efficiency demands.R-project.2. CSIRO).r-project. and of other lists that serve the R community. Note that data structure is. Novice users will notice small but occasionally important differences between the S dialect that R implements and the commercial S-PLUS implementation of S. Exposing code to the critical scrutiny of highly expert users has proved an extremely effective way to identify bugs and other inadequacies. typically. I take full responsibility for the errors that remain. Datasets are also available individually.au/~johnm/r/dsets/ Section 1. both from the University of Auckland. and provided several of the data sets and a number of the exercises. ch3-4. for copies of various versions of R.g.edu. on a Unix machine with 2GB of memory.R. to open a display file window and transfer the commands across from the that window. _________________________________________________________________________ Jeff Wood (CMIS.edu.. go to http://cran. The packages give access to up-to-date methodology from leading statistical and other researchers.anu. I am grateful. The r-help email list gives access to an informal support network that can be highly effective.au/pub/CRAN There is no official support for R. R is an "open source" system. Reported bugs are commonly fixed in the next minor-minor release.anu. Under Windows an alternative to typing the commands at the console is. written by the R Development Core Team. Readers of these notes may find it helpful to have available for reference the document: "An Introduction to R".stat. repeated smaller analyses with subsets of the total data may give insight that is not available from a single global analysis. Additionally. Australian users may wish to go directly to http://mirror. Source-code is available for inspection.aarnet. Datasets that relate to these notes Copy down the R image file usingR. The R Project The initial version of R was developed by Ross Ihaka and Robert Gentleman. hosted at http://www. as demonstrated in Section 1.ac.R.RData from http://wwwmaths. With the very large address spaces now possible. Details of the R-help list. also. go to http://wwwmaths. which will usually appear within a matter of weeks. and to elicit ideas for enhancement. and systems that are specifically designed for such manipulation may be preferable. widely drawn from different institutions worldwide. the main issue is often manipulation of data. Email archives can be searched for questions that may have been previously answered.auckland. Web Pages and Email Lists For a variety of official and contributed documentation. Like Linux. and as a result of continuing improvements in the efficiency of R's coding and memory management. and for other information. R's routines can readily process data sets that by historical standards seem large – e. etc. to various scientists named in the notes who have allowed me to use their data. .org/ Note also the Australian and New Zealand list.) using the notes as commentary. With very large datasets. a regression with 500. John Braun gave valuable help with proofreading.edu. be handled by a call to a function that is written in the C or Fortran language. Development of R is now overseen by a `core team' of about a dozen people. or for adaptation to other systems. Those who write substantial functions and (more importantly) packages will find large differences. are available from the web site for the R project at http://www. 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.org and find the nearest CRAN (Comprehensive R Archive Network) mirror site. Computerintensive components can.2 The Use of these Notes The notes are designed so that users can run the examples in the script files (ch1-2. Andreas Ruckstuhl (Technikum Winterthur Ingenieurschule.nz/r-downunder. The R language environment is designed to facilitate the development of new scientific computational tools. an even more important issue for large data sets than for small data sets.6 explains how to access the datasets. supplied with R distributions and available from CRAN sites.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.
The > indicates that R is ready for another command.exe from the relevant base directory on the nearest CRAN site.0 of R has just been started. It pays to have a separate working directory for each major project. right click|paste to place a copy on the desktop. For this demonstration I will click on my r-notes icon. and then finally go to right click|properties to set the Start in directory to be the working directory that was set up earlier. and follow instructions. choose a name that closely matches the name of the workspace directory. Then right click|copy2 to copy an existing R icon. type 2+2 and press the Enter key. the >. follow the installation instructions appropriate to the operating system. at the command line prompt. Here. Users of Microsoft Windows may wish to create a separate icon for each such working directory. Enter the name of your choice into the name field. immediately after starting up. Here is what appears on the screen: > 2+2 [1] 4 > Here the result is 4. "first requested element will follow". see the README file that is included with the R distribution. The[1] says. then left click on the copy menu item". First create the directory. Figure 1 shows the command window as it appears when version 2. for version 2. Starting Up R must be installed on your system! If it is not. For ease of remembering. Figures 1 and 2 demonstrate two of the possibilities. it. For example. choose the icon that corresponds to the project that is in hand. click on its icon to start installation.0 of R.0. 2 3 This is a shortcut for "right click. . Installation is now especially straightforward for Windows users. i. The command line prompt.1 Getting started under Windows Click on the R icon.0. For more details. Either or both of the following may be used at the user's discretion. Fig. Copy down the latest SetupR.3 1. right click|rename on the copy to rename it3. a little strangely. there is just one element. 1.e. 1: The upper left portion of the R console (command line) window. Or if there is more than one icon. In interactive use under Microsoft Windows there are several ways to input commands to R. perhaps the name itself. is an invitation to start typing in your commands. Packages that do not come with the base distribution must be downloaded and installed separately. This document mostly assumes that users will type commands into the command window.
or new text added. Both are free. to send commands to R. 3: This shows the five icons that appear when the focus is on a script file window. Highlight the commands that are intended for input to R. and Print. To get a script file window.2 Use of an Editor Script Window The screen snapshot in Figure2 shows a script file window. you will be asked for the name of a file whose contents are then displayed in the window. Fig. Look under Software|Other on the CRAN page. Under both Microsoft Windows and Linux (or Unix). 2: The focus is on an R display file window. Display file windows. go to the File menu. starting from the left: Open script. click on Open script…. a further possibility is to run R from within the emacs editor4. To load an existing file. 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. There will be a message asking whether to save the workspace image.4 For later reference. Under Unix. click on New script. Fig. which have a somewhat similar set of icons but do not allow editing. 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. are another possibility. If a new blank window is required. with the console window in the background. or that have been typed or copied into the window. and ESS which runs under emacs. Under Microsoft Windows. the standard form of input is the command line interface. . which is the middle icon of the script file editor toolbar in Figs. Click on the `Run line or selection' icon. 1. 4 This requires emacs. Return focus to console. This allows input to R of statements from a file that has been set up in advance.R. 2 and 3. The text in a script file window can be edited. or to click on the in the top right hand corner of the R window. Save script. The icons are. Run line or selection. In Figure 2 the file was firstSteps.
a data frame. Now type in austpop at the command line prompt.5 attractive options are to use either the R-WinEdt utility that is designed for use with the shareware WinEdt editor. 1. It means "is assigned to". If column headings are not included in the file. Code to plot the graph is: plot(ACT ~ Year. and total population. 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. Qld 683 873 993 SA 440 565 589 646 873 WA Tas. is a "plugin" for WinEdt. 5 4 3 8 4941 6182 1 1917 1904 1409 2 1927 2402 1727 . header=TRUE) The <. look under Software|Other on the CRAN web page." "NT" Note that in the version of the dataset that is in the DAAG package. at various times since 1917. Data frames that consist entirely of numeric data have the same form of rectangular layout as numeric matrices. in R parlance. then using this file to create a graph. The object austpop is.read. 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 on a disk in drive d: > austpop <. which is free. Use of header=TRUE causes R to use= the first line to get header information for the columns. data=austpop." "Vic. 306 392 193 211 NT ACT Aust. the argument can be omitted. at various times between 1917 and 1997. Figure 4: Population of Australian Capital Territory. displaying the object on the screen: > austpop Year NSW Vic. Here is a plot (Figure 4) of the Australian Capital Territory (ACT) population between 1917 and 1997 (Figure 4). or to use the free tinn-R editor5.txt". pch=16) We first of all remind ourselves of the column names: > names(austpop) [1] "Year" [9] "ACT" "NSW" "Aust. For links to the relevant web pages. Qld 683 873 SA 440 565 WA Tas. . .is a left diamond bracket (<) followed by a minus sign (-). for WinEdt R-WinEdt and various other editors that work with R. year has a lower case y." "Qld" "SA" "WA" "Tas. 306 392 457 502 688 879 193 211 233 257 326 375 NT ACT Aust. 5 The R-WinEdt utility.3 Year A Short R Session NSW Vic.table("d:/austpop.
141. type elasticband <.166)) 1.table() takes.edit(elasticband) Figure 5: Editor window.3 Options for read.182. data input will fail. specify na.frame(stretch=c(46. which has settings that are suitable for comma delimited (csv) files that have been generated from Excel spreadsheets.42. There are several variants of read.3. If read.2 Entry and/or editing of data in an editor window To edit the data frame elasticband in a spreadsheet-like format. and both columns will be numeric. Specify header=TRUE if there is an initial row of header names.3.1 Entry of Data at the Command Line A data frame is a rectangular array of columns of data. The default is header=FALSE.table() detects that lines in the input file have different numbers of fields.strings=". Similarly. Command alternatives to the default use of a space are sep=". optionally various parameters additional to the file name that holds the data. Note in particular read.table() The function read. This last choice makes tabs separators. It is then often useful to use the function count. The following data gives."). users can control over the choice of missing value character or characters.".48. Here we will have two columns.50.frame() can be used to input these (or other) data directly at the command line. data=austpop. showing the data frame elasticband.109. which by default is NA.csv().data. We can specify more informative axis labels. distance=c(148.44. 1.52). We will give the data frame the name elasticband: elasticband <.166. In addition users can specify the separator character or characters. 1.6 A simple way to get the plot is: > plot(ACT ~ Year. and so on.54. If the missing value character is a period (".3.fields() to report the number of fields that were identified on each separate line of the file. This plot can be improved greatly. change size of the text and of the plotting symbol.table() that differ only in having different default parameter settings.173. for each amount by which an elastic band is stretched over the end of a ruler. replace 'Year' by 'year' The option pch=16 sets the plotting character to a solid black dot. with an error message that draws attention to the discrepancy. . pch=16) # For the DAAG version." and sep="\t". 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.
.g. in order to quit from the R session we had to type q(). an alternative is to click on the help menu item.. There is also a continuation prompt.4 Options for plot() and allied functions The function plot() and related functions accept parameters that control the plotting symbol. case is significant. For file names however. with the semicolon (. Typing q on its own. Thus Austpop is different from austpop. from within the R session. and case does not distinguish file names. For the names of R objects or commands. but omit the + that will appear on the commands window if the command is typed in as we show it. Place this file in the working directory and. from the datasets or MASS packages) should then be available without need for any further action. e. without the parentheses.RData") Files that are mentioned in these notes.4 > Further Notational Details As noted earlier. used when.5 On-line Help To get a help window (under R for Windows) with a list of help topics. 6 Multiple commands may appear on the one line. plot(). On Unix systems letters that have a different case are treated as different. Note: Recall that. type: > attach("usingR. and that are not supplied with R (e. Experimentation often helps clarify the precise action of an R function.) as the separator.RData.) > apropos("matrix") (This lists all function names that include the text "matrix". displays the text of the function on the screen. Examples of their use are: > help. we often continue commands over more than one line.) The function help. Details will be given in Section 3. the continuation prompt is + In these notes.search() and apropos() can be a huge help in finding what one wants.g. and the size and colour of the plotting symbol. By default. 1. available from the author's web page. the Microsoft Windows conventions apply. the command is still not complete. Anything that follows a # on the command line is taken as comment and ignored by R.start() opens a browser window that gives access to the full range of documentation for syntax. 1.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.3.7 1. packages and functions..3. type in > help(plot) The two search functions help. type: > help() In R for Windows. Try it! 1.search("matrix") (This lists all functions whose help pages have a title or alias in which the text string "matrix" appears. the command line prompt is R commands (expressions) are typed following this prompt6. This is because q is a function. and then use key words to do a search. To get help on a specific R function. following a carriage return.
Enter the data into R. Plot snow. Input the following data. 1. blowby and total. the difference being that with attach() datasets are loaded into memory only when required for use.1).8 Users can also load (use load()) or attach (use attach()) specific files. The following ten observations. iii Use the hist() command to plot a histogram of the snow cover values. .7 Exercises 1.5 1971 12. taken during the years 1970-79. (Snow cover is in millions of square kilometers): year snow.2 i.9 1976 21.3.5 1979 9. are on October snow cover for Eurasia. on damage that had occurred in space shuttle launches prior to the disastrous launch of Jan 28 1986. for 6 launches out of 24. Repeat ii and iii after taking logarithms of snow cover.3. 3. These are the data. erosion. that were included in the pre-launch charts that were used in deciding whether to proceed with the launch. To save keystrokes.7 1975 7.9 1977 12. enter the successive years as 1970:1979] ii. from the DAAG package. The attaching of data frames will be discussed later in these notes. Plot total incidents against temperature.1 showed one way to do this. plot distance against stretch.cover versus year.0 1972 14. (Data for the 23 launches where information is available is in the data set orings.1.5 1978 14.cover 1970 6. In the data frame elasticband from section 1. with (for example) column names temperature. Distinguish between the attaching of image files and the attaching of data frames. 2. These have a similar effect.0 1974 10. iv.9 1973 10.) 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. (Refer back to Section 1. [Section 1.3.
climb.a.: 15.8660254 1.: 4.1 The Uses of R 2.0000000 # at 7.529 3rd Qu.95 (minutes) to 204.00 Median: 39. An Overview of R 2. 7 There are also versions in the DAAG package and in the Venables and Ripley MASS package.1000 # Interest on $1000.500 Median: 6.5000000 0.141593 > 2*pi*6378 #Circumference of Earth at Equator. consider the data frame hills that accompanies these notes7. Typing in summary(hills)gives summary information on these variables.Rdata is in the working directory We may for example require information on ranges of variables. Data frames are central to the way that all the more recent R routines process data. .: 300 1st Qu. with the names distance. appears on subsequent lines > 2+2 [1] 4 > sqrt(10) [1] 3. if any.6293 > > pi # R knows about pi [1] 3.: 725 Median:1000 Mean:1815 3rd Qu.: 68. think of data frames as matrices. The result.: 2. Thus the range of distances (first column) is from 2 miles to 28 miles. stores rectangular arrays in which the columns may be vectors of numbers or factors or text strings. radius is 6378 km [1] 40074. for five years 2.: 8.:204.:28. and time. As a first example.9 2.:2200 Max.2 R will provide numerical or graphical summaries of data A special class of object.: 28. For now. R evaluates and prints out the result of any expression that one types in at the command line in the console window.162278 > 2*3*4*5 [1] 120 > 1000*(1+0.6 minutes.000 1st Qu.90)*pi/180) # Convert angles to radians. while the range of times (third column) is from 15.000 climb Min.000 Mean: 7.88 3rd Qu.1 R may be used as a calculator.62 Max. compounded annually [1] 435.60 # Assumes hills.5% p. There is one column for each variable.1.95 1st Qu.075)^5 .000 Max.60. Expressions are typed following the prompt (>) on the screen.:7500 time Min. thus: > load("hills.16 > sin(c(30. in km. then take sin() [1] 0.75 Mean: 57. called a data frame.Rdata") > summary(hills) distance Min. This has three columns (variables). We will discuss graphical summaries in the next section.1. where the rows are observations and the columns are variables.
920 0.10 2. Notice that the correlations between time and distance.000 Unfortunately R was not clever enough to relabel distance as log(distance). 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). We can do all this in one step.652 0.70 1. for placing text at specified positions.890 0.920 0. and so on. Correlation: We calculate the correlation matrix for the hills data: > options(digits=3) > cor(hills) distance climb distance climb time time 1. The data are stored in the data frame elasticband (DAAG package). have reduced. .805 0.00 0. Why has this happened? Straight Line Regression: Here is a straight line regression calculation. shown in Figure 5.89 0. and time as log(time).000 Suppose we wish to calculate logarithms. The variable names are the names of columns in that data frame.000 0. for specifying tick marks and tick labels.4 R will handle a variety of specific analyses The examples that will be given are correlation and regression. for labelling axes. A helpful graphical summary for the hills data frame is the scatterplot matrix.652 1.1.700 0. type: > pairs(hills) Figure 5: Scatterplot matrix for the Scottish hill race data 2.805 1. climb as log(climb).000 0. thus: > cor(log(hills)) distance climb distance climb time time 1. In addition to plot() there are functions for adding points and lines to existing graphs.724 0.1.724 1.3 R has extensive graphical abilities The main R graphics function is plot().000 0. There are various other alternative helpful forms of graphical summary. and then calculate correlations. For this. and between time and climb.
8 80. for use in the next session in the same workspace.lm(distance~stretch. For example attach("tempscales.571 stretch 4. starting with the letter `p'8.data=elasticband) Call: lm(formula = distance ~ stretch. 26.4 84.data=elasticband)) Try it! 2. data = elasticband) Coefficients: (Intercept) -63. Some possibilities are: save. any objects that were created in the current session. can be attached.g. An alternative to ls() is objects().data.9.RData # Save into the file archive. in Section 2. 30: > celsius <.lm <. it makes sense to save the contents of the working directory from time to time.data=elasticband.RData in the working directory. Try typing q.11 > plot(distance ~ stretch. . etc.25:30 > fahrenheit <.9/5*celsius+32 > conversion <. …. This allows the retention. R offers the option of saving the workspace image.0 78. by default in the file .RData save.frame(Celsius=celsius. file="tempscales.0 2.data=elasticband) > lm(distance ~stretch.RData") save(celsius.6 82. The pattern matching conventions are those used for grep(). Type in ls() to see the names of all objects in your workspace.RData") ls(pos=2) # Check the contents of the file that has been attached The parameter pos gives the position on the search list.554 More complete information is available by typing > summary(lm(distance~stretch. They can all be operated on as data. (The search list is discussed later in this chapter. mean. fahrenheit. Careful housekeeping may be needed to distinguish between objects that are 8 Type in help(ls) and help(grep) to get details.) Important: On quitting. from the working directory or (with the path specified) from another directory. Typing the name of an object causes the printing of its contents. In a long session. pch=16) > elastic. thus making objects in the file available on request.RData") Image files.2 86. including functions and data structures. which is modelled on the Unix grep command.1. In both cases there is provision to specify a particular pattern.image(file="archive. It is also possible to save individual objects. e.2 R Objects All R entities. or collections of objects into a named image file. exist as objects.5 R is an Interactive Programming Language We calculate the Fahrenheit temperatures that correspond to Celsius temperatures 25.image() # Save contents of workspace. Fahrenheit=fahrenheit) > print(conversion) Celsius Fahrenheit 1 2 3 4 5 6 25 26 27 28 29 30 77. into the file .
0 84.4 [1] 29. and the fourth is a string vector (i.91)){answer <. Initially. j=31.0 > for (j in c(31.0 78.TRUE. Then j=51. 4.FALSE) c("Canberra".4 Vectors Examples of vectors are c(2.5.e.51. and answer is assigned the value 91 + 81 = 173..FALSE.8 [1] 27.e.2.3. to do in a complicated way what we did very simply in section 2. The first two vectors above are numeric. Finally. 51 and 91: > answer <.1.TRUE.3.0 82.0 80."Newcastle".e. The workspace image will be automatically loaded upon starting another session in that directory. and answer is assigned the value 51 + 31 = 82. There is a more straightforward way to do this calculation: > sum(c(31.j+answer} > answer [1] 173 The calculation iteratively builds up the object answer. Saving the workspace image will then save everything remains. a vector with elements of mode logical). using the successive values of j listed in the vector (31.. There are often.91)) [1] 173 Skilled R users have limited recourse to loops.91).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. the third is logical (i. or a sequence of statements that are enclosed within braces . and the contents of answer can be examined by typing in answer and pressing the Enter key. as in this and earlier examples. numeric or character11.6 [1] 28.51. *92. 2. j=91.12 to be kept and objects that will not be used again.5: > # Celsius to Fahrenheit > for (celsius in 25:30) + print(c(celsius.7.1 More on looping Here is a long-winded way to sum the three numbers 31. 9/5*celsius + 32)) [1] 25 77 [1] 26. better alternatives. 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. and answer is assigned the value 31 + 0 = 31.FALSE.51."Darwin") Vectors may have mode logical."Sydney". .FALSE.2 [1] 30 86 2.1) 3:10 # The numbers 3. Before typing q() to quit. Then the procedure ends. a vector with elements of mode character). 10 c(TRUE. i. use rm() to remove objects that are no longer required.
the meaning is: "Join these numbers together in to a vector.2 Subsets of Vectors There are two common ways to extract subsets of vectors12.c(2. 1) above is an acronym for "concatenate".c(1. be important to know the "class" of such objects.3 The Use of NA in Vector Subscripts Note that any arithmetic operation or relation that involves NA generates an NA. 7.12) > y [1] 10 15 12 > z <. which we then concatenate to form a vector z: > x <.4)] [1] 11 15 # Assign to x the values 3. > 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 <. which is NA. One can then use a vector of names to extract the elements. <=."Jeff")] John Jeff 185 183 .g. i.13 The missing value symbol. thus: > c(Andreas=178.1) > x [1] 2 3 5 2 7 1 > y <. e. 15.3. The first four compare magnitudes.c(3.3)] [1] 3 15 12 2.c(3. 5.15. can be included as an element of a vector.e. >=. NA) 11 It will. == tests for equality.12) > x[c(2. The elements that are extracted are those for which the logical value is T. 2. Set y <.4. Existing vectors may be included among the elements that are to be concatenated. >. ==.8. 3. 1. > x <.7. 12 # Extract elements (rows) 2 and 4 One can use negative numbers to omit elements: > x <.12) > x[-c(126.96.36.199. 3. Specify a vector of logical values. plot() and summary(). John=185.c(10.15. Jeff=183)[c("John". y) > z [1] 2 3 5 2 7 1 10 15 12 The concatenate function c() may also be used to join lists.5.15. This determines how the method used by such generic functions as print(). Use the function class() to determine the class of an object. later in these notes. and != tests for inequality. In the following we form vectors x and y.1 Joining (concatenating) vectors The c in c(2.2.11. Specify the numbers of the elements that are to be extracted.c(x. Thus suppose we want to extract values of x that are greater than 10.8. 11. 8. 2. NA. and !=. 0. 12 A third more subtle method is available when vectors have named elements. 2.
rep("male". use y[is. Try gender <. the levels are in alphanumeric order.4 Factors A factor is stored internally as a numeric vector with values 1. Incorrect spelling of the level names will generate an error message. Hence: > levels(gender) # Assumes gender is a factor. k. "female")) This last syntax is available both when the factor is first created. This does not select an element of y.factor(gender.factor(c(rep("female".0 leaves y unchanged. The reason is that all elements of y==NA evaluate to NA.692)) (The usage is that rep("female".692))) table(gender) gender <. For example attributes(factor(1:3)) The function levels() gives a better way to inspect factor levels. by entering: gender <. To replace all NAs by 0. The class attribute of a factor has.relevel(gender.factor(gender. It has stored with it the table: 1 2 female male Once stored as a factor.0 2. or later when one wishes to change the order of levels in an existing factor.) We can change the vector to a factor. where k is the number of levels.factor(gender."male" rows now hold missing values table(gender) rm(gender) # Remove gender 13 The attributes() function makes it possible to inspect attributes. followed by 692 2's. not surprisingly. They are crucial in the representation of categorical effects in model and graphics formulae. 3. By default. . we can create a vector of strings that that holds the values thus: gender <.factor(gender) Internally the factor gender is stored as 691 1's. An attributes table gives the 'level' for each integer value13. so that "female" precedes "male". The values "female" and "male" are the levels of the factor. the 1 is translated into "female" and the 2 into "male". If the first 691 are females and the next 692 males. To cause "male" to come before "female". and there is no assignment. and in tables. the value "factor". Consider a survey that has data on 691 females and 692 males. In most cases where the context seems to demand a character string.4.691).691). "female")) table(gender) gender <. levels=c("male". 2. the space required for storage is reduced. Factors provide a compact way to store character strings.c(rep("female". levels=c("Male". and similarly for the creation of 692 copies of "male". 691) creates 691 copies of the character string "female".na(y)] <.14 Be warned that y[y==NA] <. use gender <. 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. "female")) # Erroneous . ref="male") An alternative is gender <. rep("male". levels=c("male".
cars.4] type <. This gives a data frame with the single column Type. .1 Data frames as lists A data frame is a list15 of column vectors.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. in which different columns may have different modes.of. Large.Cars93. subscripting extracts a list.4] or Cars93.summary)) are Min.Cars93. Among the data sets in the DAAG package is Cars93. vectors.names(Cars93. No.summary[. functions. A data frame is a generalisation of a matrix. available both with read. The column abbrev could equally well be stored as a factor. Max.15 2.. vectors of character strings are by default turned into factors. or all character. height and volume for 31 Black Cherry Trees. created from information in the Cars93 data set in the Venables and Ripley MASS package.passengers. then storing it in the vector type. e.frame().summary)) Compact. Thus Cars93. takes advantage of the rectangular structure of data frames. which gives girth. 2. Cars93. which is the fourth column vector of Cars93.summary[2].summary[[4]] or Cars93.5 Data Frames Data frames are fundamental to the use of the R modelling and graphics functions.is=TRUE.Cars93.summary$abbrev type <.5. . the minimum number of passengers for cars in this category).summary.4] to extract the column vector.summary[. all numeric or all factor.2 Inclusion of character string vectors in data frames When data are input using read. In general forms of list. Here it is: > Cars93. Columns can be vectors of any mode. As noted above.of.Cars93.table() and with data. use Cars93.passengers No. all of equal length. an alternative is as. The parameter setting stringsAsFactors=TRUE.e. The use of matrix-like subscripting. and the fourth has mode character. in the datasets package. is trees.frame() function is used to create data frames.5."abbrev"] type <.3 Built-in data sets We will often use data sets that accompany one of the R packages. . usually stored as data frames. All elements of any column must however have the same mode.summary[[4]] # in the fourth list element. One such data frame.5. # Take the object that is stored 2.g. . The column names (access with names(Cars93.passengers (i.summary[1.summary[. 4]. or when the data. 2. etc.table(). and abbrev. elements that are of arbitrary type.table(). The first three columns have mode numeric. i. For read. They may be any mixture of scalars.summary[. Just as with any other list.summary Min.summary[4] is a data frame with a single column. type <.summary. Any of the following14 will pick out the fourth column of the data frame Cars93.summary. Here is summary information on this data frame 14 15 Also legal is Cars93. will if needed ensure that character strings are input without such conversion.passengers Max.e.
22) 2. and a number of other functions. and Southern. C. . The function order() places NAs last. Edith & Joy London Foundation.25 :20. :13. diff(x) has one less element than x # Sort elements into order.1 Applying a function to all columns of a data frame The function sapply() takes as arguments the the data frame. J.factor() to all columns of the supplied data frame rainforest17.:11. and Helman. with NAs last # Prints a single R object # Prints multiple objects. New South Wales. one after the other # Number of elements in a vector or of a list The functions mean().:15.20 Mean Max.] 4 56 NA NA NA NA NA NA NA NA # The final column (7) is a factor NA NA dbh wood bark root rootsk branch 16 17 The list include all packages that are in the current environment. 1990: Floristics and vegetation biomass of a forest catchment. :30.6.factor) dbh FALSE wood FALSE bark FALSE root FALSE rootsk FALSE branch species FALSE TRUE > sapply(rainforest[.30 Min. Hence:ƒ > x <.25 Height :63 Min. remove NAs.W. south coastal N. > sapply(rainforest. :76 :87 3rd Qu.-7]. Cunninghamia. 20. range) [1.:80 Volume :10.e.90 Mean Max. > order(x) [1] 1 3 2 5 4 > x[order(x)] [1] [1] 1 1 2 20 22 NA 2 20 22 > sort(x) 2. 1982: Forest biomass at Butler's Creek. 1st Qu.60 3rd Qu. sort() omits any NAs.05 Median :12. take the argument na. i.00 3rd Qu.17 :77.c(1. See also Ash.16 > summary(trees) Girth Min.20 1st Qu.40 Median :24. The following applies the function is.:19. median(). B.6 Common Useful Functions print() cat() length() mean() median() range() unique() diff() sort() order() cumsum() cumprod() rev() # reverse the order of vector elements # Gives the vector of distinct values # Replace a vector by the vector of first differences # N. Kioloa.S. NA. : 8. range().:37. then proceed with the calculation.] [2. W.:72 Median :76 Mean Max. Unpublished manuscript. 2(2): 167-182. 2.rm=T.30 Type data() to get a list of built-in data sets in the packages that have been attached16. Source: Ash. but omitting NAs # x[order(x)] orders elements of x. By default. and the function that is to be applied. is. J. 1st Qu.
2.-7].7 Making Tables table() makes a table of counts.] 4 3 8 105 2 135 0.exclude=NULL) > x <.0 4 120 56 1530 Chapter 8 has further details on the use of sapply(). na. This argument is then automatically passed to the function that is specified in the second argument position. type: > search() [1] ".factor(x. number of branches over 2cm in diameter).5. For example > range(rainforest$branch.8) > x <.rm=TRUE as a third argument to the function sapply.rm=TRUE) [1] 4 120 # Omit NAs.exclude=NULL) 2. For example: > library(lattice) # The data frame barley accompanies lattice > table(barley$year. There is an example that shows how to use it to count the number of missing values in each column of data. myrtifolia 6 0 15 1 10 12 11 10 Thus for Acacia mabellae there are 6 NAs for the variable branch (i. on whether or not the vector is a factor. fraseri Acmena smithii B.] [2. specify exclude=NULL.NA. 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.rm=TRUE) dbh wood bark root rootsk branch [1. If the vector is not a factor.GlobalEnv" "package:methods" "package:stats" .e. !is. This can be changed in the course of a session.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. Specify one vector of values (often a factor) for each table margin that is required.rm. range. na. If the vector is a factor then it is necessary to generate a new factor that includes "NA" as a level.na(rainforest$branch)) FALSE TRUE Acacia mabellae C.17 The functions mean() and range().3 24. To get a full list of these directories. Specify x <. and a number of other functions. annoyingly. then determine the range One can specify na. take the parameter na. out of a total of 16 data values. called databases. The action needed to get NAs tabulated under a separate NA category depends.factor(x) > x [1] 1 Levels: [1] 1 Levels: 5 5 NA 8 1 5 8 NA 8 1 5 8 NA > factor(x.c(1. For example: > sapply(rainforest[.8 The Search List R has a search list where it looks for objects.7. 2.
9 Functions in R We give two simple examples of R functions. xlab="Bush". detach(primates).0 62.km(175) [1] 280 # Approximate distance to Sydney.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. This is however an uncommon construction .0 6.km <. The following plots the vote for Buchanan against the vote for Bush.to.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. > library(MASS) > search() [1] ". Use the function thus > miles.9. 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. its columns can be accessed by giving their names.. 2.1 An Approximate Miles to Kilometers Conversion miles.8 52. This function can be used to attach data frames or lists (use the name.0 207. For example: > av <.g.2 # R will now know where to find Bodywt "Brainwt" Once the data frame primates has been attached. specify: > miles.to. specify detach("florida") 18 Alternatively a return value may be given using an explicit return() statement.with(primates. In technical terms. when it is no longer required. attach(florida) plot(BUSH. Note also the function with(). county by county in the state of Florida. in miles The function will do the conversion for several distances all at once. BUCHANAN. without quotes) or image (.300)) [1] 160 320 480 2.2 A Plotting function The data set florida has the votes in the 2000 election for the various US Presidential candidates.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.200. e. To convert a vector of the three distances 100.km(c(100. which is searched as required for objects that the user may specify. mean(Bodywt)) 2. ylab="Buchanan") detach(florida) # In S-PLUS. the data frame becomes a database.RData) files (the file name is placed in quotes). 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. Remember to detach an attached data frame. without further reference to the name of the data frame.to.9.
11 Exercises 1. Figure 6: Election night count of votes received. and use prod() to do the calculation in 1(c) above.0 for (j in 3:5){ answer <.ylab=yvar) mtext(side=3.florida(xvar="GORE".10 More Detailed Information Chapters 7 and 8 have a more detailed coverage of the topics in this chapter. Now apply the function to the sequence 1:100.florida(). conversion of columns of numeric data into factors) from errors in data (7. additional to those covered above. "Votes in Florida. Topics from chapter 7. xlab=xvar.g.function(xvar="BUSH".4.75.1).1.g. in the US 2000 Presidential election. at this point. Add up all the numbers from 1 to 100 in two different ways: using for and using sum. Alternatively. Look up the help for the function prod(). e.j*answer } 2.florida[.10 for (j in 3:5){ answer <. by county. For each of the following code sequences.j+answer } c) answer <. that may be important for relatively elementary uses of R include: The entry of patterned data (7. What is its action? . by county. predict the result.yvar] plot(x.2) Unexpected consequences (e. yvar="NADER") 2.florida[. in \nthe 2000 US Presidential election") } Note that the function body is enclosed in braces ({ }).3) The handling of missing values in subscripts when vectors are assigned (7. As well as plot.florida(yvar="NADER") # yvar="NADER" over-rides the default plot. i.florida <. how would you expect prod() to work? Try it! 3. line=1.e.10 for (j in 3:5){ answer <. the function allows. plot.florida(). yvar="BUCHANAN"){ x <.19 Here is a function that makes it possible to plot the figures for any pair of candidates. 2. to glance through chapters 7 and 8. parameter settings are left at their defaults. y. It may pay. Then do the computation: a) answer <. Figure 6 shows the graph produced by plot. Remember also to use R's help pages and functions.j+answer } b) answer<.xvar] y <. plot.
6. Multiply all the numbers from 1 to 50 in two different ways: using for and using prod. 5. For spheres having radii 3. The volume of a sphere of radius r is given by 4 r3/3. Use sapply() to apply the function is.5 to construct a data frame with columns radius and volume. For each of the columns that are identified as factors. 20 find the corresponding volumes and print the results out in a table. Use the technique of section 2. .factor to each column of the supplied data frame tinting. 4.ordered()]. Which columns are ordered factors? [Use is. determine the levels.1.20 4. …. 5.
Explicitly setting type = "p" causes either function to plot points. Here is a further possibility attach(austpop) plot(spline(Year. until it becomes a doublesided arror.92. the par() function does this. 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. data=austpop. The mtext() function places text in one of the margins.2 Fine control – Parameter settings The default settings of parameters. type="b") The points() function adds points to a plot. For example. points().92)) Comment on the appearance that these graphs present.21 3. mtext(). and for lines() is type = "l". The lines() function adds lines to a plot19. specify detach("austpop") 3.1 plot () and allied functions The following both plot y against x: plot(y ~ x) plot(x. type="l") detach(austpop) # Fit smooth curve through points # In S-PLUS.) Here are two further examples. Try plot((0:20)*pi/10. lines(). axis(). When it is necessary to change parameter settings for a subsequent plot. text(). 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. 3. y) # Use a formula to specify the graph # Obviously x and y must be the same length. . data=austpop. identify() etc. are often adequate.25) # character expansion increases the text and plot symbol size 25% above the default. type = "l" gives lines.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. The axis() function gives fine control over axis ticks and labels. For example. Try plot(hills) # Has the same effect as pairs(hills) 3. Plotting The functions plot(). sin((0:20)*pi/10)) plot((1:30)*0. a normal probability plot. To see some of the possibilities that R offers. attach(elasticband) # R now knows where to find distance & stretch plot(distance ~ stretch) plot(ACT ~ Year. 19 Actually these functions differ only in the default setting for the parameter type. Drag the border in towards the top border. for each numeric variable. plotting a data frame gives. The default setting for points() is type = "p". type="l") plot(ACT ~ Year. par(cex=1. making the graph sheet short and wide. enter demo(graphics) Press the Enter key to move to each new graph. ACT). form a suite that plots points. The text() function adds text at specified locations. lines and text. sin((1:30)*0. such as character size.1.
They can be assigned directly. mex=1.c("Potar monkey". pch=1) # Restore to 1 figure per page # This dataset is in the MASS package. To restore the original parameter settings at some later time. "Rhesus monkey".25) # mex=1. Graphs can be enlarged or shrunk by pointing at one corner. row. on the one page.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. Here is an example: attach(elasticband) oldpar <.25.1) plot(log(body). height (in inches) and pointsize (in 1/72 of an inch). It is the relative sizes of these parameters that matter for screen display or for incorporation into Word and similar programs. e. e."Chimp") . 20 Row names can be created in several different ways.0 62.2.22 On the first use of par() to make changes to the current device. 3.2).1).par(cex=1.1. e. use mfcol instead. holding down the left mouse button.names(primates) <.8 52. The R for Windows functions win. in a two by two layout: par(mfrow=c(2. and pulling.3 Adding points. brain) plot(sqrt(body). one after the other in a rectangular layout. which must be attached 3.par(cex=1. pch=16) attach(Animals) plot(body. then changes parameters (here cex and mex) as requested.graph() or x11() that set up the Windows screen take the parameters width (in inches). > primates Bodywt Brainwt Potar monkey Gorilla Human Rhesus monkey Chimp 10.g.log(brain)) detach(Animals) par(mfrow=c(1. so that the individual plots are rectangular rather than square. The setting of pointsize (default =12) determines character heights."Human".exit(par(oldpar)) Type in help(par) to get details of all the parameter settings that are available with par(). specify oldpar <.2 The shape of the graph sheet Often it is desirable to exercise control over the shape of the graph page.2. (brain)^0."Gorilla".g. sqrt(brain)) plot((body)^0.2 115 406 1320 179 440 Observe that the row names store labels for each row20. For a column by column layout.25 expands the margin by 25% This stores the existing settings in oldpar. oldpar <. it is often useful to store existing settings.25) on.par(cex=1. enter par(oldpar). 3.0 207. For this.0 6.g. In the example below we present four different transformations of the primates data.5) plot(distance ~ stretch) par(oldpar) detach(elasticband) # Restores the earlier settings Inside a function specify. 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. so that they can be restored later.
5.rep(3. Figure 7A would be adequate for identifying points. Try plot(1.rep(2.5).5.names is available to specify.7). It sets pch=16 to make the plot character a heavy black dot. pch=(0:6)+7) text((1:7). The parameters cex and col may be used in a similar way with text(). by number or name. pch=16. ylim=c(1. 1. cex=1:7. y=Brainwt. pch=(0:6)+14) # Plot symbols 7 to 13 # Label with symbol number # Plot symbols 14 to 20 text((1:7). pos=4) detach(primates) To place the text to the left of the points. xlim=c(0. cex=1:7.table() to input data. xlab="Body weight (kg)".7). Here is the R code for Figure 7B: plot(x=Bodywt.7). Figure 7B uses the xlab (x-axis) and ylab (y-axis) parameters to specify meaningful axis titles. paste((0:6)+14). labels=row.1350)) # Specify xlim so that there is room for the labels text(x=Bodywt. with labels on points. 7).rep(2. axes=F. labels=row. ylim=c(0.3.75. 7.names(primates). This helps make the points stand out against the labelling. colour and choice of plotting symbol For plot() and points() the parameter cex ("character expansion") controls the size. points(1:7. col=1:7.5. rep(2. added to the plot that results from the above three statements.280).5). We now show how to improve it. paste((0:6)+7). Brainwt) text(x=Bodywt. pch=0:6) text(1:7. rep(4. specify text(x=Bodywt. pos=4) points(1:7. xlim=c(1. y=Brainwt. The parameter pch controls the choice of plotting symbol. y=Brainwt.names(primates).23 attach(primates) # Needed if primates is not already attached.7). while col ("colour") controls the colour of the plotting symbol.names(primates). Figure 8B improves on Figure 8A. .1 Size. ylab="") box() points(1:7.5. pos=4) # Labels with symbol number When using read. It uses the parameter setting pos=4 to move the labelling to the right of the points. but is not a presentation quality graph. plot(Bodywt. y=Brainwt. pos=2) 3. labels=paste(0:6). col=1:7) # Do not plot points The following. ylab="Brain weight (g)". demonstrates other choices of pch. the parameter row. rep(2. type="n". Figure 7: Plot of the primates data. pos=4) # pos=4 positions text to the right of the point Figure 7A shows the result. 7). xlab="". labels=row. a column that holds the row names.
4=right)."Y". rep(1. identify() labels points. 1. xlim=c(0."B". enter view.) adds text in the margin of the current plot. rep(2. rainchars.24 Note the use of pos=4 to position the text to the right of the point (1=below. col=palette()[1:6]. The data set florida has the votes for the various Presidential candidates.colors" into its 9 characters text(1:9.5) text(10.1 identify() This function requires specification of a vector x. 0. "rainbow(7)". col=cm.ylab="") text(1:6. adj=0) } To run the function.4.3. adj=0) cmtxt <. axes=F.14).3).function(){ plot(1. Here (Figure 8) is the plot: Figure 8: Different plot symbols. colours and sizes A variety of color palettes are available. then call one or other of these functions. county by county in . 2=left.5. paste(1:6). 1:9.7)."O".9). rep(0.5.colors(9).substring("cm.6). while for locator() the default setting is n = 500.2 Adding Text in the Margin mtext(side. cex=2. 3(top) and 4.1:9) # Split "cm.5.colors". For identify() the default setting of n is the number of data points. . "cm.colours() 3. One positions the cursor at the location for which coordinates are required. unless the setting of the parameter n is reached first. line. 2.c("R".5. The sides are numbered 1(x- axis). "Default palette". adj=0) rainchars <. locator() prints out the co-ordinates of points. cmtxt.5. One positions the cursor near the point that is to be identified. ylim=c(0.5) text(10. cex=2. type="n".5.4 Identification and Location on the Figure Region Two functions are available for this purpose. xlab="". Draw the graph first. cex=3) text(10. 3=top. 3. and a vector of text strings that are available for use a labels. 3."I". and clicks the left mouse button. 1. Here is a function that displays some of the possibilities: view."G". A click with the right mouse button signifies that the identification or location task is complete. a vector y. text."V") text(1:7. 2(y-axis).colors(9)". and clicks the left mouse button. col=rainbow(7).colours <..
attach(florida) plot(BUSH. xlab="Total length". Panel C has a different choice of breakpoints. xlab="Bush". 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").5.1 Histograms and density plots The shapes of histograms depend on the placement of the breaks.5 Plots that show the distribution of data values We discuss histograms. boxplots and normal probability plots. breaks = 72. . so that a density curve can be readily superimposed. a very crude form of density plot.5. attach(possum) here <. are often a preferred alternative to a histogram. ylim = c(0. We plot the vote for Buchanan against the vote for Bush.5 + (0:5) * 5. County) detach(florida) Click to the left or right.sex == "f" hist(totlngth[here]..25 the state of Florida. superimposed. Panel B is drawn with a density scale. A histogarm is. density plots. main ="A: Breaks at 72. as Figure 9 illustrates: Figure 9: Panel A shows a histogram with a frequency scale. in fact. The density curve is. 3.5. ylab="Buchanan") identify(BUSH.. and slightly above or below a point. When labelling is terminated (click with the right mouse button). then invoking identify() so that we can label selected points on the plot. 77. 3. depending on the preferred positioning of the label. again. in order. BUCHANAN.") detach(possum) Density plots. now that they are available. the row numbers of the observations that have been labelled are printed on the screen.2 locator() Left click at the locations whose coordinates are required attach(florida) locator() detach(florida) # if not already attached plot(BUSH. 3. . BUCHANAN. so that the histogram gives a rather different impression of the distribution of the data. 22).4. xlab="Bush". BUCHANAN. Here is the code used to plot the histogram in Figure 10A.
range(dens$x) ylim <. so that the histogram gives a rather different impression of the distribution of the data. the frequency for the rectangle is divided by the width of the rectangle. xlim = xlim. x11(width=8. In order to calibrate the eye to recognise plots that indicate nonnormal variation. height=6) # This is a better shape for this plot .. This gives a scale that is appropriate for superimposing a density curve estimate. The following will give a density plot: attach(possum) plot(density(totlngth[here]). i.26 Density plots do not depend on a choice of breakpoints. ylim = ylim. side=1) detach(possum) Figure 12 adds information that is intended to assist in the interpretation of boxplots. The code for the boxplot is attach(possum) boxplot(totlngth[here].3 Boxplots We now (Figure 10) make a boxplot of the above data: Figure 10: Boxplot of female possum lengths.density(totlngth[here]) xlim <.range(dens$y) hist(totlngth[here].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. breaks = 72. The code for Figure 9B is: attach(possum) here <. 3.e. main="") lines(dens) detach(possum) 3. The choice of width and type of window.sex == "f" dens <. The only difference between Figure 9C and Figure 9B is that a different choice of breakpoints is used for the histogram. does affect the appearance of the plot. with additional labelling information. The points of this plot will lie approximately on a straight line if the distribution is Normal. The main effect is to make it more or less smooth. it is helpful to do several normal probability plots for random samples of the relevant size from a normal distribution.4 Normal probability plots qqnorm(y) gives a normal probability plot of the elements of y. xlab="Total length".5. probability = TRUE. horizontal=TRUE) rug(totlngth[here].5 + (0:5) * 5.5. controlling the nature and amount of smoothing.
3. Figure 11: Normal probability plots. The other plots are for independent normal random samples of size 43. In order to judge whether data are normally distributed. along a line. It is a way to train the eye. If data are from a normal distribution then points should fall.D. By default. sport and body-size dependency of hematology in highly trained athletes. approximately. There is one unusually small value. The idea is an important one. xlab = "Height".sex == "f" par(mfrow=c(2. main="Simulated") detach(possum) # Before continuing. type dev. 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. we use data on the heights of 100 female athletes21. The plot in the top left hand corner shows the 43 lengths of female possums.B.4)) y <. then adds a smooth curve through the points. 3. .totlngth[here] qqnorm(y.smooth(ht[here].pcBfat[here]) detach(ais) 21 Data relate to the paper: Telford. R.1 Scatterplot smoothing panel.off() # A 2 by 4 layout of plots Figure 11 shows the plots.sex=="f" plot(pcBfat[here]~ht[here]. ylab="Simulated lengths". rnorm() generates random samples from a distribution with mean 0 and standard deviation 1. Medicine and Science in Sports and Exercise 23: 788-794. ylab="Length". main="Possums") for(i in 1:7)qqnorm(rnorm(43).xlab="".6. examine a number of randomly generated samples of the same size from a normal distribution.smooth() plots points. 1991: Sex.xlab="".27 attach(possum) here <. For example: attach(ais) here<.6 Other Useful Plotting Functions For the functions demonstrated here. ylab = "% Body fat") panel. and Cunningham. R.
The parameters may be an intercept and slope.par(mar=c(4. attach(ais) here<. horizontal=TRUE) rug(ht[here].5) .1). or a vertical line (v = <ordinate>). 3. but shrinks the sizes of the points so that they almost disappear: dotchart(islands.28 3. with a rugplot added on the y-axis. ylab = "Height". ylab = "% Body fat") abline(lm(pcBfat[here] ~ ht[here])) detach(ais) 3. 0. Figure 12 shows a boxplot of the distribution of height of female athletes. mgp=c(2.5 Dotcharts These can be a good alternative to barcharts.1.1.6.sex=="f" plot(pcBfat[here] ~ ht[here]. boxwex = 0. They have a much higher information to ink ratio! Try dotchart(islands) # vector of named numeric values.1. The following is better.1. It can however be particularly useful for showing the actual values along the side of a boxplot. 3.sex == "f" oldpar <. or a vector that holds the intercept and slope.4 Scatterplot matrices Section 2.6. or an lm object.6. and there is substantial overlap. side = 1) par(oldpar) detach(ais) The parameter boxwex controls the width of the boxplot. vertical bars showing the distribution of values of x.0)) boxplot(ht[here]. Alternatively it is possible to draw a horizontal line (h = <height>).1.2 Adding lines to plots Use the function abline() for this. Figure 12: Distribution of heights of female athletes. along the x-axis of the current plot. xlab = "Height".35. 1.3 Rugplots By default rug(x) adds. Here is the code attach(ais) here <. cex=0. The bars on the left plot show actual data values.6.1. in the datasets base package Unfortunately there are many names.75.3 demonstrated the use of the pairs() function.25.
29 3. In scatterplots the graph should attract the eye's attention to the points that are plotted. Use scatterplots in preference to bar or related graphs whenever the horizontal axis represents a quantitative effect. Draw graphs so that reduction and reproduction will not interfere with visual clarity. 95% confidence interval. The data set huron that accompanies these notes has mean July average water surface elevations.8 Guidelines for Graphs Design graphs to make their point tersely and clearly. When there is extensive overlap of plotting symbols. for example analytical error when the relevant source of `error' for comparison of treatments is between fruit. See help(plotmath) for detailed information on the plotting of mathematical expressions. Use graphs from which information can be read directly and easily in preference to those that rely on visual impression and perspective. In plot(). in feet. Use colour or different plotting symbols to distinguish different groups. y. 3. 3. xlab="Radius". Figure 13 was produced with x <. there is a double equals sign ("=="). or whatever. ylab=expression(Area == pi*r^2)) Figure 13: The y-axis label is a mathematical expression. Explain clearly how error bars should be interpreted — SE limits. use open plotting symbols.1:15 plot(x. SD limits. which is exactly what one wants. Notice that in expression(Area == pi*r^2). with a single equals sign. . overlapping points will give a high ink density. Thus in scientific papers contour plots are much preferable to surface plots or twodimensional bar graphs. with a minimum waste of ink. It is pointless to present information on a source of error that is of little or no interest. Where points are dense. either or both of xlab and ylab can be an expression. Use solid points. Take care to use colours that contrast. and to important grouping in the data. Explain what source of `error(s)' is represented. when there is little or no overlap. Label as necessary to identify important features. There is a further example in chapter 12. The final plot from demo(graphics) shows some of the possibilities for plotting mathematical symbols. although what will appear on the plot is Area = pi*r^2.9 Exercises 1. large enough to stand out relative to other features.7 Plotting Mathematical Symbols Both text() and mtext() will take an expression rather than a text string.
rt(10. Australian Journal of Zoology 43: 449-458.2) for four separate `random' samples of size 10. U. 3. National Oceanic and AtmosphericAdministration. [To access this data frame. Viggers. C. all from a normal distribution. you might like to try it for some further distributions. For example x <rchisq(10.math. the mean levels are correlated from year to year. Then repeat exercise 6 above.4. that has mean heights for 1875-1772 only. 8. D. In the bottom four rows. For the first two columns of the data frame hills. display plots for samples of size 100. Dept. 22 Source: Great Lakes Water Levels... Try x <. specify library(MASS)] 3. If you find exercise 6 interesting. . 1995. Use the row labels to label the points with the three largest body weight values. National Ocean Survey. examine the distribution using: (a) histogram. Try x <.height against year.1) will generate 10 random values from a t distribution with degrees of freedom one. To see how each year's mean level is related to the previous year's mean level.rnorm(10). That is. of Commerce.html#DataVis has links to many different collections of information on statistical graphics. now working with the logarithms of the data values.huron$mean. use lag. 6.height.runif(10). What are the advantages and disadvantages of these different forms of display? 4. and d) a density plot (plot(density(possum$hdlngth)). press both mouse buttons simultaneously. Plot the graph of brain weight (brain) versus body weight (body) for the data set Animals from the MASS package. The statement x <. (Alternatively work with the vector LakeHuron from the datasets package. Trichosurus caninus Ogilby (Phalangeridae: Marsupialia). Label the axes in untransformed units.30 IGLD (1955) for Harbor Beach.yorku.plot(huron$mean. 5. Now generate a sample of size 10 from a normal distribution with mean 170 and standard deviation 4. B.ca/SCS/StatResource. What shape do the points follow? *7. but taking samples from a uniform distribution rather than from a normal distribution. R. and Donnelly. 2. type identify(huron$year. In the middle 4 rows. c) a normal probability plot (qqnorm(possum$hdlngth)). b) a stem and leaf plot (stem(qqnorm(possum$hdlngth)). 1860-1986.labels=huron$year) and use the left mouse button to click on the lowest point and highest point on the plot. F. log(brain weight) versus log(body weight). Michigan. Use mfrow() to set up the layout for a 3 by 4 array of plots. for 1860-198622. K.height) This plots the mean level at year i against the mean level at year i-1. show normal probability plots (section 3. (c) normal probability plots. c) As in the case of many time series. Repeat (a). Alongside on the same graphics page. B. Make normal probability plots for samples of various sizes from these distributions. Cunningham. by default on the interval 0 to 1. (b) density plots. Station 5014. Check the distributions of head lengths (hdlngth) in the possum23 data set that accompanies these notes. display plots for samples of size 1000. Comment on how the appearance of the plots changes as the sample size changes. and print out the numbers you get. Label the axes appropriately.10 References The web page http://www.1) will generate 10 random values from a chi-squared distribution with degrees of freedom 1. The function runif() generates a sample from a uniform distribution. b) Use the identify function to label points for the years that correspond to the lowest and highest mean levels.) a) Plot mean.. To quit. Print out the numbers that you get.S. Compare the following forms of display: a) a histogram (hist(possum$hdlngth)). L. 23 Data relate to the paper: Lindenmayer. In the top 4 rows. Morphological variation among populations of the mountain brush tail possum. Look up the help for rnorm. on Lake Huron. (b) and (c).
1 Examples that Present Panels of Scatterplots – Using xyplot() The basic function for drawing panels of scatterplots is xyplot(). and Willson. Figure 14 shows the style of graph that one can get from xyplot(). The older coplot() function that is in the base package has some of same abilities as xyplot( ). the time in milliseconds required to recognise an alphanumeric target). both these packages must be installed. Nettlebeck. Figure 14: csoa versus it. the time required for a simple discrimination task) and age are variables. Lattice graphics Lattice plots allow the use of the layout on the page to reflect meaningful aspects of data structure. N. i. for each combination of females/males and elderly/young. . data=tinting) # Simple use of xyplot() Here is the statement used to get Figure 16. Providing it is installed. The different symbols are different contrasts. J. it would be enough to type: xyplot(csoa ~ it | sex * agegp. 2 = an older participant. To use lattice graphics. csoa (critical stimulus onset asynchrony.. 4. R. Ergonomics 42: 428-443. In this data frame. The authors were mainly interested in visual recognition tasks that would be performed when looking through side windows. 2 = female) and agegp (1 = young. in the early 70s) are factors..e. e. hicon) are ordered factors. tint (level of tinting: no. we might plot csoa against it for each combination of sex and agegp. the grid package will be loaded automatically when lattice is loaded. in the early 20s. They offer abilities similar to those in the S-PLUS trellis library. Effects of car window tinting on visual performance: a comparison of elderly and young drivers. T. In a simplified version of Figure 16 above. We will use the data frame tinting (supplied) to demonstrate the use of xyplot(). For this simplified version. sex (1 = male. 1999. The two different symbols distinguish between low contrast and high contrast targets. M. + = high contrast) are shown with different symbols. but with a limitation to two conditioning factors or variables only. hi) and target (contrast: locon. i. 24 Data relate to the paper: Burns. 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. lo.31 4. The two targets (low. White. it (inspection time.
groups=tint. outer=FALSE."smooth") The relationship between csoa and it seems much the same for both levels of contrast.key=list(columns=2)) If colour is available. cex=2) ) .settings=simpleTheme(pch=16. data=tinting. +=low. It is apparent that the long response times for some of the elderly males occur.32 xyplot(csoa~it|sex*agegp. separately for the two target types: xyplot(csoa~it|sex*agegp. groups=target. The longest times are for the high level of tinting. A striking feature is that the very high values. The following puts smooth curves through the data. The following use the dataset grog from the DAAGxtras package: library(DAAGxtras) xyplot(Beer+Spirit+Wine ~ Year | Country. 4. type=c("p". data=tinting. If on the same panel.2. with the low contrast target.2 Some further examples of lattice plots These are given with a minimum of explanation. different colours will be used for the different groups.key=list(columns=3). as we might have expected. 4. it is desirable to use auto. groups=target. outer=TRUE. for both csoa and it. Finally. groups=Country. we do a plot (Figure 15) that uses different symbols (in black and white) for different levels of tinting. data=tinting. data=grog) xyplot(Beer+Spirit+Wine ~ Year. data=grog. auto. data=grog) xyplot(Beer+Spirit+Wine ~ Year | Country.key to give a simple key. outer=TRUE. The different levels of tinting (no. par. >=high) are shown with different symbols.1 Plotting columns in parallel Use the parameter outer to control whether the columns appear on the same or separate panels. panel=panel. auto.superpose. occur only for elderly males. xyplot(csoa~it|sex*agegp. for each combination of females/males and elderly/young. auto.key=list(columns=3)) Figure 15: csoa versus it.
930. Investigate the effect of varying the density bandwidth (bw). 961 Present the data (i) in side by side boxplots (use bwplot()).33 In the final plot.) qqmath(numeric ~ numeric . 973. . 2. . generate the following plots: a) histograms of hdlngth – use histogram(). 945.different slices of same scale xyplot(BC+Alberta ~ Date.) dotplot(factor ~ numeric .) In each instance. 899. 4. data=jobs.3 An incomplete list of lattice Functions splom( ~ data.frame) bwplot(factor ~ numeric . . 25 Data are from the paper ``Smoking During Pregnancy and Lactation and Its Effects on Breast Milk Volume'' (Amer. 981.) qqnorm(numeric .2.. . . data=jobs. unless reset. Use of the parameter par. .) barchart(character ~ numeric . 1086. In most cases. . Display contourplot(numeric ~ numeric * numeric. . . .) parallel( ~ dataframe. the plot is repeated for the groupings within the one panel. 1202.set() to make the changes for the duration of the current device. .) histogram( ~ numeric . . 655. 693 Nonsmoking Mothers: 947. .) stripplot(factor ~ numeric . . 767. .) # Smoothed version of histogram # QQ plot # Scatterplot matrix # Parallel coordinate plots # 3-D plot # Contour plot # Variation on a contour plot # Scatterplot matrix # Box and whisker plot # normal probability plots # 1-dim. a groups parameter can be specified. sliced and free scales library(DAAG) ## scale="fixed" xyplot(BC+Alberta ~ Date. . i. conditioning variables can be added. (ii) using a dotplot form of display (use dotplot()). . 895. The following exercises relate to the data frame possum that accompanies these notes: (a) Using the xyplot function. .independent scales xyplot(BC+Alberta ~ Date. . outer=TRUE. . taking into account sex and Pop. scales=list(y=list(relation="free")) ) 4. 745. explore the relation between hdlngth and totlngth. 4. Use trellis. . c) density plots of hdlngth – use densityplot(). . . 793. 753. 714.) cloud(numeric ~ numeric * numeric.e. 545. outer=TRUE) ## scale="sliced" . . 593. scales=list(y=list(relation="sliced")) ) ## scale="free" . 12. Display # 1-dim.4 Exercises 1.par.) splom( ~ dataframe. data=jobs.2 Fixed. note the use of simpleTheme() as a simple mechanism for controlling common parameter settings. . . b) normal probability plots of hdlngth – use qqmath(). The following data gives milk volume (g/day) for smoking and nonsmoking mothers25: Smoking Mothers: 621.) densityplot( ~ numeric . For the possum data set. 903.settings makes the change for the current plot only. 598.) levelplot(numeric ~ numeric * numeric. of Clinical Nutrition). outer=TRUE. . J.
(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.
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)
The code for the regression calculation is:
elastic.lm <- lm(distance ~ stretch, data=elasticband)
Here distance ~ stretch is a model formula. Other model formulae will appear in the course of this chapter. Figure 16 shows the plot:
Figure 16: Plot of distance versus stretch for the elastic band data, with fitted least squares line The output from the regression is an lm object, which we have called elastic.lm . Now examine a summary of the regression results. Notice that the output documents the model formula that was used:
> options(digits=4) > summary(elastic.lm) Call: lm(formula = distance ~ stretch, data = elasticband) Residuals: 1 2.107 2 -0.321 3 18.000 4 5 6 13.321 7 -7.214 1.893 -27.786
Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) stretch -63.57 4.55 74.33 1.54 -0.86 2.95 0.431 0.032
Residual standard error: 16.3 on 5 degrees of freedom Multiple R-Squared: 0.635, Adjusted R-squared: 0.562 p-value: 0.0319 F-statistic: 8.71 on 1 and 5 degrees of freedom,
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"
[5] "fitted.values" "assign" [9] "xlevels" "call" "qr" "terms" "df.residual" "model"
Various functions are available for extracting information that you might want from the list. This is better than manipulating the list directly. Examples are:
> coef(elastic.lm) (Intercept) -63.571 1 2.1071 2 -0.3214 stretch 4.554 3 18.0000 4 5 6 13.3214 7 -7.2143
> resid(elastic.lm) 1.8929 -27.7857
The function most often used to inspect regression output is summary(). It extracts the information that users are most likely to want. For example, in section 5.1, we had
There is a plot method for lm objects that gives the diagnostic information shown in Figure 17.
Figure 17: Diagnostic plot of lm object, obtained by plot(elastic.lm). To get Figure 17, type:
x11(width=7, height=2, pointsize=10) par(mfrow = c(1, 4), mar=c(5.1,4.1,2.1,1.1)) plot(elastic.lm) par(mfrow=c(1,1))
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
etc.6 + 4. The aim is to reproduce.1 173-154. attr(. the first column by a and the second column by b.3. aov.37 The parameters that are to be estimated are a and b.matrix() prints out the model matrix. etc.55 … Stretch 46 = 145. fac:x allows a different slope for each different level of fac. and adding.6) and b ( = 4.lm) The following are the fitted values and residuals that we get with the estimates of a (= -63. data=elasticband) (Intercept) stretch 1 2 3 4 .1 … ! 46 54 48 50 … Note the use of the symbol ˆ y [pronounced y-hat] for predicted values. .8 = 18.6 + 4.6 1 1 1 1 … 4.1 48 = 154. Fitted values are given by multiplying each column of the model matrix by its corresponding parameter. chooses a and b so that the sum of squares of the residuals is as small as possible. y~x + fac + fac:x : lm.) Model formulae are widely used to set up most of the model calculations in R.8 50 = 163.6 + 4. . Another name is predicted values. For this we have y=x b+e ! where e is a random variable with mean 0. The X matrix then consists of a single column.3 182-182. the values in the y-column.6 + 4. glm. 5.55 -63. the x's. i.1.9 y (Observed) Distance (mm) ! 148 182 173 166 … ˆ y"y (Residual) Observed – Fitted 148-145. which is the form of regression that we describe in this course. is: model.7 = 2. as closely as possible. Notice the similarity between model formulae and the formulae that are used for specifying coplots."assign") [1] 0 1 1 1 1 1 46 54 48 50 Another possibility. Thus.55 -63. The function model. glm.matrix(distance ~ stretch.1 Model Formulae in General Model formulae take a form such as: y~x+z : lm.7 54 = 182. 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 . (If fac is a factor and x is a variable.e. with elastic. We might alternatively fit the simpler (no intercept) model.55 -63.lm as in section 5. Thus: > model.55 ˆ y (Fitted) -63.. The residuals are the differences between the values in the y-column and the fitted values.2 166-163.matrix(elastic.6 + 4. Least squares regression.55 -63.55) that result from least squares regression X Stretch (mm) -63.9 = 2.1 = -0.
The variables are loss (the abrasion loss in gm/hr).1395 Note that graphics formulae can be manipulated in exactly the same way as model formulae.1 The data frame Rubber The data set Rubber from the MASS package is from the accelerated testing of tyre rubber26. Table 6. 5.3. Davies (1947) Statistical Methods in Research and Production.g. be a text string.formula("y~x+z") The argument to formula() can. as just demonstrated."~". formyxz <. e.3780 distance 0.formula(paste(nam[1]. hard (hardness in `Shore' units). . For example > names(elasticband) [1] "stretch" "distance" > nam <.nam[2])) > lm(formds.4. data=elasticband) Call: lm(formula = formds.names(elasticband) > formds <. data = elasticband) Coefficients: (Intercept) 26. Oliver and Boyd. 26 The original source is O.formula(y~x+z) or formyxz <. We first obtain a scatterplot matrix (Figure 18) : Figure 18: Scatterplot matrix for the Rubber data frame from the MASS package.L.2 Manipulating Model Formulae Model formulae can be assigned. This makes it straightforward to paste the argument together from components that are stored in text strings. and tens (tensile strength in kg/sq m). 119.4 Multiple Linear Regression Models 5.1 p.38 *5.
at the upper end of the range. data = Rubber) Residuals: Min 1Q Median 3.1)) This plot raises interesting questions.82 3Q 19.77e-011 71 on 2 and 27 degrees of freedom.61 Coefficients: Estimate Std. F-statistic: Adjusted R-squared: 0. obtained with termplot(). We proceed to regress loss on hard and tens. showing the contribution of each of the two terms in the model.583 0.84. Error t value Pr(>|t|) (Intercept) hard tens 885. A smooth curve has.smooth) par(mfrow=c(1.75 Max 65. Rubber.lm) Call: lm(formula = loss ~ hard + tens.571 -1. been fitted through the partial residuals. note the use of termplot().07 3.27 -7. data=Rubber) > options(digits=3) > summary(Rubber.38 -14. the response is not linear with tensile strength.3e-07 Residual standard error: 36. .752 0. Figure 19: Plot.8e-14 1. Figure 19) used the following code: par(mfrow=c(1. There is a clear suggestion that.828 p-value: 1.194 14.lm. at the mean of the contributions for the other term.2)) termplot(Rubber.0e-11 1.98 -79.lm(). smooth=panel.33 -11.5 on 27 degrees of freedom Multiple R-Squared: 0. partial=TRUE. In addition to the use of plot.lm(loss~hard+tens.374 61.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. in each panel.lm <.161 -6.
5 15.877 3. The design of the data collection really is important for the interpretation of coefficients from a regression equation. Error t value Pr(>|t|) (Intercept) thick height -1. Even though regression equations from observational data may work quite well for predictive purposes.2 21.5 23.755 0. the individual coefficients may be misleading.552 0.9 6.219 13.117 0. 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.263 0.719 0.data=logbooks) > logbooks.5 12.7 8.26e-04 # proportional to thick * height * width > logbooks.7 19.0 20.3 22.07 0.829 0.678 -0.lm1<-lm(weight~thick. This is more than an academic issue.lm3<-lm(weight~thick+height+width.316 0.070 0.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. The regression estimate of the treatment effect.0124 > logbooks. height and width are so strong that if one tries to use more than one of them as a explanatory variables. as the analyses in Lalonde (1986) demonstrate. the coefficients are ill-determined.216 0.708 0.lm3)$coef Estimate Std.5243 0.40 5.6 10.8 13.lm2)$coef Estimate Std. The regressions on height and width give plausible results.5 23.6 23.5 19.434 1.69 -1.4.070 -0. weight reduces. while the coefficient of the regression on thick is entirely an artefact of the way that the books were selected.117 So is weight proportional to thick * height * width? The correlations between thick.1 27.8 17.lm1)$coef Estimate Std. Error t value Pr(>|t|) (Intercept) thick 9.6 12.5 29.data=logbooks) > summary(logbooks. Error t value Pr(>|t|) (Intercept) thick height width -0.35e-08 -4.5 18.lm2<-lm(weight~thick+height. They contain very similar information.2 1075 940 625 400 550 600 450 450 300 690 400 250 Notice that as thickness increases. They had data from experimental "treatment" and "control" groups.2 14.7303 0. .0 9.356 0.907 0.465 0.313 2.114 3. as is evident from the scatterplot matrix. 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.662 3.log(oddbooks) # We might expect weight to be > > summary(logbooks.273 1.8 17.4 11.154 1. and also from two comparable non-experimental "controls".472 0. > logbooks <.121 1.0 15. when comparison was with one of the non-experimental controls.224 1.5 15.data=logbooks) > summary(logbooks.
80 -6. with fitted quadratic curve: Figure 20: Number of grain per head versus seeding rate. Spline curves. plot(grain ~ rate.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.41 Dehejia and Wahba demonstrate the use of scores ("propensities") to identify subsets that are defensibly comparable. with fitted quadratic curve.09429 -0.1 Polynomial Terms in Linear Models The data frame seedrates27 (DAAG package) gvies. 5. New Zealand Journal of Agriculture 10: 133-136.066686 0. C. It is impossible to be sure that any method is giving the right answer.000171 0. For this. constructed by joining low order polynomial curves (typically cubics) in such a way that the slope changes smoothly. Summary details are in Maindonald.02138 0.000049 52.5.00036 0. Note that polynomial curves of high degree are in general unsatisfactory. the number of barley grain per head. for the barley seeding rate data.07294 0.50 0. for each of a number of different seeding rates.12286 Coefficients: Estimate Std. > seedrates. C.lm2 <.009911 0.lm2) Call: lm(formula = grain ~ rate + I(rate^2). data=seedrates) # Plot the data Figure 20 shows the data. (1992). Error t value Pr(>|t|) (Intercept) 24.455694 0. The propensity is then the only covariate in the equation that estimates the treatment effect. are in general preferable. . data=seedrates) > summary(seedrates.01429 27 Data are from McLeod.73 3. 5. We will need an X-matrix with a column of ones.lm(grain ~ rate+I(rate^2). Curves are constructed from linear combinations of transformed values. (1982) Effect of rates of seeding on barley grown for grain.04571 -0.00286 -0. both rate and I(rate^2) must be included in the model formula.060000 rate I(rate^2) -0. a column of values of rate. J. H. and a column of values of rate2.
2 What order of polynomial? A polynomial of degree 2.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(."assign") This was a (small) extension of linear models.g. However the paper from which these data come gives an independent estimate of the error mean square (0.matrix(grain~rate+I(rate^2). The increase in the number of degrees of freedom more than compensates for the reduction in the F-statistic.) based on 8 replicate results that were averaged to give each value for number of grains per head. This may seem good. How does one check? One way is to fit polynomials.lm2) > lines(spline(seedrates$rate. The estimate is 0.lm2. 35) [1] 0. an x3 column might be added. Type: > model. (R2 is larger when there is more variation to be explained. of each of degrees 1 and 2. . not when R2 achieves some magic threshold. Once the model matrix has been formed. e.data=seedrates) > anova(seedrates. Adjusted R-squared: 0. to handle a specific form of non-linear relationship. F-statistic: > hat <.026286 2 3 0.17 on 35 d. but given the accuracy of these data it was not good enough! The statistic is an inadequate guide to whether a model is adequate. > 1-pf(0.992 p-value: 0.187000 -1 -0.07294 The F-value is large. 5. Thus. on 35 df.) A predictive model is adequate when the standard errors of predicted values are acceptably small. Even for any one context. a quadratic curve.0244 Finally note that R2 was 0.160714 12. i. and we have p=0. but on this evidence there are too few degrees of freedom to make a totally convincing case for preferring a quadratic to a line.f.data=seedrates) > seedrates.5.4 on 1 and 35 degrees of freedom.16/0.17^2.predict(seedrates. 1.228 0. check the form of the model matrix. the F-value is now 5. Any transformation can be used to form columns of the model matrix.Sum Sq Df Sum Sq F value Pr(>F) 1 2 0.lm1<-lm(grain~rate. we are limited to taking linear combinations of columns.115 on 2 degrees of freedom Multiple R-Squared: 0.972 for the straight line model.Df Res.lm1) Analysis of Variance Table Model 1: grain ~ rate + I(rate^2) Model 2: grain ~ rate Res. and compare them thus: > seedrates. R2 will in general increase as the range of the values of the dependent variable increases.1607. Again. > # However we have an independent estimate of the error mean > # square.42 Residual standard error: 0.lm2<-lm(grain~rate+I(rate^2).172 (35 df).e.024 . on 1 df) with a mean square of 0. looked about right for the above data.0039 256 on 2 and 2 degrees of freedom.996. If we compare the change in the sum of squares (0. hat)) > # Placing the spline fit through the fitted points allows a smooth curve. > # For this to work the values of seedrates$rate must be ordered.seedrates.
Spline functions variables extend this idea further. 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. once the knots are fixed..lty=2) > lines(cars$speed. Different choices.fit=TRUE) > names(ci5) [1] "fit" "se.data.lm(grain ~ rate + I(rate^2). data = seedrates. 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. (1994). The places where the cubics join are known as `knots'.scale" > lines(cars$speed.lm5.frame(fit=pred2[. hat2$upper. where each basis function is a transformation of the variable."fit"]. Confidence bounds for a fitted line spread out more slowly..data=cars) > ci5<-predict(cars.df <.5. 28 Data are from Hand.ci5$fit[.lty=3) # NB assumes values of speed are sorted # try for a closer fit (1 knot) # By default. Lunn. newdata = new.xlab="Seeding rate".frame(rate = c((4:14) * 12.6 Using Factors in R Models Factors are crucial for specifying R models that include categorical or "factor" variables. The fitting of polynomial functions was a simple example of this. give output in which some of the numbers are different and must be interpreted differently.lm(dist~bs(speed. > plot(dist~speed.5)) seedrates. xlim = c(50. eds.df.lm<-lm(dist~bs(speed). but are even less believable! 5.lty=2) lines(rate. interval="confidence") hat2 <.4 Spline Terms in Linear Models By now. in order to show how the confidence bounds spread out.lty=2) 5. The data frame cars is in the datasets package..lm2.ylab="Grains per head") new.interval="confidence". D.hat2$lower.df) The extrapolation has deliberately been taken beyond the range of the data. upper=pred2[. It turns out that. Daly.5.se. ylim = c(15. A. . A Handbook of Small Data Sets. 22). Consider data from an experiment that compared houses with and without cavity insulation28. data = seedrates) pred2 <.data=cars) > library(splines) > cars. E.5)."upr"].lm) > lines(cars$speed.."upr"]) attach(new.5.predict(seedrates. and a set of contrasts.df) lines(rate.fit" "df" "residual. pch = 16. F."lwr"].ci5$fit[.lty=2) detach(new. spline functions of a variable can be constructed as a linear combination of basis functions. 175). it makes a simple example to illustrate the choice of a baseline level. Ostrowski.lm5 <. While one would not usually handle these calculations using an lm model.data. hat2$fit) lines(rate."fit"]) > lines(cars$speed. and depending on the class of spline curves that are used.43 5.ci5$fit[.. although they fit equivalent models."lwr"]. The splines that I demonstrate are constructed by joining together cubic curves. Chapman and Hall.3 Pointwise confidence bounds for the fitted curve Here is code that gives pointwise 95% confidence bounds.lm2 <. J.hat.data=cars) > hat<-predict(cars. D.lower=pred2[. in such a way the joins are smooth. there are no knots > cars.
14683. 14683.. 12669.8e-07 0. Here are the X and the y that are used for the calculations. 12086. 0 0 ... corr=F) Call: lm(formula = kWh ~ insulation) Residuals: Min -4409 1Q Median -979 132 3Q 1575 Max 3690 Coefficients: Estimate Std..lm <.0223 F-statistic: 6. which may be taken to indicate (p < 0. 10689.05) that we can distinguish between the two types of houses. 9708 6700 Residual 10225-10993 10689-10993 .29 p-value: 0... 7))) > # 8 without. 8). 8162.022 Residual standard error: 2310 on 13 degrees of freedom Multiple R-Squared: 0. 6584.44 We begin by entering the data from the command line: insulation <. > insulation <.. 9708-7980 6700-7980 . 9708.lm. and that the initial level is taken as the baseline...6.. 8017. So with comes before without.. 7980+0 7980+0 . 8162. 1 1 ..lm(kWh ~ insulation) > summary(insulation. 12467. The p-value is 0. Error t value Pr(>|t|) (Intercept) insulation 7890 3103 874 1196 9.factor(c(rep("without". Add to get 7980+3103=10993 7980+3103=10993 kWh Compare with 10225 10689 . 8541. and the factor insulation as the explanatory variable. 6584.. Hence: Average for Insulated Houses = 7980 To get the estimate for uninsulated houses take 7980 + 3103 = 10993. Adjusted R-squared: 0. 8541.1 The Model Matrix It often helps to keep in mind the model matrix or X matrix. 10315.59 5. we need to know that the factor levels are. and with is the baseline. 3103 1 1 . we take kWh as the dependent variable. .factor(c(rep("without".. But what does the "intercept" of 7890 mean. rep("with". taken in alphabetical order. 5. 12467. by default. 10689. Note that the first eight data values were all withouts: Contrast 7980 1 1 . 4307... 6700. The standard error of the difference is 1196.. 8022) > insulation. 9708. + 12669.03 2... 4307.. 8).. 12086. and what does the value for "insulation" of 3103 mean? To interpret this. 7))) # 8 without. 8017.73 on 1 and 13 degrees of freedom. then 7 with > kWh <.c(10225. & is the baseline kWh <.022.c(10225. 8022) To formulate this as a regression model.. 6700. then 7 with # `with' precedes `without' in alphanumeric order..341. 10315. rep("with"..
matrix(kWh~insulation) and check that it gives the above model matrix.73 on 1 and 13 degrees of freedom. It is possible to set the choice of contrasts for each factor separately.45 Type model.lm <. i. use either the treatment contrasts or the sum contrasts.29 p-value: 0. "contr. The sum contrasts are sometimes called "analysis of variance" contrasts. so that it becomes the baseline. Adjusted R-squared: 0.e. "mean for with") = 9442 To get the estimate for uninsulated houses (the first level). contr=treatment) Also available are the helmert contrasts.factor(insulation. specify: insulation <.relevel(insulation. These are not at all intuitive and rarely helpful. Here is the output from use of the `sum' contrasts29: > options(contrasts = c("contr. with a statement such as: insulation <.poly".lm(kWh ~ insulation) > summary(insulation.59 7. Here is the interpretation: average of (mean for "without". digits = 2) # Try the `sum' contrasts > insulation <.C(insulation. "with")) # Make `without' the baseline > insulation. corr=F) Call: lm(formula = kWh ~ insulation) Residuals: Min -4409 1Q Median -979 132 3Q 1575 Max 3690 Coefficients: Estimate Std. So the effect for the second level (`with') is -1551. For this.4e-10 0.022 Residual standard error: 2310 on 13 degrees of freedom Multiple R-Squared: 0. Error t value Pr(>|t|) (Intercept) insulation 9442 1551 598 598 15. *5.sum". 29 The second string element. .0223 F-statistic: 6.6. With the "sum" contrasts the baseline is the mean over all factor levels. Novices should avoid them30. even though S-PLUS uses them as the default. "contr.2 Other Choices of Contrasts There are other ways to set up the X matrix.lm. Thus to get the estimate for insulated houses (the first level). take 9442 .78 2. One obvious alternative is to make without the first factor level. there are other choices of contrasts. the user has to calculate it as minus the sum of the remaining effects. is the default setting for factors with ordered levels. The effect for the first level is omitted.1551 = 7980. With the sum contrasts the baseline is the overall mean. levels=c("without".] 30 In general. baseline="without") # Make `without' the baseline Another possibility is to use what are called the "sum" contrasts. [Use the function ordered() to create ordered factors. In technical jargon.341. take 9442 + 1551 = 10993 The `effects' sum to one.poly").
024 6.6 on 1 and 14 degrees of freedom.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.] C: Two separate lines: y = a1 + a2 x1 + b1 x2 + b2 x1 x2 [For the first group (Delphinus delphis. . In the example that follows. x1 = 0) the constant term is a1 and the slope is b1. we had weights for a porpoise species (Stellena styx) and for a dolphin species (Delphinus delphis). . 16 [1] 0 1 2 attr(.lm1.5e-07 # Plot the data Residual standard error: 0. We take x2 to be body weight. data = dolphins) Residuals: Min 1Q Median 0. data=dolphins) > options(digits=4) > cet.827 p-value: 6. data=dolphins) Check the model matrix: > model.] We show results from fitting the first two of these models. A and B: > plot(logheart ~ logweight.matrix(cet.989 1 1 1 1 3.111 on 14 degrees of freedom Multiple R-Squared: 0. and 1 for Stellena styx. data = dolphins) > summary(cet.555 3. and the slope is b1 + b2. while for the second group (Delphinus delphis.08249 Coefficients: Estimate Std. Here we take advantage of this to fit different lines to different subsets of the data.951 attr(. 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.738 .21858 -0.133 0.325 1.treatment" 1 0 3.lm2) (Intercept) factor(species) logweight 1 2 .838. x1 = 0) the constant term is a1. . x1 = 1) the constant term is a1 + a2.lm(logheart ~ factor(species) + logweight. or interactions between variables and factors. For model B (parallel lines) we have > cet.04981 Max 0.lm2 <. x1 = 1) the constant term is a1 + a2.133 2. 8 . .lm(logheart ~ logweight.lm1 <. i.51e-007 F-statistic: 72. We take x1 to be a variable that has the value 0 for Delphinus delphis.00274 3Q 0."assign") 1 0 3.e. or interactions between factors. corr=F) Call: lm(formula = logheart ~ logweight. ."contrasts") [1] "contr.46 5.522 0. Adjusted R-squared: 0. For the second group (Stellena styx.52 0.15874 -0. .54 8. Error t value Pr(>|t|) (Intercept) logweight 1.
for data where all combinations of factor levels have the same number of observations. R uses the treatment contrasts for factors.8461 0.0959 0. .lm(PlantGrowth.Df Res.cet.4921 > summary. .e.1717 0. The parameter ref can be used to set the level that you want as the reference level. aov provides. and plot(cet.lm2) to get an output summary. In essence."contrasts")$"factor(species)" [1] "contr.lm(logheart ~ factor(species) + logweight + factor(species):logweight.8.lm3) (Intercept) factor(species) logweight factor(species).Sum Sq Df Sum Sq F value Pr(>F) 1 2 3 14 13 12 0. 8 .0069 0.aov <.7346 5. See section 5.0758 1 0.3886 Pr(>F) 4. .7663 1.0010 10.01591 > PlantGrowth. data=dolphins) Check what the model matrix looks like: > model.logweight 1 .28 0.000 1 1 3.aov(weight~group) 27 10.lm1. 16 [1] 0 1 2 3 attr(. .989 0.2 below.cet.1 Plant Growth Example Here is a simple randomised block design: > attach(PlantGrowth) # PlantGrowth is from the base datasets package # Looks OK > boxplot(split(weight.555 3.8.8 aov models (Analysis of Variance) The class of models that can be directly fitted as aov models is quite limited.group)) > summary(PlantGrowth. another view of an lm model.951 0.matrix(cet.treatment" 1 0 3.0949 1 0. For model C.aov) Df group Residuals 2 Sum Sq Mean Sq F value 3. .lm2) to plot diagnostic information for this model. A useful function is relevel().lm3) Analysis of Variance Table Model 1: logheart ~ logweight Model 2: logheart ~ factor(species) + logweight Model 3: logheart ~ factor(species) + logweight + factor(species):logweight Res.8832 0.aov) Call: aov(formula = weight ~ group) . > anova(cet."assign") 1 0 3. the statement is: > cet. By default. i. the first level is taken as the baseline or reference level.000 attr(.555 Now see why one should not waste time on model C.47 Enter summary(cet. 5. One can however specify the error term that is to be used in testing for treatment effects.12 0. .lm2.lm3 <.
J. Influence of time of shading on flowering and yield of kiwifruit vines.93 Pr(>F) 4.aov) Df Cult Date Residuals 2 Sum Sq Mean Sq F value 909. and shading from February to May.. in one case from the east and in the other case from the west.aov) Error: block:shade Df block shade Residuals 2 6 Sum Sq Mean Sq F value 172.179e-09 9. W..0710 -0. They are from an experiment31 where there were four treatments .4940 0.48 Residuals: Min 1Q Median 3Q 0.4180 -0. The two papers have different shorthands (e.772 <2e-16 0.65 47.data=cabbages) # cabbages is from the MASS package Examination of the plot suggests that cultivars differ greatly in the variability in head weight.0877 Residual standard error: 0.35 125.84 22. P.3690 -1.2112 0. Journal of Horticultural Science 67: 481-487.527 -1.2788 0. (If this is not specified. were thought more important.aov<-aov(yield~block+shade+Error(block:shade). Further details. are in Maindonald (1992). including a diagram showing the layout of plots and vines and details of shelter.2 Shading of Kiwifruit Vines These data (yields in kilograms) are in the data frame kiwishade that accompanies these notes.331 1. Manson. shading from December to February.2788 25.aov<-aov(VitC~Cult+Date.no shading. Error t value Pr(>|t|) (Intercept) grouptrt1 grouptrt2 5. In experimental design parlance.data=cabbages) > summary(VitC. > help(cabbages) > names(cabbages) [1] "Cult" "Date" "HeadWt" "VitC" > coplot(HeadWt~VitC|Cult+Date.1176 0.074879 3 1394. For the remaining two blocks shelter effects.0411 1.2627 Max 1.01591 F-statistic: 4. Martin.J.57 86. But the F-statistics and p-values will be wrong.0320 -0.846 on 2 and 27 degrees of freedom.1971 0. one still gets a correct analysis of variance breakdown.30 454. 1992.g.15 53.0002486 56 2635.51 464. The northernmost plots were grouped in one block because they were similarly affected by shading from the sun. Each treatment appeared once in each of the three blocks.3710 0. the four vines within a plot constitute subplots.8.06 Pr(>F) 1 2496.40 *5.0060 Coefficients: Estimate Std.) > kiwishade$shade <. Adjusted R-squared: 0.1944 0.6234 on 27 degrees of freedom Multiple R-Squared: 0. Results are given for each of the four vines in each plot. . shading from August to December. Variation in the vitamin C levels seems relatively consistent between cultivars.data=kiwishade) > summary(kiwishade. ref="none") > ## Make sure that the level "none" (no shade) is used as reference > kiwishade.001194 31 Data relate to the paper: Snelgar.2641.15 2496.17 20.relevel(kiwishade$shade. > VitC. The block:shade mean square (sum of squares divided by degrees of freedom) provides the error term.P.2096 p-value: 0. Sept-Nov versus Aug-Dec) for describing the time periods for which the shading was applied. P.6609 0.
58 > coef(kiwishade. 52 distance (cm): 183.e. 50 30. regress time on distance and climb. Using the data frame hills (in package MASS).3."Mater").18 5. 60 distance (cm): 71. Does the quadratic curve give a useful improvement to the fit? If you have studied the dynamics of particles.c(2. 6. Which of these regression equations would you prefer? Use the method of section 5. 249. and sum contrasts. formally. 50.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. 187. Here are two sets of data that were obtained the same apparatus. can you find a theory that would tell you how stopping distance might change with speed?5. Hunter. 150. and plot the line. 189. determine the regression of stretch on distance. 54. whether one needs different regression lines for the two data frames elastic1 and elastic2.aov) (Intercept) : (Intercept) 96.5327 block:shade : blocknorth 0.030833 -10. Do the two sets of results appear consistent. 196. For each of the data sets elastic1 and elastic2. 228.8. this time keeping the factor levels in the order RNC. 45. Using a different symbol and/or a different colour. Then try fitting and plotting a quadratic curve. treatment contrasts.5. In section 5. For the data set elastic2. 127. In each case determine (i) fitted values and standard errors of fitted values and (ii) the R2 statistic. plot distance (i. obtained by supplying the name of the lm object as the first parameter to plot(). 35.2) hosp fhosp<-factor(hosp) levels(fhosp) Now repeat the steps involved in forming the factor fhosp.7 to determine. stopping distance) versus speed. Mater. 48. 114.428333 12. Using an outcome variable y <. 178. Fit a line to this relationship. What is the key difference between the two sets of data? Using the data frame beams (in the data sets accompanying these notes). For the data set elastic1.9 Exercises 1. 42. including the same rubber band.49 Error: Within Df Sum Sq Mean Sq F value Pr(>F) Residuals 36 438.281667 -7.10. Do this using helmert contrasts. Repeat the exercise for the helmert contrasts. Carefully examine the regression diagnostic plot. as the data frame elasticband. the values are: stretch (mm): 46. 50. 217. 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). plot the data from the two data frames elastic1 and elastic2 on the same graph. Type hosp<-rep(c("RNC". Compare the two sets of results. 55. 40. 217. Use contrasts(fhosp) to form and print out the matrix of contrasts. the values are: stretch (mm): 25. What does this indicate? Using the data frame cars (in the datasets package). 291."Hunter".9) . carry out a regression of strength on SpecificGravity and Moisture.430000 3. 208. 44.993125 Within : numeric(0) blockwest shadeAug2Dec shadeDec2Feb shadeFeb2May -3. 249.
x2. To what extent is it potentially misleading? Also do the analysis where the block:shade term is omitted altogether. 1988. quadratic. 1999. E. Geology 11: 195-1976. N. Venables. with accompanying 95% pointwise confidence bounds. do they require transformation. and Mark. C. Introductory Statistics with R. D. Sage Books. A. Statistical design.H.2). Williams. American Economic Review 76: 604-620. As the explanatory variables are proportions. Maindonald J H and Braun W J 2007. Modern Applied Statistics with S. Harrell.8. Multivariable Prognostic Models: Issues in Developing Models. Atkinson. Dehejia. x3 and x4 (which are proportions of four constituents). New York.10 References Atkinson. Multivariate and Tree-Based Methods . A logarithmic transformation is too extreme. Springer. repeating the fit for each of the three different choices of contrasts. Applied Regression including Computing and Graphics. Comment on that analysis. Journal of the American Statistical Association 94: 1053-1062. and Wahba. Statistics in Medicine 15: 361-387.3 to fit: (a) a line. Causal effects in non-experimental studies: re-evaluating the evaluation of training programs. 2nd edn. Applied Linear Regression. Fox. 5. Springer.. New Zealand Journal of Agricultural Research 35: 121-141. Cook. S. and Measuring and Reducing Errors. C.50 fit the model lm(y~fhosp). R. and Weisberg. Evaluating Assumptions and Adequacy.. 1983. Statistical Science 1. and (b) a cubic curve. Tutorial in Biostatistics. F. S. Data Analysis and Graphics Using R – An Example-Based Approach. [Transformation of temperature makes sense only if one first converts to degrees Kelvin. Evaluating the economic evaluations of training programs. P 2002. the direction of the curvature changes. Wiley. 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). Weisberg.. Begin by examining the scatterplot matrix. B. Technometrics 30: 311-318. What family of transformations might one try? Modify the code in section 5. 1986. P. Wiley. Comment on the output that you get from summary(). Dalgaard. but replacing Error(block:shade) with block:shade. Lee. 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). Cambridge University Press. J 2002. 1986. analysis and presentation issues. S. A. 1985. D. and Ripley.5. examine the dependence of pressure on temperature. New York. Maindonald J H 1992. Comment on what you get. 4th edn 2002. Which of the three possibilities (line. curve) is most plausible? Can any of them be trusted? *Repeat the analysis of the kiwishade data (section 5.. Transformations Unmasked. Improper use of regression equations in the earth sciences. Comment: Aspects of diagnostic regression analysis. L. W. Lalonde. D. K. R. B. G. examine the dependence of y (amount of heat produced) on x1. 397–402. 1999. An R and S-PLUS Companion to Applied Regression. R. 2nd edn. Consider transformation of pressure. 1996. with accompanying 95% confidence bounds.
6. 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. An analysis that works with the unscaled variables.e. Multivariate and Tree-based Methods 6. L. C. The second principal component is the component that. B. col=palette()[as. Cunningham. This is because the ratio of largest to smallest value is relatively small. F. never more than 1. K. It is good practice to begin by examining relevant scatterplot matrices. 1995. A plot in which the sites and/or the sexes are identified will draw attention to any very strong structure in the data. Because all variables measure much the same quantity (i. With biological measurement data. groups=possum$site. and so on.1 Multivariate EDA. Australian Journal of Zoology 43: 449-458. and because the standard deviations are typically fairly comparable.integer(possum$site)]) here<-!is. gives a greater weight to variables that have a large variance. This may draw attention to gross errors in the data. for some or all of the variables. The measure of variation used is the sum of the variances of variables. for all variables. For example one site may be quite different from the others.princomp(log(possum[here. and Principal Components Analysis Principal components analysis is often a useful exploratory tool for multivariate data.key=list(columns=3)) # Principal components # Print scores on second pc versus scores on first pc.6:14]. The first principal component is the component (linear combination of the initial variables) that explains the greatest part of the variation.. scaling to give equal variances is unnecessary. B. Here are some of the scatterplot matrix possibilities: pairs(possum[. it is usually desirable to begin by taking logarithms.6:14])) # by populations and sex.2] ~ possum.51 6. .prc <. The common alternative –scaling variables so that they each have variance equal to one – is equivalent to working with the correlation matrix. The data set possum has nine morphometric measurements on each of 102 mountain brushtail possums. Taking logarithms of these data does not make much difference to the appearance of the plots. auto. Trichosurus caninus Ogilby (Phalangeridae: Marsupiala). explains the greatest part of the remaining variation.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. Morphological variation among columns of the mountain brushtail possum. D.6:14]. Viggers. among linear combinations of the variables that are uncorrelated with the first principal component.. col=palette()[as. 32 Data relate to the paper: Lindenmayer. R. and Donnelly. perhaps after scaling so that they each have variance one.prc$scores[.integer(possum$sex)]) pairs(possum[.prc$scores[.. and hence with the variance-covariance matrix. identified by site xyplot(possum. trapped at seven sites from southern Victoria to central Queensland32. relative variability).1]|possum$Pop[here]+possum$sex[here]. The standard deviations then measure the logarithm of relative change.
The function hclust() does hierarchical agglomerative clustering. and algorithms based on iterative relocation. which future patients will remain free of disease symptoms for twelve months or more. A cruder distinction is between populations. we may wish to predict. I have not thought it necessary to take logarithms. our interest is in supervised classification. For example. cluster analysis is a form of unsupervised classification. The mva package has the cluster analysis routines. and want a rule that will allow us to predict the group to which a new data value will belong. based on prognostic measurements and outcome information for previous patients. It is "unsupervised" because the clusters are not known in advance. The function kmeans() (k-means clustering) implements iterative relocation. by a prior use of a hierarchical agglomeration algorithm.lda(site ~ hdlngth+skullw+totlngth+ # Only if not already attached. with a choice of methods available. In iterative relocation. the algorithm starts with an initial classification. A judgement is then needed on the point at which further merging is unwarranted. sites in Victoria (an Australian state) as opposed to sites in other states (New South Wales or Queensland). 33 References are at the end of the chapter. 6. > library(MASS) > here<. using the lda() function from the Venables & Ripley MASS package. on the basis of morphometric measurements. The output yields a hierarchical clustering tree that shows the relationships between observations and between the clusters into which they are successively merged. for the possum data. I discuss this further below.na(possum$footlgth) > possum. In the language of Ripley (1996). that it then tries to improve. Because it has little on the distribution of variable values. i. Here are calculations for the possum data frame. 6. by population and by sex. There are two types of algorithms – algorithms based on hierachical agglomeration.e. . How does one get the initial classification? Typically. to distinguish animals from different sites.3 Discriminant Analysis We start with data that are classified into several groups. The function dist() calculates distances.52 Figure 21: Second principal component versus first principal component. Our interest is in whether it is possible.2 Cluster Analysis In the language of Ripley (1996)33. Groups that are "close" to one another are then successively merged.!is. In hierarchical agglomeration each observation starts as a separate group.lda <.
from MASS package glass.53 + taillgth+footlgth+earconch+eye+chest+belly. One advantage of such methods is that they automatically handle non-linearity and interactions. but working with the logarithms of measurements.5078 1.lda. for the canonical variates in turn. data=fgl) plot(glass.rpart(type ~ RI+Na+Mg+Al+Si+K+Ca+Ba+Fe.4 Decision Tree models (Tree-based models) We include tree-based classification here because it is a multivariate supervised classification. The singular values are the ratio of between to within group sums of squares. also known as "Classification and Regression Trees" (CART). The reader may wish to repeat the above analysis. Output includes a "decision tree" that is immediately useful for prediction. which implies a linear relationship between the logarithms of the measurements. text(glass. library(rpart) # Use fgl: Forensic glass fragment data. Clearly canonical variates after the third will have little if any discriminatory power.7772 Figure 22: Scatterplot matrix of first three canonical variates. or discrimination.data=possum. subset=here) > options(digits=4) > possum.lda() to get (among other information) scores on the first few canonical variates.9372 # Examine the singular values 3. as in Figure22 3.1860 1. may be suitable for regression and classification problems when there are extensive data. One can use predict.1420 0. Note that there may be interpretative advantages in taking logarithms of biological measurement data.tree). 6.7578 > > plot(possum.tree) . The standard against which patterns of measurement are commonly compared is that of allometric growth. Differences between different sites are then indicative of different patterns of allometric growth.lda$svd [1] 15. Tree-based models. dimen=3) > # Scatterplot matrix for scores on 1st 3 canonical variates. method.tree <. Tree-based methods seem more suited to binary regression and classification than to regression with an ordinal or continuous dependent variable. A tree-based regression approach is available for use for regression problems.
R. (1998). New York. Friedman. Maindonald J H and Braun W J 2008. An Introduction to Recursive Partitioning Using the RPART Routines. as in section 6. Using the data set painters (MASS package).. The MASS package has the Aids2 data set. J. 6.6 References Chambers. The data set Cars93 is in the MASS package.ac. B. Therneau. Pacific Grove CA. J. Data Analysis and Graphics Using R – An Example-Based Approach. Modern Applied Statistics with S-Plus. It integrates cross-validation with the algorithm for forming trees. Investigate the comparison between (i) USA and non-USA cars. Use tree-based classification (rpart()) to identify major influences on survival. and include these in the principal components analysis. Investigate discrimination between plagiotropic and orthotropic species in the data set leafshape34. 2nd edn 1997. and Ripley. and Maindonald. 1999. This is one of two documents included in: http://www. D. 3. and about crossvalidation. but this time using the function lda() from the MASS package to derive canonical discriminant scores. J. Ripley. Plot a scatterplot matrix of the first three principal components. D. T. you also need to know about approaches to pruning trees. 2nd edn. 6. and Hastie. and (ii) the six different types (Type) of car. Repeat the calculations of exercises 1 and 2. Hastie. 4. Available from the internet.3. and Expression. W.uk/pub/SWin/rpartdoc. N. Methods for reduction of tree complexity that are based on significance tests at each individual node (i. branching point) typically choose trees that over-predict. Tree architecture in relation to leaf dimensions and tree stature in temperate and tropical rain forests.54 summary(glass. 1996. Using the columns of continuous or ordinal data. Journal of Ecology 87: 1012-1024. Wadsworth and Brooks Cole Advanced Books and Software. Cambridge University Press. Pattern Recognition and Neural Networks. 1997. 1992.ox. B. T. M. Cambridge University Press.stats. and Atkinson. Colour. 5.e. 2. The Atkinson and Therneau rpart (recursive partitioning) package is closer to CART than is the S-PLUS tree library. determine scores on the first and second principal components. Drawing. . E.zip Venables. Springer. and Tibshirani. D. Examine the loadings on the first three principal components.. apply principal components analysis to the scores for Composition. Now create a new data set in which binary factors become columns of 0/1 data. using different colours or symbols to identify the different schools.H. J. Cambridge UK.5 Exercises 1. T. Statistical Models in S. Additive logistic regression: A statistical view of boosting. M.A. 34 Data relate to the paper: King.tree) To use these models effectively. containing de-identified data on the survival status of patients diagnosed with AIDS before July 1 1991. J.
55 *7.2.4)). c(4. Entering 15:5 will generate the sequence in the reverse order. in place of c(4. One can use negative numbers to omit elements. >.4) [1] 2 3 5 2 3 5 2 3 5 2 3 5 > If instead one wants four 2s.6. R Data Structures 7.5).3.4. 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] <. 15. for vectors that have named elements: > c(Andreas=178.2 Patterned Data Use 5:15 to generate the numbers 5.3. Be sure to use is. 7.NA) > is. >=. enter rep(c(2. As x==NA gives a vector of NAs.3.2) two common ways to extract subsets of vectors: Specify the numbers of the elements that are to be extracted.4.4. Users who do not carefully consider implications for expressions that include Nas may be puzzled by the results.1 Vectors Recall that vectors may have mode logical. note that the function rep() has an argument length. 3. 7. ==.5).out." 7. rep(4.2 Missing Values In R. the missing value symbol is NA.3. In addition to the above. you get no information at all about x.1.6.na(x) > x==NA > NA==NA [1] NA # TRUE for when NA appears. <=.4)) we could enter rep(c(2.na(x) to test which values of x are NA. Specifically.x[x>2] generates the error message "NAs are not allowed in subscripted assignments".3.5). Jeff=183)[c("John".5). To repeat the sequence (2. .3.c(1. 6. then four 5s.4)) [1] 2 2 2 2 3 3 3 3 5 5 5 5 # An alternative is rep(c(2. !=. This applies also to the relations <. c(4.3.1 Subsets of Vectors Recall (section 2.4. Thus suppose we want to extract values of x that are greater than 10. 4) thus: > rep(c(2. Any arithmetic operation or relation that involves NA generates an NA.1. then four 3s. The first four compare magnitudes. each=4) Note further that. ….3)).5). note that x==NA generates NA. John=185. numeric or character. So a further possibility is that in place of rep(c(2.out. For example > x <.c(4. enter rep(c(2. == tests for equality. > rep(c(2.4) we could write rep(4. The following demonstrates a third possibility. The elements that are extracted are those for which the logical value is T. Specify a vector of logical values. 5) four times over.5). and != tests for inequality. meaning "keep on repeating the sequence until the length is length."Jeff")] John Jeff 185 183 A vector of names has been used to extract the elements.5).3).
3. and several other packages. . use the library() command.36667 No. in which different columns may have different modes."yield")] variety 66 72 78 84 90 96 102 108 114 yield Manchuria 22.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. A data frame is a generalisation of a matrix. specify data(package="datasets") and similarly for any other package. In place of c("variety".50000 Peatland 31. We will have more to say on missing values in the section on data frames that now follows. 38 29.g. To attach any other installed package. use of e.3 Data frames The concept of a data frame is fundamental to the use of most of the R modelling and graphics functions.70000 No."yield") we could have written. 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. To turn a data frame of numeric data into a matrix of numeric data. of the properties of matrices.36667 120 Wisconsin No.barley[barley$year=="1932" & barley$site=="Duluth".33333 The first column holds the row labels. or all character.46667 Trebi 30.3.23333 Velvet 22.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. more simply. which in this case are the numbers of the rows that have been extracted. to those elements of x that are not NAs. data from an attached package are automatically available. i. The base package. 462 22. in section 7. Lists are discussed below.6. Other packages must be explicitly installed. In most packages.56667 Glabron 25. There are important differences that arise because data frames are implemented as lists. Data frames where all columns hold numeric data have some.4). 475 27. All elements of any column must however have the same mode. 7. + c("variety".60000 No. For remaining sections of these notes. will be used from time to time. on one or both sides as necessary.86667 Svansota 22. 7. 457 22. which comes with the default distribution. To get information on the data sets that are included in the datasets package..matrix(). c(2. the MASS package.e.56 Users are advised to use !is. use as. > Duluth1932 <. are automatically attached at the beginning of the session. at the Duluth site. but not all. all numeric or all factor.na(x) to limit the selection. 7.
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. If you use any missing value symbols other than the default (NA). 36 Specifying as. 7.is = T prevents columns of (intended or unintended) character strings from being converted into factors. The "" will ensure that empty cells are entered as NAs. the first column (country) has the name of the country. The row names are the city names.table("a:/myfile."*".4 million or more.". in millions.57 7. they have a central role in the incorporation of qualitative effects into model and graphics formulae. 37 One way to get mixed text and numeric data across from Excel is to save the worksheet in a .4. one of the columns contains text strings.g. They have an attribute levels that holds the level names.".) 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.strings=c(".csv text file with comma as the separator. Otherwise any appearance of such symbols as *.4. To set tab as the separator. [But watch that none of the cell entries include commas. There may be multiple missing value indicators.") to read the data into R. For example there may be an O (oh) somewhere where there should be a 0 (zero). 7. They provide an economical way to store vectors of character strings in which there are many multiple occurrences of the same strings.strings to recognise other characters as missing value indicators.""). the column is by default stored as a factor with as many different levels as there are unique text strings35.strings=c("NA". it may be necessary to make it into a factor at that time.] 38 Factors are column objects that have mode numeric and class "factor". use read. it is sometimes necessary to use tab ("\t") or comma as the separator. you will probably want to set na. 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.1 Idiosyncrasies The function read.2 below.table() expects missing values to be coded as NA.5 Factors and Ordered Factors We discussed factors in section 2. sep=". as in the above example.table() The function read.2 Missing values when using read.4. and the second column (population) has the urban centre population.csv and is on drive a:. you need to make this explicit see section 7.4. They are stored as integer vectors. Factors have a dual identity. If you have a text file that has been output from SAS.csv". 36 If the column is later required for use as a factor in a model or graphics formula. which are automatically turned into factors.table() is straightforward for reading in rectangular arrays of data that are entirely numeric. If for example file name is myfile. unless you set na. The default separator is white space. Problems may arise when small mistakes in the data cause R to interpret a column of supposedly numeric data as character strings. More crucially. na. Some functions do this conversion automatically.3 Separators when using read."). This copes with any spaces which may appear in text strings. period(. e.is = TRUE. .6. specify sep="\t". or an el (l) where there should be a one (1). 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.4 Data Entry Issues 7.table() confusing may wish to use the parameter setting as. When.3.table() With data from spreadsheets37. 7.
We might prefer countries to appear in order of latitude. Factors have the potential to cause a few surprises. 3. there are inequalities that relate factor levels.4. Ordered factors inherit the attributes of factors. specify e.2) > ordf.6 Ordered Factors Actually.g. If we form a table that has the number of times that each country appears. you get details both of the class of the object. i.3. As in most operations with factors."high").3. To extract the numeric levels 1.6. specify as. as. factors with ordered levels.58 [There are 19 cities in all. 7. Try > stress.5.levels(islandcities$country) > lev[c(7. R sorts the level names in alphabetical order. it is their levels that are ordered. To create an ordered factor. 2. not the integer values. R by default turns it into a factor."medium". To be sure of getting the vector of text strings. specify as. or to turn a factor into an ordered factor.stress>="medium" [1] FALSE TRUE FALSE Later we will meet the notion of inheritance. The levels of an ordered factor are assumed to specify positions on an ordinal scale.character(islandcities$country) To get the integer values.stress) [1] "ordered" "factor" . so be careful! Here are two points to note: When a vector of character strings becomes a column of a data frame. from North to South. R does the translation invisibly.level. use the function ordered(). We can change the order of the level names to reflect this desired order: > lev <. and of any classes from which it inherits.stress<-ordered(stress. ….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 <. Thus: > class(ordf. To be sure of getting the country names.character(country).stress [1] low Levels: [1] medium high low medium high low < medium < high TRUE FALSE FALSE TRUE TRUE > ordf.] Printing the contents of the column with the name country gives the names. levels=lev[c(7.e. Enclose the vector of character strings in the wrapper function I() if it is to remain character.stress<"medium" TRUE FALSE FALSE TRUE > ordf. There are though annoying exceptions that can make the use of factors tricky.level<-rep(c("low".1)]) In ordered factors.numeric(country). and have a further ordering attribute. There are some contexts in which factors become numeric vectors. specify unclass(islandcities$country) By default.5.2.2.factor(islandcities$country."high")) > ordf. When you ask for the class of an object."medium". 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.4. levels=c("low".6.
The matrix construct generalises to array.lm[1].000 4 5 6 13.3] [1. data=elasticband) It is readily verified that elastic. > elastic.893 -27. which may have more than two dimensions.values" "assign" The first list element is: > elastic.] [2. You can get the names of these objects by typing in > names(elastic. or all character.7 Lists Lists make it possible to collect an arbitrary set of R objects together under a single name.lm["coefficients"] or elastic. specify elastic. matrices or more general arrays. Thus consider > xx <.321 3 18.nrow=2) . and often are. You might for example collect together vectors of several different modes and lengths.553571 Alternative ways to extract this first list element are: elastic. scalars. f.553571 The second list element is a vector of length 7 > options(digits=3) > elastic. Note that matrices are stored columnwise. etc. a rag-tag of different objects.1] [.lm(distance~stretch.107 2 -0.e. There is a subtle difference in the result that is printed out.59 7. Thus a matrix is a more restricted structure than a data frame.lm consists of a variety of different kinds of objects. We will use for illustration the list object that R creates as output from an lm calculation. Lists can be. consider the linear model (lm) object elastic. meaning "list element with name coefficients". sections 1.lm$call) [1] "call" *7. Matrices are likely to be important for those users who wish to implement new regression and multivariate methods.2] [.1.lm[1] $coefficients (Intercept) -63.8 Matrices and Arrays All elements of a matrix have the same mode.matrix(1:6.lm$residuals 1 2.lm[["coefficients"]] elastic.lm (c.214 1.lm$coefficients. functions.4 and 2.786 The tenth list element documents the function call: > elastic.lm <.lm) [1] "coefficients" [9] "xlevels" "residuals" "call" "effects" "qr" "terms" "rank" "df. all numeric. The information is preceded by $coefficients.571429 stretch 4.] 1 2 3 4 5 6 # Equivalently.lm[[1]] We can alternatively ask for the sublist whose only element is the vector elastic.321 7 -7.ncol=3) > xx [. For example.571429 stretch 4.residual" "model" [5] "fitted. enter matrix(1:6. data = elasticband) > mode(elastic. i.lm$call lm(formula = distance ~ stretch. For this.4) created thus: elastic. stored as a list. One reason for numeric matrices is that they allow a variety of mathematical operations that are not available for data frames.1.lm$coefficients (Intercept) -63.
4] [1.1] [.4. Consider a numeric vector of length 24. we will make them 1. in order.7] [.12) turns this into a 2 x 12 matrix. So that we can easily keep track of the elements.] [2.4)] x34[2. . This generalises in the obvious way for use with arrays.1] [.] [2.2] [.] x34[-2. The dimnames() function gives a list.8] [. Thus > dim(xx) [1] 2 3 In fact a matrix is a vector (numeric or character) whose dimension attribute has length 2.2] [. . 24.ncol=4) > x34 [.4] [.3] [.4] [1. Thus x <.c(1. .2] [.11] [.] .matrix(1:12.vector(xx) places columns of xx.12] [1. .4] [1..60 If xx is any matrix.8.] [2.] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 Now try > dim(x) <-c(3. > x [. into the vector x.3] [.as.1] [. Now set > x34 <.6] [.3] [.1] [. in which the first list element is the vector of row names. . 2 [. which we now discuss. 7. we get back the elements 1.2] [.9] [. Matrices have the attribute "dimension".] 1 2 3 4 5 6 7 8 9 10 11 12 Here are examples of the extraction of columns or rows or submatrices x34[2:3.10] [.c(2. 1 [. .] [2. 2. the assignment x <.] x34[-2. . A matrix is a 2dimensional array.3] [. 2. and the second list element is the vector of column names.] [3.1:24 Then dim(x) <.] 13 14 16 17 19 20 22 23 1 2 3 4 5 6 7 8 9 10 11 12 ..5] [.] [3.1 Arrays The generalisation from a matrix (2 dimensions) to allow > 2 dimensions gives an array. 6.2) > x .-3] # Extract rows 2 & 3 & columns 1 & 4 # Extract the second row # Extract all rows except the second # Extract the matrix obtained by omitting row 2 & column 3 The dimnames() function assigns and/or extracts matrix row and column names. In the example above.
c(2.2 Conversion of Numeric Data frames into Matrices There are various manipulations that are available for matrices. 102. and cpus (MASS package). on what this reveals. 7. 12. etc. three 3's. …. Create a vector consisting of one 1.matrix() to handle any conversion that may be necessary. in the columns of a 6 by 4 matrix.sum(answer[j]. 3. for each of these data sets. 7. 6.answer[j-1])} answer <. Store the numbers obtained . 3). For each of the following calculations. 4) for (j in 2:length(answer)){ answer[j] <. and store the result in the vector x. and ending with nine 9's. 7. Which of these are ordered factors? Use summary() to get information about data in the data frames attitude (both in the datasets package). then two 2's.. Determine which columns of the data frame Cars93 (MASS package) are factors. but not for data frames. Find the mean of the snow cover (a) for the oddnumbered years and (b) for the even-numbered years. in order. 4) for (j in 2:length(answer)){ answer[j] <. print out the levels vector. Write brief notes. Generate the sequence consisting of eight 4s. 12.c(2. then seven 6s. 1. 3. upper and lower quartiles.9 Exercises Generate the numbers 101. 112. mean.] 15 18 21 24 7.answer[j-1])} In the built-in data frame airquality (datasets package): (a) Determine. (b) Extract the row or rows for which Ozone has its maximum value. Refer to the Eurasian snow data that is given in Exercise 1. the median. Generate four repeats of the sequence of numbers (4. what you would expect? Check to see if you were right! a) b) answer <. extract a data frame which holds only information for small and sporty cars. Use as. .61 [3.6 . For each of these factor columns. and range. 1. From the data frame Cars93 (MASS package). From the data frame mtcars (MASS package) extract a data frame mtcars6 that holds only the information for cars with 6 cylinders. 5. and finally nine 3s. (c) extract the vector of values of Wind for values of Ozone that are above the upper quartile.max(answer[j].8. 5. for each of the columns of the data frame airquality (datasets package).
brandnames[1:5] [1] "Acura" "Acura" "Audi" "Audi" "BMW" . Functions 8. I note two of the simpler functions. " ".1.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. The operator %in% can be highly useful in picking out subsets of data.brandnames <. Numeric vectors will be sorted in numerical order.62 8. For example. This test assumes that counts enter independently into the cells of a table. <vec2>) > order(<vector>) ## For each element of <vec1>. *8. Below.rep(1:5. fixed=TRUE).1 Functions for Confidence Intervals and Tests Use the help to get complete information. <last position>) nchar(<vector of text strings>) ## Returns vector of number of characters in each element. 8.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.test(). This allows both a one-sample and a two-sample test. 8.x %in% c(2. Character vectors will be sorted in alphanumeric order. library(MASS) # if needed To find the position at which the first space appears in the information on make of car.1 The t-test and associated confidence interval Use t. <first position>.rep(3. Alternatively.3 String Functions substring(<vector of text strings>.5)) > x [1] 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 > two4 <. > rank(<vector>) ## Returns the ranks of the successive elements.2 Chi-Square tests for two-way tables Use chisq.test() for a test for no association between rows and columns in the output from table().sapply(strsplit(as. the argument may be a matrix. the test is invalid if there is clustering in the data. 8.character(Cars93$Make).2 Matching and Ordering > match(<vec1>. For example: > x <.1.3. + function(x)x[1]) > car. 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. we might do the following: > car.
8. Output is a list.2 3. > sapply(airquality.3 2.3 3.mean) Ozone Solar. a list of factors. for these data. FUN=mean) Cult Date 1 2 3 4 c39 c52 c39 c52 x d16 3. and a function that operates on a vector to return a single value.80 Day 15.5 aggregate() and tapply() The arguments are in each case a variable. useful information! 8.. : Factor w/ 3 levels "d16".8 4.R 37 7 Wind 0 Temp 0 Month 0 Day 0 Here are several further examples that use the data frame moths that accompanies these notes: > sapply(moths.3 1. The function aggregate() returns a data frame.2 sapply() The function sapply() can be useful for getting information about the columns of a data frame..na(x))) Ozone Solar.93 The use of apply(airquality.R NA NA Wind 9.96 # All elements must be numeric! Temp 77.R 42.1 apply() The function apply() can be used on data frames as well as matrices. <function>) lapply(<list>.is.13 185.11 . For each combination of factor levels.. Here is an example: > apply(airquality. of 4 variables: : Factor w/ 2 levels "c39".99 Month 6.80 > apply(airquality.."c52": 1 1 1 1 1 1 1 1 1 1 .2. but simplify (e.99 Day 15.frame': $ Cult $ Date $ VitC 60 obs. to a vector ## or matrix).96 Wind 9.88 Month 6. : int $ HeadWt: num > attach(cabbages) > aggregate(HeadWt.mean) will give means for each row..18 d16 2.26 d20 2. <function>) ## As lapply().5 2. A dataframe is a list.88 Temp 77..mean. function(x)if(!is. the function is applied to corresponding values of the variable..80 d20 3. 2. by=list(Cult=Cult.. <function>) ## N. <dimension>.4 Application of a Function to the Columns of an Array or Data Frame apply(<array>. B. sapply(<list>.g."d20". function(x)sum(is.na. if possible.2. Here we use it to count that number of missing values in each column of the built-in data frame airquality.factor(x))return(0) else length(levels(x))) meters 0 A 0 P habitat 0 8 *8.1 4. These are not. 51 55 45 42 53 50 50 52 56 49 .7 3.1 .5 4."d21": 1 1 1 1 1 1 1 1 1 1 .4.factor) meters FALSE A FALSE P habitat FALSE TRUE # Determine which columns are factors > # How many levels does each factor have? > sapply(moths.rm=TRUE) Ozone Solar.1.4. Date=Date). For example: > str(cabbages) `data.63 8.
Date) and help(format.summary Min. these would have been omitted from the new data frame. format="%d/%m/%Y") . and convert it into a dates object. except that the name of the second argument is INDEX rather than by."2004/02/22".passengers No.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. and then the day of the month. 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".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. then the month.7 Dates Since version 1.drop=F]. If there had been rows with missing values of Type. The function as. 91 days Use format() to set or change the way that a date is formatted. suitable for use in plotting. By default. The output is an array with as many dimensions as there are factors. This can be avoided by making sure that Type has NA as one of its levels. See help(Dates). dates are stored using January 1 1970 as origin.integer() is used to convert a date into an integer value.as.0.merge(x=Cars93. NA is returned.names") This creates a data frame that has the abbreviations in the additional column with name "abbrev". the date package has been superseded by functions for working with dates that are in R base.Date) for detailed information. while a later column holds two character abbreviations of each of the car types. stored as a factor. I have created a data frame Cars93.summary[. Here are examples: > as. The default is that the year comes first. Here however our demands are simple.y=Cars93.47 The syntax for tapply() is similar. Where there are no data values for a particular combination of factor levels.by."2003/11/23".Date(c("2003/08/24". and we can proceed thus: new.Date() will take a vector of character strings that has an appropriate format. *8.64 5 6 c39 c52 d21 2.Date() to convert text strings into dates. 8.y="row.of. thus: > # Electricity Billing Dates > dd <.9.74 d21 1.Cars93 <. is Type.4."2004/05/23")) > diff(dd) Time differences of 91. 91. in both data frames. The following are a selection of the symbols used: %d: day. One of the variables. > Cars93. in which the row names are the distinct values of Type. by. This becomes apparent when as.summary. Use as. help(as.passengers Max.Date("1/1/1960".
Date("1960-1-1") Time difference of 335 days > as. .Date("2004-12-1") > format(dec1.8.as. In the example above."%d/%m/%Y") [1] "1960-12-31" > as. format="%a %b %d %Y") [1] "Wed Dec 01 2004" 8.000000 15.Date). this was the vector consisting of the two named elements mean and sd.314 1619. so that users can run the function without specifying a parameter. if you wish.function(x=1:10){ + av <.integer(as.Date("1/1/2000".555556 The function returns the value (fahrenheit-32)*5/9. this is enclosed between curly braces ({ }).1 Syntax and Semantics A function is created using an assignment. just to see what it does. Here is a function that prints out the mean and standard deviation of a set of numbers: > mean.integer(as.mean(x) + sd <.60)) [1] 4. SD=sd) + } > > # Now invoke the function > mean.sd(hills$climb) mean SD 1815. format="%b %d %Y") [1] "Dec 01 2004" > format(dec1."%d/%m/%Y") [1] 0 > as. the result is printed."%d/%m/%Y")) [1] 10957 The function format() allows control of the formatting of dates.format="%d:%m:%Y") [1] "1960-12-01" > as. In the example above the default was x = 1:10.444444 10. Following the closing ")" the function body appears. You can.Date("1/1/1970". See help(format.500000 3.Date("1960-12-1")-as. a function returns the value of the last statement of the function. On the right hand side.65 [1] "1960-01-01" > as. > dec1 <. Writing Functions and other Code We have already met several functions.8.sd() mean SD 5. The return value usually appears on the final line of the function body. Unless the result from the function is assigned to a name.sqrt(var(x)) + c(mean=av.and. the parameters appear within round brackets.151 8.Date("31/12/1960".Date("1:12:1960". Except where the function body consists of just one statement.sd <.027650 > mean. Here is a function to convert Fahrenheit to Celsius: > fahrenheit2celsius <.50. give a default.and. More generally.and.function(fahrenheit=32:40)(fahrenheit-32)*5/9 > # Now invoke the function > fahrenheit2celsius(c(40.
3 Compare Working Directory Data Sets with a Reference Set At the beginning of a new session.2]. Settings that may need to change in later use of the function should appear as default settings for parameters. the first argument of sapply() may be a list.8. and if a factor. R allows the use of names for elements of vectors and lists. check. As far as possible. provide a demonstration mode for functions. Remember that they can be abbreviated in actual use. even if this means that they are longer than one would like. So we create a function check. . It will allow us to determine whether x is a factor. Choose their names carefully. [More generally. and may be useful for debugging. data structures and code. where this seems appropriate. thus: dsetnames <.tot[.8. and otherwise the number of levels of x. so that they are meaningful.factor() to test whether x is a factor.66 8. and for rows and columns of arrays and dataframes. function(x)if(!is. make code self-documenting. how many levels it has.factor() will return T if x is a factor. Consider the use of names rather than numbers when you pull out individual elements.8. we might store the names of the objects in the working directory in the vector dsetnames. Choose meaningful names for arguments. defined thus: additions <. we can enter additions(dsetnames) to get the names of objects that have been added since the start of the session. Such a mode will print out summary information on the data and/or on the results of manipulations prior to analysis.df <. Use lists. 8. columns etc. we encountered the function sapply() that can be used to repeat a calculation on all columns of a data frame.4 Issues for the Writing and Use of Functions There can be many functions. Use meaningful names for R objects. with appropriate labelling."dead"] is more meaningful and safer than dead. Where appropriate. The code needed to implement this feature has the side-effect of showing by example what the function does.as. Thus dead. we may want to do similar calculations on a number of different data frames.df() that encapsulates the calculations. function(x)if(!is.2 A Function that gives Data Frame Details First we will define a function that accepts a vector x as its only argument.objects(pos=1) existing <. faclev) We can alternatively give the definition of faclev directly as the second argument of sapply. to group together parameters that belong together conceptually.factor(x))return(0) else length(levels(x))) Finally.factor(x))return(0) else length(levels(x))) 8.factor(x))return(0) else length(levels(x)) Earlier. objnames.] To apply faclev() to all columns of the data frame moths we can specify > sapply(moths. The following function faclev() uses is. thus > sapply(moths. nomatch = 0)) newnames[!existing] } At some later point in the session.function(objnames = dsetnames) { newnames <. It prints out 0 if x is not a factor. and otherwise F.function(x)if(!is.df().logical(match(newnames.function(df=moths) sapply(df. Ensure that the names used reflect the hierarchies of files. > faclev <. The built-in function is. Here is the definition of check.tot[.objects() Now suppose that we have a function additions().
Thus. one might specify attributes(elasticband)$title <"Extent of stretch of band. which is the same as the probability that the student guesses incorrectly.answers thus contains the results of the student's guesses. Structure computations so that it is easy to retrace them. we say that the student guessed correctly.2) correct. We can simulate the correctness of the student for each question by generating an independent uniform random number. R can do this as follows: guesses <.2). but where the parameter if it appears may be any one of several types of data structure. The vector correct.answers <. We can get an idea by using simulation.2 is exactly . 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. Concentrate in one function the task of pulling together data and labelling information. also. This helps ensure that functions contain welltested and well-understood components. If this number is less than .3) for useful functions for routine tasks. Watch the r-help electronic mail list (section 13. which is calculated as TRUE or FALSE. and especially where you may want to retrace your steps some months later.67 Break functions up into a small number of sub-functions or "primitives". The test if(!is. Re-use existing functions wherever possible. labelling etc must be pulled together from a number of sources. A 1 is recorded each time the student correctly guesses the answer.8. to be coerced to 1 (TRUE) or 0 (FALSE).8. while the probability that a uniform random variable exceeds . Bear in mind. 8. from a number of separate files.8. otherwise. with the identification code used to determine what switch() should pick out. You can then generate the file names automatically. to pull out specific information and data that is required for a particular run. the use of switch().2. using paste() to glue the separate portions of the name together.runif(100) correct. give parameters sensible defaults. and Resulting Distance" 8. Each question corresponds to an independent Bernoulli trial with probability of success equal to 0. choose file names that reflect it. while a 0 is recorded each time the student is wrong. Write any new "primitives" so that they can be re-used.1*(guesses < . Use user-defined data frame attributes to document data. For example. Use as the name of the list a unique and meaningful identification code. Often a good strategy is to use as defaults parameters that will serve for a demonstration run of the function. For this reason substantial chunks of code should be incorporated into functions sooner rather than later. This will work. Thus if there is a factorial structure to the data files.6 Graphs Use graphs freely to shed light both on computations and on data. or whether identification by name is preferable.5 Functions as aids to Data Management Where data. Consider whether you should include objects as list items.2 is exactly . take the same care over structuring data as over structuring code. Wherever possible. One of R's big pluses is its tight integration of computation and graphics. . the uniform random number generator is simulating the student. perhaps with some subsequent manipulation.null()) then determines whether one needs to investigate that parameter further.answers The multiplication by 1 causes (guesses<.2.7 A Simulation Example We would like to know how well such a student might do by random guessing. because the probability that a uniform random variable is less than . Structure code to avoid multiple entry of information. NULL is a useful default where the parameter mostly is not required. and makes the function a source of documentation for the data. given the data frame elasticband giving the amount of stretch and resulting distance of movement of a rubber band. we say that the student guessed incorrectly.2. on a multiple choice test consisting of 100 questions each with five alternatives. This structures the code. 8.
Seventeen people rated the sweetness of each of two samples of a milk product on a continuous scale from 1 to 7.function(){ xyrange <. [However the total number of people injured is unlikely to follow a Poisson distribution.7 per year.0833. pch=16) Rewrite this function so that. Repeat for a moving average of order 3.) 11. on [0.8 Poisson Random Numbers One can think of the Poisson distribution as the distribution of the total for occurrences of rare events. a data frame. For example. given the name of a data frame and of any two of its columns. xlim=xyrange. for each patient. 7. the four result against the one result. Write a function that will take as its arguments a list of response variables. 9. generate 50 random integers between 0 and 9. for example the mean or the median. Apply the function to the data (in the data set huron that accompanies these notes) for the levels of Lake Huron. Now modify the function so that you can specify an arbitrary interval. It will return a data frame in which each value for each combination of factor levels is summarised in a single statistic. To simulate the annual number of accidents for a 10-year period. Write a function that generates 100 independent observations on a uniformly distributed random variable on the interval [3. and the variance of such a random variable is 0.) Now. Here is a function that plots. and a function such as mean or median. Write a function that computes a moving average of order 2 of the values in a given vector. 8.75)) abline(0. 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. If no response variable is specified. Why?] The function using rpois() generates Poisson random numbers. We leave this as an exercise. . We pursue the Poisson distribution in an exercise below. Find a way of computing the moving averages in exercise 3 that does not involve the use of a for loop. 2. Suppose for example traffic accidents occur at an intersection with a Poisson distribution that has a mean rate of 3. # Line where four = one plot(four.75 in. with their row and column labels.1]. Create a function to compute the average. 12. an accident at an intersection on any one day should be a rare event. (Compare your results with the theoretical results: the expected value of a uniform random variable on [0. but insisting on the same range for the x and y axes. The supplied data frame milk has columns four and one. variance and standard deviation of such a uniform random variable. The total number of accidents over the course of a year may well follow a distribution that is close to Poisson. Use it to generate 50 random integers between 0 and 99. data=milk. with replacement.7. and January 1 in the year 2000 4. and use it for the same purpose. one. 10.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. 5.75.5. 6. the function will generate random normal data (no difference between groups) and provide the analysis of variance and boxplot information for that.75 in. one sample with four units of additive and the other with one unit of additive. 3. (This means that we do not allow any number to be sampled a second time. Write a function that prints. Determine the number of days. according to R. between the following dates: January 1 in the year 1700. variance and standard deviation of 1000 randomly generated uniform random numbers.9.9 Exercises 1. and January 1 in the year 1800 January 1 in the year 1998.3. showing also the line y=x. Use the round function together with runif() to generate 100 random integers between 0 and 99. 6.1] is 0.8.8].one <. Now look up the help for sample(). plot. it will plot the second named column against the first named column. we can specify rpois(10.7).range(milk) par(pin=c(6.1) } # Calculates the range of all values in the data frame # Set plotting area = 6. Look up the help for the sample() function. 8. sampled without replacement. 5. ylim=xyrange. only those elements of a correlation matrix for which abs(correlation) >= 0. Find the mean. by 6.
Write a function that takes any such vector x as its input. calculates the median of the absolute values of the deviations from the sample median. What are the theoretical values? 16. for any sample.16. 2. 3. and is otherwise zero. Form a vector in which 1's appear whenever the factor level is incremented. Run the function assuming an average rate of 2. 4. Hence it is reasonable to calculate the average w dt dt log w 2 " log w1 over the interval from t1 to t2 as . where there is chance of 1 in 5 of getting the right answer to each question. Compare with the theoretical values of . which is 1 if the light bulb works and 0 if the light bulb does not work. Write a function that. has been calculated from data on heights of trees at times 1.. Find the mean and variance of the student's answers. and the value of the function at the minimum. . for independent random samples from a normal distribution. . 14.5 and . Write an R function that simulates the number of working light bulbs out of 500. . . Write an R function that simulates a student guessing at a True-False test consisting of 40 questions. Write a function that calculates the average of the correlations for any given lag.). df=2) b) xe<-rexp(100) d) xt2<-rt(100. . Write an R function that simulates a student guessing at a multiple choice test consisting of 40 questions. and outputs the vector of factor levels.25. where each bulb has a probability . 17. *19. plotting the result in a suitable way. A "between times" correlation matrix. Generate random samples from normal. Compare with the theoretical values of . f. Compare with the standard deviation in each case. . 4. 6 The first element is the number of occurrences of level 1. 1.8 accidents per year. 4. 1 dw d log w = . Write a function that calculates the minimum of a quadratic. Write a function that does an arbitrary number n of repeated simulations of the number of accidents in a year. The vector x consists of the frequencies 5. thus: a) xn<-rnorm(100) c) xt2<-rt(100. Using simulation.] *21. *20.). and t (1 d. df=1) Apply the function from exercise 17 to each sample. estimate the expected value and variance of the random variable X. Write a function that simulates the repeated calculation of the coefficient of variation (= the ratio of the mean to the standard deviation). 3. exponential. . . *23. .69 13. here 1 1 1 1 1 2 2 2 3 . .2 and . f. t (2 d.] t 2 " t1 [The relative growth rate may be defined as ! ! . 3. 15. the second is the number of occurrences of level 2. Given data on trees at times 1. [You'll need the information that is provided by cumsum(x). 2. and so on. . Find the mean and variance of the student's answers. 18. write a function that calculates the matrix of "average" relative growth rates over the several intervals. Assume that the number of accidents in a year follows a Poisson distribution.99 of working. *22.
Models which have this form may be nested within other models which have this basic form. Thus there may be `predictions' and `errors' at different levels within the total model. In 1 dimension y = "1 (x1 ) + # ! ! z1 = "1 (x1 ). 1# " We can turn this model around. ! Use glm() to fit generalized linear models.z p = " p (x p ) may be smoothing functions.. and General Non-linear Models GLM models are Generalized Linear Models. 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. For logit models. generalizations of this model. and " log( ) = logit(" ) is log(odds). where " log( ) = a + b1 x1 1# " Here ! !is an expected proportion.. In this chapter we describe the alternative functional forms. We can add more explanatory variables: a + b1x1 + .. and we leave until later sections the discussion of different possibilities for the "error" distribution. Generalized Linear Models... allow a variety of non-normal distributions.+ pxp + Use lm() to fit multiple regression models. Generalized Linear Model (e.. GLM.z2 = " 2 (x 2 ). and write ! y = g(a + b1 x1 ) + " = exp(a + b1 x1 ) +" 1+ exp(a + b1 x1 ) Here g(. our focus is on the form of the model prediction. 9. ! Some of Generalized Additive Model ! 39 There are various generalizations. y = " + # .70 *9. in essence.g. .+ " p (x p ) + # Additive models are a generalization of lm models. The various other models we describe are. logit model) y = g(a + b1x1) + Here g(.. and the other extensions we describe. They extend the multiple regression model. + bpxp.) undoes the logit transformation. Multiple regression model y= + 1x1 + 2x2 + . while others may be the usual linear model terms..1 A Taxonomy of Extensions to the Linear Model R allows a variety of extensions to the multiple linear regression model. . The GAM (Generalized Additive Model) model is a further extension. In the discussion of this section.) is selected from one of a small number of options. Additive Model y = "1 (x1 ) + " 2 (x 2 ) + . . The constant term gets absorbed into one or more of the " s..
It is however in order to use a linear model to predict logit(proportion).1 to 0.. We can transform to get the model y = " p (x p ) may be smoothing functions. such as the inverse of the logit function. It is not satisfactory to use a linear model to predict proportions. Some such transformation (`link' function) as the logit is required.. we may still want to retain both "1(. The logit or log(odds) function turns expected proportions into values that may range from . For example. Shown here is a plot of log(odds) versus proportion. it is not reasonable to expect that the expected proportion will be a linear function of x.e.. but log( ! ! Figure 23: The logit or log(odds) transformation. With proportions that range from less than 0.99. . 1" p The linear model predicts.. 9. The function g(. The logit function is an example of a link function. as in a logistic regression model... not p.). i.. then the corresponding odds are p/(1-p). There are various other link functions that we can use with proportions. The values from the linear model may well lie outside the range from 0 to 1. A good way to think about logit models is that they work on a log(odds) scale.. y = g("1 (x1 )) + # . and log(odds) = log( p ) = log(p) -log(1-p) 1" p p ).to + .71 y = g("1 (x1 ) + " 2 (x 2 ) + . One of the commonest is the complementary log-log function.g. 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.) and g(.) doing anything more that seems necessary. that horse A will win the race). If p is a probability (e.z p usual linear model terms. ! The reason is that g(.) may be the function that undoes the logit transformation.+ " p (x p )) + # Generalized Additive Models are a generalisation of Generalized Linear Models.+ z p ) + " ! Notice that even if p = 1. with "1 (.) does as much as it can of the task of transformation. Notice how the range is stretched out at both ends. g(.z 2 = " 2 (x 2 ). while others may be the = g(z1 + z2 + ..) is a specific function. ! Some of z1 = "1 (x1 ).2 Logistic Regression ! We will use a logistic regression model as a starting point for discussing Generalized Linear Models. Figure 23 shows the logit transformation.
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. There is a small number of concentrations. for each of six different alveolar concentrations. A deviance has a role very similar to a sum of squares in regression. or to the proportions in Table 1. .72 9. To understand the output. Logistic regression degrees of freedom deviance mean deviance (divide by d.2. for reasons that will emerge later. Australian National University) for allowing me use of these data.) We prefer models with a small mean residual sum of squares. We can fit the models either directly to the 0/1 data. of proportion of patients not moving. Figure 24 then displays a plot of these proportions. It was then noted whether the patient moved. The horizontal line estimates the expected proportion of nomoves if the concentration had no effect.5 0 2 2 _____________________________________________ Table 1: Patients moving (0) and not moving (1). Figure 24: Plot. jerked or twisted. i. Thus we have: Regression degrees of freedom sum of squares mean sum of squares (divide by d.f. University of Chicago) and to Alan Welsh (Centre for Mathematics & its Applications. the logit model and the complementary log-log model.4 2 4 6 1. 40 I am grateful to John Erickson (Anesthesia and Critical Care.6 0 4 4 2. Alveolar Concentration Nomove 0 1 Total 0.f.2 2 4 6 1.8 6 1 7 1 4 1 5 1. so we begin by tabulating proportion that have the nomove outcome against concentration.e. The response is best taken as nomove. The interest is in estimating how the probability of jerking or twisting varies with increasing concentration of the anesthetic agent.) We prefer models with a small mean deviance. you need to know about "deviances". versus concentration. We fit two models.
8 on 28 degrees of freedom Number of Fisher Scoring Iterations: 5 Correlation of Coefficients: (Intercept) conc -0. The line is the estimate of the proportion of moves.07 Coefficients: Value Std.744 0.logit). they are drawn at random from some larger population. data = anesthetic) Deviance Residuals: Min 1Q Median 3Q Max -1.logit) Call: glm(formula = nomove ~ conc.glm(nomove ~ conc. of log(odds) [= logit(proportion)] of patients not moving.04 -2. 25 is a graphical summary of the results: Figure 25: Plot.47 conc 5.42 2. With such a small sample size it is impossible to do much that is useful to check the adequacy of the model.73 If individuals respond independently. family = binomial(link = logit). family = binomial(link = logit). each time a new individual is taken.68 2. . + data = anesthetic) The output summary is: > summary(anaes. Justification for assuming the same probability will arise from the way in which individuals are sampled. with the same probability. based on the fitted logit model.981 Fig.0341 0.57 2. Error t value (Intercept) -6. Try also plot(anaes. Here is the R code: > anaes.72 (Dispersion Parameter for Binomial family taken to be 1 ) Null Deviance: 41.5 on 29 degrees of freedom Residual Deviance: 27.logit <. versus concentration. then we have Bernoulli trials. While individuals will certainly be different in their response the notion is that.77 -0.687 2.
Also it provides a more helpful form of output.4. maxtemp and dewpoint. for each combination of mintemp and maxtemp. Engineering or business failures can be modelled using this same methodology. data=anesthetic) The family parameter specifies the distribution for the dependent variable. The log link transforms from positive numbers to numbers in the range .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. Australian National University. You need to install the splines package. Geography Department. In place of x one specifies bs(x)or ns(x). of monthly data from a number of different times and locations. 9.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. family = binomial(link = logit).3. but often the default works quite well. This way of formulating an lm type model does however have the advantage that one is not restricted to the identity link.lm <.glm(nomove ~ conc. i. There is an optional argument that allows us to specify the link function.74 9. The default link is then the identity. from datasets package air. R does not at present have a facility for plots that show the contribution of each term to the model. 9.logit <.3 The gaussian family If no family is specified. .to + that a linear model may predict. The R survival package has state of the art abilities for survival analysis. We fit the model: dewpoint = mean of dewpoint + smooth(mintemp) + smooth(maxtemp) Taking out the mean is a computational convenience. Here are details of the calculations: dewpoint.4 Models that Include Smooth Spline Terms These make it possible to fit spline and other smooth transformations of explanatory variables.3 glm models (Generalized Linear Regression Modelling) In the above we had anaes. One can control the smoothness of the curve.3.lm(dewpoint ~ bs(mintemp) + bs(maxtemp). normal errors model summary(air.glm<-glm(Ozone^(1/3) ~ Solar.glm) 9. for making these data available. # Dataset airquality.R + Wind + Temp. data = airquality) # Assumes gaussian family. then the family is taken to be gaussian. 9.e. The dewpoint values are averages. The link that is commonly used here is log. Visiting Fellow. One can request a `smooth' b-spline or n-spline transformation of a column of the X matrix.lm) 9.1 Dewpoint Data The data set dewpoint41 has columns mintemp. Below we give further examples. data = dewpoint) options(digits=3) summary(dewpoint. as for an lm model. 41 I am grateful to Dr Edward Linacre.
The output may not be what you might expect.7 Model Summaries Type in methods(summary) to get a list of the summary methods that are available. They have now become a stock-in-trade of statistical analysts. Modern Applied Statistics with S-Plus. summary.. N. and Tibshirani. An important elaboration is to the incorporation of more than one term in the error structure. 2nd edn 1997.g.. 9. McCullagh. A. Generalized Additive Models. London.. This is fortunate. Maindonald J H and Braun W J 2003. T. Allow different intercepts for different habitats. by the 1970s. R. 9. Most professional users of R will regularly encounter data where the methodology that the data ideally demands is not yet available. Chapman and Hall. An Introduction to Statistical Modelling. W. So be careful! 9. J. They unified a huge range of diverse methodology. Chapman and Hall. e. had been developed for handling linear model calculations.8 Further Elaborations Generalised Linear Models were developed in the 1970s. J. J.lm() on an aov or glm object. D. and Ripley. Data Analysis and Graphics Using R – An Example-Based Approach.75 9. The R nlme package implements such extensions.6 Non-linear Models You can use nls() (non-linear least squares) to obtain a least squares fit to a non-linear function. Use log(meters) as a covariate. Each such new development builds on the theoretical and computational tools that have arisen from earlier developments. Cambridge University Press. Their practical implementation built on the powerful computational abilities that. Practical data analysis demands further elaborations. 2nd edn. Venables. 1983. 1989. P. both for linear models and for a wide class of nonlinear models. J. London. You may want to mix and match. Hastie. 9. Springer. B. Generalized Linear Models.10 References Dobson. and Nelder. A. Exciting new analysis tools will continue to appear for a long time yet.9 Exercises 1. 1990. Chapman and Hall. Fit a Poisson regression model to the data in the data frame moths that Accompanies these notes. New York. .
Multi-level Models.20250 1. which may however be observed at a large number of time points.867629 6 -3.8. 10.76 *10.1558 shadeDec2Feb -10. If we treat laboratories and observers as random. In the time series context there is usually just one realisation of the series.9663 278.53 shadeDec2Feb -0. sep=".53 0. and between multiple determinations made by the same observer on different samples. Here is the analysis: > library(nlme) Loading required package: nls > kiwishade$plot<-factor(paste(kiwishade$block.1 Multi-Level Models.03083 1.761621 36 56.019373 Formula: ~1 | plot %in% block (Intercept) Residual StdDev: 1.53 0. Including Repeated Measures Models Models have both a fixed effects structure and an error structure.50 Standardized Within-Group Residuals: .lme<-lme(yield~shade. between observers within laboratories.42833 1. from the Pinheiro and Bates nlme package. Repeated Measures and Time Series 10. kiwishade$shade.490378 Fixed effects: yield ~ shade Value Std.867629 6 1.478639 3.0073 Correlation: (Intr) shdA2D shdD2F shadeAug2Dec -0. Version 3 of lme is broadly comparable to Proc Mixed in the widely used SAS statistical package. The random effects are block (though making block a random effect is optional. data=kiwishade) > summary(kiwishade.0015 shadeFeb2May -7.867629 6 -5. Again Refer back to section 5. the only fixed effect is the mean.")) > kiwishade. For example. handle models in which a repeated measures error structure is superimposed on a linear (lme) or non-linear (nlme) model.4556 -125.88086 <. in an inter-laboratory comparison there may be variation between laboratories.0001 shadeAug2Dec 3.62282 0.1 The Kiwifruit Shading Data. The function lme has associated with it highly useful abilities for diagnostic checking and for various insightful plots.2 for details of these data.lme) Linear mixed-effects model fit by REML Data: kiwishade AIC BIC logLik 265. In the repeated measures context there may be a large number of realisations of a series that is typically quite short. The fixed effects are block and treatment (shade). There is a strong link between a wide class of repeated measures models and time series models.Error DF t-value p-value (Intercept) 100. plot within block.50 0. The functions lme() and nlme().50520 0.50 shadeFeb2May -0.28167 1.9831 Random effects: Formula: ~1 | block (Intercept) StdDev: 2.1. and units within each block/plot combination. for purposes of comparing treatments).random=~1|block/plot.97741 0.
030833 -7. est.478639 5. By squaring the standard deviations and converting them to variances we get the information in the following table: Variance component block plot residual (within group) 2.lme) Approximate 95% confidence intervals Fixed effects: lower (Intercept) shadeAug2Dec -1.14 12.019373 7.211 36 5190.85159 -10.076 d.490378 4.4902=12.180 Mean square for anova table 12.186 2.0689948 Number of Observations: 48 Number of Groups: block plot %in% block 3 12 > anova(kiwishade.905037 Within-group standard error: lower est. upper sd((Intercept)) 0.7804597 Max 1.428333 upper 7. 2 (3-1) 6 (3-1) (2-1) 3 4 (4-1) . We can get an output summary that is helpful for showing how the error mean squares match up with standard deviation information given above thus: > intervals(kiwishade.711743 -2.858410 96.202500 103. upper 2. We have: Variance component block plot residual (within gp) 4.479 = 2.180 2 2 Notes Three blocks 4 plots per block 4 vines (subplots) per plot The above gives the information for an analysis of variance table.552 Q3 0.62977 100.186 12.f.3702555 1.45086 upper sd((Intercept)) 0.019 = 4.92 12.4153887 -0.18 2.5890938 This was a balanced design.lme) numDF denDF (Intercept) shade 1 3 6 F-value p-value <.0001 0.186 3.397024 We are interested in the three sd estimates.5473014 2. 3.770678 3.186 + 16 4.2 could use aov().180 + 4 = 20.076 1.99826 Random Effects: Level: block lower Level: plot lower est.775232 shadeDec2Feb -14.5981415 -0.77 Min Q1 Med -2.0012 22.180 + 4 = 86. which is why section 5.53909 est.076 2.600757 -5.8.281667 shadeFeb2May -11.
93 Pr(>F) 4. N.lme<-lme(log(it)~(tint+target+agegp+sex)^2.1.074879 3 1394.17 20.data=kiwishade) > summary(kiwishade.926906 1. The variable sex is coded 1 for males and 2 for females. Data are in the data frame tinting.lme) Model df itstar. to get variances that are approximately homogeneous. 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). J.lme. Effects of car window tinting on visual performance: a comparison of elderly and young drivers. In this data frame.lme.lme<-lme(log(it)~(tint+target+agegp+sex). Here is what we get: > anova(itstar.lme it1. the time in milliseconds required to recognise an alphanumeric target).001194 Error: Within Df Sum Sq Mean Sq F value Pr(>F) Residuals 36 438.45036 21.it1.. csoa (critical stimulus onset asynchrony. i. The authors are mainly interested in effects on side window vision.2 The Tinting of Car Windows In section 4. We start with a model that is likely to be more complex than we need (it has all possible interactions): itstar.57 86.51 464.method="ML") A reasonable guess is that first order interactions may be all we need.e.aov<-aov(yield~block+shade+Error(block:shade). the time required for a simple discrimination task) and age are variables. we need to work with log(csoa) and log(it).lme<-lme(log(it)~tint*target*agegp*sex.430915 2 vs 3 22. Ergonomics 42: 428-443.it2. random=~1|id.2112 0.77022 7.e.84 22. data=tinting.1 make it clear that. White. random=~1|id.aov) Error: block:shade Df block shade Residuals 2 6 Sum Sq Mean Sq F value 172.8. and Willson.58 12. T.18 10.lme 1 26 3 8 AIC BIC logLik Test L. i. This is so that we can do the equivalent of an analysis of variance comparison. there is the very simple model. We have two levels of variation – within individuals (who were each tested on each combination of tint and target). it (inspection time.146187 91.1 we encountered data from an experiment that aimed to model the effects of the tinting of car windows on visual performance42.78 Now fsee where these same pieces of information appeared in the analysis of variance table of section 5.1176 0. allowing only for main effects: it1.35 125.138171 26. Nettlebeck.Ratio p-value 6. M. Plots such as we examined in section 4.742883 50.2: > kiwishade. . data=tinting. i.871441 1 vs 2 42 Data relate to the paper: Burns. and between individuals.88105 2 17 -3.e. 1999.method="ML") Note that we have fitted all these models by maximum likelihood. while tint (3 levels) and target (2 levels) are ordered factors.7288 0.. Here we examine the analysis for log(it).11093 0.method="ML") Finally. and hence in visual recognition tasks that would be performed when looking through side windows.0065 8. random=~1|id. R. it2.lme it2.72523 18. data=tinting. So we need to specify id (identifying the individual) as a random effect.
agegp tint.lme1) Linear mixed-effects model fit by REML Data: michelson AIC BIC logLik -546 1113 1142 Random effects: Formula: ~Run | Expt Structure: General positive-definite StdDev Corr (Intercept) 46. with the slope varying randomly between the five experiments of 20 runs each.17e-13 2.3 The Michelson Speed of Light Data Here is an example.7589 145 0.196 3.L. using the Michelson speed of light data from the Venables and Ripley MASS package.0520 145 0.Q.74542 DF t-value p-value 0.88e-02 -1.lme. but now with method="REML".79 The model that limits attention to first order interactions is adequate.lme.63e-01 -1.13075 -0.0584 145 0.method="REML") We now examine the estimated effects: > options(digits=3) > summary(it2. 10.sex -0.0890 145 0.378 1.017 3.sex targethicon.78e-02 -0.90e-02 -1. correlation = corAR1(form = ~ 1 | Expt)) ~ 1 | Expt).104 9.0520 145 0.33164 0. a first order autoregressive model would probably be adequate. setting method="ML" in each case.376 1.19e-02 1.17126 -0.81e-01 -2.lme1 <.agegp tint.as.341 1.targethicon -0.reml)$tTable Value Std.00722 tint. one needs to fit models with and without these effects. > library(MASS) # if needed > michelson$Run <.82e-01 1.83e-02 -0. (I leave this as an exercise!) For purposes of doing this test.443 1. data = michelson.numeric(michelson$Run) # Ensure Run is a variable > mich.13449 -0. A model that ignores the sequential dependence entirely does give misleading results.0492 145 0.L tint.5167 22 22 7.74e-03 -1. A model that has neither fixed nor random Run effects seems all that is justified statistically. Maximum likelihood tests suggest that one needs at least this complexity in the variance and dependence structure to represent the data accurately.0482 145 0.85e-01 1.991 4.3261 22 Because tint is an ordered factor.agegp targethicon.84e-02 -2. polynomial contrasts are used.L. it2.0584 145 0.24012 -1.150 8.49 (Intr) .Q targethicon agegp sex 6.658 8.17e-01 2.5167 0.975 4. We will need to examine the first order interactions individually.1.lme(fixed = Speed ~ Run. We allow the variance to change from one experiment to another. We assume an autoregressive dependence structure of order 1.06972 0.87e-02 -2.13887 -0.Q.0461 145 0.Q. To test this.22658 0. weights = varIdent(form = > summary(mich.20e-01 tint.332 1. and compare the likelihoods.05231 0.07785 0.reml<-update(it2.sex tint. The model allows the determination to vary linearly with Run (from 1 to 20).0933 145 0.547 1.996 4.L. random = ~ Run| Expt.0492 145 0.1010 145 0.sex agegp.09193 tint.Error (Intercept) tint. For this we re-fit the model used for it2.00542 0.targethicon -0.836 6.09794 -0.
Run Residual 3.62 -1 121.29
Correlation Structure: AR(1) Formula: ~1 | Expt Parameter estimate(s): Phi 0.527 Variance function: Structure: Different standard deviations per stratum Formula: ~1 | Expt Parameter estimates: 1 2 3 4 5 1.000 0.340 0.646 0.543 0.501 Fixed effects: Speed ~ Run Value Std.Error DF t-value p-value (Intercept) Run Correlation: (Intr) Run -0.934 Standardized Within-Group Residuals: Min Q1 Med 0.109 Q3 0.740 Max 1.810 -2.912 -0.606 868 -2 30.51 94 2.42 94 28.46 -0.88 <.0001 0.381
Number of Observations: 100 Number of Groups: 5
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".
10.4 References
Chambers, J. M. and Hastie, T. J. 1992. Statistical Models in S. Wadsworth and Brooks Cole Advanced Books and Software, Pacific Grove CA. Diggle, Liang & Zeger 1996. Analysis of Longitudinal Data. Clarendon Press, Oxford.
Everitt, B. S. and Dunn, G. 1992. Applied Multivariate Data Analysis. Arnold, London.
Hand, D. J. & Crowder, M. J. 1996. Practical longitudinal data analysis. Chapman and Hall, London.
Little, R. C., Milliken, G. A., Stroup, W. W. and Wolfinger, R. D. (1996). SAS Systems for Mixed Models. SAS Institute Inc., Cary, New Carolina.
Maindonald J H and Braun W J 2003. Data Analysis and Graphics Using R – An Example-Based Approach. Cambridge University Press.
Pinheiro, J. C. and Bates, D. M. 2000. Mixed effects models in S and S-PLUS. Springer, New York. Venables, W. N. and Ripley, B. D., 2nd edn 1997. Modern Applied Statistics with S-Plus. Springer, New York.
*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"){
Let's see how this works in practice > x <. just the arguments.exp <. Thus > parse(text="mean(x+y)") expression(mean(x + y)) We store the expression in my. function (x) paste (deparse(x).exp2. the syntax is checked and it is turned into code that the R computing engine can more immediately evaluate.expression(mean(x+y)) > my. Setting my. collapse = "\n")) } else stop ("argument is not a function call") } ".exp .expression(mean(x+y)) stores this unevaluated expression in my.y).83 # the first element of a 'call' is the function called # so we don't deparse that.expression("mean(x+y)") > eval(my. R gives you as much information as it thinks helpful.3 Parsing and Evaluation of Expressions When R encounters an expression such as mean(x+y) or cbind(x.exp) [1] 131 > eval(my.parse(text="mean(x+y)") > eval(my. and then evaluate it > my. there are two steps: The text string is parsed and turned into an expression.args(list(x+y. as is obvious when the expression is evaluated."()". unless one explicitly changes it into an expression or part of an expression. print(paste("The function is: lapply (s[-1].exp <. but indicate that the parameter is text rather than a file name. Note that expression(mean(x+y)) is different from expression("mean(x+y)"). Upon typing in expression(mean(x+y)) the output is the unevaluated expression expression(mean(x+y)). and want to turn it into an expression? The answer is to use the function parse(). The expression is evaluated.exp are a little different from what is printed out.e. s[1]. The actual contents of my. foo(bar))) [1] "The function is: [[1]] [1] "x + y" [[2]] [1] "foo(bar)" list ()" 11.txt <.101:110 > y <. collapse="")) For example: > deparse.txt) [1] "mean(x+y)" What if we already have "mean(x+y)" stored in a text string. A text string is a text string is a text string.exp2) [1] 131 .exp2 <.21:30 > my. i.
illustrates some of the possibilities.exit(detach(old. ")". and use only the column names make. .names(list(.leftright[2] xname <. argtxt.call is used it is only necessary to use parse() in generating the argument list. given without explanation. 1:10) else { nam <.frame(".parse(text=right) . colnames = c("NSW".)) if(nam!=xname)stop("Clash of variable names") # The part to the left of the "=" # The part to the right of the "=" expr <.new.". Once in the function.df() NSW ACT 1 1904 . sep = "") expr <. ")".eval(expr) names(df) <. sep = "")) df <.call() makes it possible to supply the function name and the argument list in separate text strings.colnames df } 11.all.paste(colnames..strsplit(equation. i.84 Here is a function that creates a new data frame from an arbitrary set of columns of an existing data frame.call("data.df)) argtxt <.parse(text = exprtxt) df <.function(old. collapse=" & ").paste(xname. "ACT")) { attach(old.df) on. argtxt..4 Plotting a mathematical expression The following. ..do.paste("data. colnames = c("NSW".colnames df } To verify that the function does what it should. eval(listexpr)) names(df) <. It needs better error checking than it has at present: plotcurve <function(equation = "y = sqrt(1/(1+x^2))".. "ACT")) { attach(old. "on the right of the equation")) if(length(list(.leftright[1] right <.") listexpr <.df) on.df <function(old.new. When do.e. 3 The function do.") exprtxt <. For example make.df = austpop.frame". type in > make. collapse = ".df)) argtxt <. .df <. split = "=")[[1]] left <. collapse = ".vars(expr) if(length(xname) > 1)stop(paste("There are multiple variables.))==0)assign(xname.){ leftright <. .new.parse(text=paste("list(".df = austpop.exit(detach(old.. we attach the data frame so that we can leave off the name of the data frame..paste(colnames.
mygrep <function(str) { ## ## ## tempobj <. For example. y.y)) mainexpr <. list(. ylab = yexpr.get(i) if(is.frame(-1)) objstring <. for example function names are tokens. xlab = xexpr. deparse(myfunc)))) objstring <.ls(envir=sys. i) } return(objstring) } Assign the names of all objects in current R working directory to the string vector tempobj mygrep("for") # Find all functions that include a for loop . right)) title(main = mainexpr) } Try plotcurve() plotcurve("ang=asin(sqrt(p))".5 Searching R functions for a specified token.)[[1]]) assign(xname. p=(1:49)/50) 11..85 assign("x".function(myfunc)) if(length(grep(str.character(0) for(i in tempobj) { myfunc <.c(objstring. "==". we search all functions in the working directory.parse(text=paste(left. A token is a syntactic entity.eval(expr) yexpr <..parse(text=xname)[[1]] plot(x. The purpose of using unlist() in the code below is to change myfunc from a list into a vector of character strings. x) } y <.parse(text=left)[[1]] xexpr <. type="n") lines(spline(x.
au/pub/CRAN/windows/windows-9x/ Look in the directory contrib for packages. but often terse. [This has become a text book for the use of S-PLUS and R for applied statistical analysis. Venables. Appendix 1 12. Springer.org/doc/bib/R-books. New York.edu. go to: http://cran.2 Contributed Documents and Published Literature Dalgaard. Data Analysis and Graphics Using R – An Example-Based Approach. D.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. 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 . [This is an intermediate level text. Together with the 'Complements' it gives brief introductions to extensive packages of functions that have been written or adapted by Ripley. Modern Applied Statistics with S. For other books and literature. 4th edn 2002. 2nd edn.aarnet. see http://www. [This is an R-based introductory text.html 12. Springer. B. So it pays to check the CRAN site from time to time. P.au/pub/CRAN/ For Windows 95 etc binaries. look in http://mirror. huron and islandcities. 2002. Also. N. New packages are being added all the time.] Venables. florida. are in the DAAG package. watch for announcements on the electronic mailing lists r-help and r-announce. Cambridge University Press. except beams..html http://www.r-project.r-project. Explanation is careful.1 R Packages for Windows To get information on R packages. Introductory Statistics with R. 12. It assumes a fair level of statistical sophistication.aarnet.3 Data Sets Referred to in these Notes Data sets included in the image file usingR.edu. The units of time for hills are however hours (not minutes) in DAAG.org The Australian link (accessible only to users in Australia) is: http://mirror.org/doc/bib/R-other.] Maindonald J H and Braun W J 2008.86 12. W. NY. and a number of other statisticians. with a biostatistical emphasis.RData Cars93.r-project. and Ripley.
lab=brainaxis) box() identify(log(Animals$body). lab=bodyaxis) axis(2.9)) or rep(c(4. log(Animals$brain). nrow=6.ylab="Brain weight (g)") plot(log(Animals$body).rep(3.names(Animals)) Section 7. (iv) plot(snow.rep(3.rep(6.cover ~ year.4) 3.cover)) Section 2.6 1.1.101:112 2. levels) sapply(tinting[. sapply(tinting. (iii). at=log(brainaxis).112) or x <. (c) 600.4) bodyaxis <-10^seq(-2. xlab="Body weight (kg)". pch=1. at=log(bodyaxis). 2. prod(c(10. volume=volume) for (i in 1:50) {bigprod <.7 1.4) axis(1. 4:6].bigsum+i }. axes=F) brainaxis <.bigprod*i }. is.data=elasticband) 2.2)) plot(Animals$body.3).factor) sapply(tinting[. ncol=4) .matrix(c(rep(4. rep(c(4. x <. 4(ii) prod(1:50) 5.log(Animals$brain).frame(radius=radius. xlab="Body weight (kg)". (ii). is. Animals$brain.7).cover) hist(log(snow$snow. (b) 22.ordered) Section 3.8). bigprod for (i in 1:100) {bigsum <. 3(ii) sum(1:100) 4(i) bigprod <.4 Answers to Selected Exercises Section 1. c(rep(4. labels=row.data. radius <. ylab="Brain weight (g)".3:5)) 3(i) bigsum <.seq(101.c(8. bigsum 6.4*pi*radius^3/3 sphere.9)) mat64 <.8).7). The value of answer is (a) 12.7.6. plot(distance~stretch. volume <.pch=1. par(mfrow = c(1.6. 4:6].3).87 painters pressure ships Data Set from Package lattice barley 12.data <.9 1. data = snow) hist(snow$snow.0.10^seq(-1.9 2.3:20.9)).rep(6.
summary(cpus) Comment on ranges of values. Cars93[Cars93$Type=="Small"|Cars93$Type=="Sporty". (a) Use summary(airquality) to get this information.] (c) airquality$Wind[airquality$Ozone > quantile(airquality$Ozone. whether distributions seem skew. mean(snow$snow. .75)] 7. 1:9) 6. mtcars6<-mtcars[mtcars$cyl==6.] . etc. (b) airquality[airquality$Ozone == max(airquality$Ozone).88 4.seq(1. 10.10.cover[seq(2.cover[seq(1.2)]) mean(snow$snow.9).] 11.9. summary(attitude).9)) or rep(1:9. rep(seq(1.2)]) 9.