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