statistics

Nicholas J. Horton Randall Pruim Daniel T. Kaplan

A Student's Guide to

R

Project MOSAIC

2 horton, kaplan, pruim

Copyright (c) 2015 by Nicholas J. Horton, Randall Pruim, & Daniel Kaplan.

Edition 1.2, November 2015

This material is copyrighted by the authors under a Creative Commons Attribution 3.0 Unported License. You are free to Share (to copy, distribute and transmit the work) and to Remix (to adapt the work) if you attribute our work. More detailed information about the licensing is available at this web page: http: //www.mosaic-web.org/go/teachingRlicense.html.

Cover Photo: Maya Hanna.

Contents

1 Introduction 13

2 Getting Started with RStudio 15

3 One Quantitative Variable 27

4 One Categorical Variable 39

5 Two Quantitative Variables 45

6 Two Categorical Variables 55

7 Quantitative Response, Categorical Predictor 61

8 Categorical Response, Quantitative Predictor 69

9 Survival Time Outcomes 73

4 horton, kaplan, pruim

10 More than Two Variables 75

11 Probability Distributions & Random Variables 83

12 Power Calculations 89

13 Data Management 93

14 Health Evaluation (HELP) Study 107

15 Exercises and Problems 111

16 Bibliography 115

17 Index 117

About These Notes

We present an approach to teaching introductory and in- termediate statistics courses that is tightly coupled with computing generally and with R and RStudio in particular. These activities and examples are intended to highlight a modern approach to statistical education that focuses on modeling, resampling based inference, and multivari- ate graphical techniques. A secondary goal is to facilitate computing with data through use of small simulation studies and appropriate statistical analysis workflow. This follows the philosophy outlined by Nolan and Temple Lang1. The importance of modern computation in statis- 1 D. Nolan and D. Temple Lang.

Computing in the statistics curriculum. The American Statistician, 64(2):97–107, 2010

tics education is a principal component of the recently adopted American Statistical Association’s curriculum guidelines2.

2 ASA Undergraduate Guide- lines Workgroup. 2014 cur- riculum guidelines for under- graduate programs in statisti- cal science. Technical report, American Statistical Associa- tion, November 2014. http: //www.amstat.org/education/ curriculumguidelines.cfm

Throughout this book (and its companion volumes), we introduce multiple activities, some appropriate for an introductory course, others suitable for higher levels, that demonstrate key concepts in statistics and modeling while also supporting the core material of more tradi- tional courses.

A Work in Progress Caution!

Despite our best efforts, you WILL find bugs both in this document and in our code. Please let us know when you encounter them so we can call in the exterminators.

These materials were developed for a workshop entitled Teaching Statistics Using R prior to the 2011 United States Conference on Teaching Statistics and revised for US- COTS 2011, USCOTS 2013, eCOTS 2014, ICOTS 9, and USCOTS 2015. We organized these workshops to help instructors integrate R (as well as some related technolo- gies) into statistics courses at all levels. We received great feedback and many wonderful ideas from the participants and those that we’ve shared this with since the work- shops.

6 horton, kaplan, pruim

Consider these notes to be a work in progress. We ap- preciate any feedback you are willing to share as we con- tinue to work on these materials and the accompanying mosaic package. Drop us an email at pis@mosaic-web. org with any comments, suggestions, corrections, etc.

Updated versions will be posted at http://mosaic-web. org.

Two Audiences

We initially developed these materials for instructors of statistics at the college or university level. Another audi- ence is the students these instructors teach. Some of the sections, examples, and exercises are written with one or the other of these audiences more clearly at the forefront. This means that

1. Some of the materials can be used essentially as is with students.

2. Some of the materials aim to equip instructors to de- velop their own expertise in R and RStudio to develop their own teaching materials.

Although the distinction can get blurry, and what works “as is" in one setting may not work “as is" in an- other, we’ll try to indicate which parts fit into each cate- gory as we go along.

R, RStudio and R Packages

R can be obtained from http://cran.r-project.org/. Download and installation are quite straightforward for Mac, PC, or linux machines.

RStudio is an integrated development environment (IDE) that facilitates use of R for both novice and expert users. We have adopted it as our standard teaching en- vironment because it dramatically simplifies the use of R for instructors and for students. RStudio can be installed

More Info

Several things we use that can be done only in RStudio, for instance manipulate() or RStudio’s integrated support for reproducible research).

as a desktop (laptop) application or as a server applica- tion that is accessible to users via the Internet. RStudio server version works

well with starting students. All they need is a web browser, avoiding any potential prob- lems with oddities of students’ individual computers.

In addition to R and RStudio, we will make use of sev- eral packages that need to be installed and loaded sep- arately. The mosaic package (and its dependencies) will

a student’s guide to r 7

be used throughout. Other packages appear from time to time as well.

Marginal Notes

Marginal notes appear here and there. Sometimes these Have a great suggestion for a marginal note? Pass it along.are side comments that we wanted to say, but we didn’t

want to interrupt the flow to mention them in the main text. Others provide teaching tips or caution about traps, pitfalls and gotchas.

What’s Ours Is Yours – To a Point

This material is copyrighted by the authors under a Cre- ative Commons Attribution 3.0 Unported License. You are free to Share (to copy, distribute and transmit the work) and to Remix (to adapt the work) if you attribute our work. More detailed information about the licensing is available at this web page: http://www.mosaic-web. org/go/teachingRlicense.html. Digging Deeper

If you know LATEX as well as R, then knitr provides a nice solution for mixing the two. We used this system to produce this book. We also use it for our own research and to intro- duce upper level students to reproducible analysis methods. For beginners, we introduce knitr with RMarkdown, which produces PDF, HTML, or Word files using a simpler syntax.

Document Creation

This document was created on November 15, 2015, using

• knitr, version 1.11

• mosaic, version 0.12.9003

• mosaicData, version 0.12.9003

• R version 3.2.2 (2015-08-14)

Inevitably, each of these will be updated from time to time. If you find that things look different on your com- puter, make sure that your version of R and your pack- ages are up to date and check for a newer version of this document.

Project MOSAIC

This book is a product of Project MOSAIC, a community of educators working to develop new ways to introduce mathematics, statistics, computation, and modeling to students in colleges and universities.

The goal of the MOSAIC project is to help share ideas and resources to improve teaching, and to develop a cur- ricular and assessment infrastructure to support the dis- semination and evaluation of these approaches. Our goal is to provide a broader approach to quantitative stud- ies that provides better support for work in science and technology. The project highlights and integrates diverse aspects of quantitative work that students in science, tech- nology, and engineering will need in their professional lives, but which are today usually taught in isolation, if at all.

In particular, we focus on:

Modeling The ability to create, manipulate and investigate useful and informative mathematical representations of a real-world situations.

Statistics The analysis of variability that draws on our ability to quantify uncertainty and to draw logical in- ferences from observations and experiment.

Computation The capacity to think algorithmically, to manage data on large scales, to visualize and inter- act with models, and to automate tasks for efficiency, accuracy, and reproducibility.

Calculus The traditional mathematical entry point for col- lege and university students and a subject that still has the potential to provide important insights to today’s students.

10 horton, kaplan, pruim

Drawing on support from the US National Science Foundation (NSF DUE-0920350), Project MOSAIC sup- ports a number of initiatives to help achieve these goals, including:

Faculty development and training opportunities, such as the USCOTS 2011, USCOTS 2013, eCOTS 2014, and ICOTS 9 workshops on Teaching Statistics Using R and RStu- dio, our 2010 Project MOSAIC kickoff workshop at the Institute for Mathematics and its Applications, and our Modeling: Early and Often in Undergraduate Calculus AMS PREP workshops offered in 2012, 2013, and 2015.

M-casts, a series of regularly scheduled webinars, de- livered via the Internet, that provide a forum for in- structors to share their insights and innovations and to develop collaborations to refine and develop them. Recordings of M-casts are available at the Project MO- SAIC web site, http://mosaic-web.org.

The construction of syllabi and materials for courses that teach MOSAIC topics in a better integrated way. Such courses and materials might be wholly new construc- tions, or they might be incremental modifications of existing resources that draw on the connections be- tween the MOSAIC topics.

More details can be found at http://www.mosaic-web. org. We welcome and encourage your participation in all of these initiatives.

Computational Statistics

There are at least two ways in which statistical software can be introduced into a statistics course. In the first ap- proach, the course is taught essentially as it was before the introduction of statistical software, but using a com- puter to speed up some of the calculations and to prepare higher quality graphical displays. Perhaps the size of the data sets will also be increased. We will refer to this approach as statistical computation since the computer serves primarily as a computational tool to replace pencil- and-paper calculations and drawing plots manually.

In the second approach, more fundamental changes in the course result from the introduction of the computer. Some new topics are covered, some old topics are omit- ted. Some old topics are treated in very different ways, and perhaps at different points in the course. We will re- fer to this approach as computational statistics because the availability of computation is shaping how statistics is done and taught. Computational statistics is a key com- ponent of data science, defined as the ability to use data to answer questions and communicate those results.

Students need to see aspects of computation and data science early and often to develop deeper skills. Establishing precursors in introductory courses help them get started.

In practice, most courses will incorporate elements of both statistical computation and computational statistics, but the relative proportions may differ dramatically from course to course. Where on the spectrum a course lies will be depend on many factors including the goals of the course, the availability of technology for student use, the perspective of the text book used, and the comfort-level of the instructor with both statistics and computation.

Among the various statistical software packages avail- able, R is becoming increasingly popular. The recent addi- tion of RStudio has made R both more powerful and more accessible. Because R and RStudio are free, they have be- come widely used in research and industry. Training in R

12 horton, kaplan, pruim

and RStudio is often seen as an important additional skill that a statistics course can develop. Furthermore, an in- creasing number of instructors are using R for their own statistical work, so it is natural for them to use it in their teaching as well. At the same time, the development of R and of RStudio (an optional interface and integrated de- velopment environment for R) are making it easier and easier to get started with R.

Information about the mosaic package, including vignettes demonstrating features and supplementary materials (such as this book) can be found at https://cran.r-project.org/ web/packages/mosaic.

We developed the mosaic R package (available on CRAN) to make certain aspects of statistical computation and computational statistics simpler for beginners, with- out limiting their ability to use more advanced features of the language. The mosaic package includes a modelling approach that uses the same general syntax to calculate descriptive statistics, create graphics, and fit linear mod- els.

1 Introduction

In this reference book, we briefly review the commands and functions needed to analyze data from introductory and second courses in statistics. This is intended to com- plement the Start Teaching with R and Start Modeling with R books.

Most of our examples will use data from the HELP (Health Evaluation and Linkage to Primary Care) study: a randomized clinical trial of a novel way to link at-risk subjects with primary care. More information on the dataset can be found in chapter 14.

Since the selection and order of topics can vary greatly from textbook to textbook and instructor to instructor, we have chosen to organize this material by the kind of data being analyzed. This should make it straightforward to find what you are looking for. Some data management skills are needed by students1. A basic introduction to 1 N.J. Horton, B.S. Baumer, and

H. Wickham. Setting the stage for data science: integration of data management skills in introductory and second courses in statistics (http: //arxiv.org/abs/1401.3269). CHANCE, 28(2):40–50, 2015

key idioms is provided in Chapter 13. This work leverages initiatives undertaken by Project

MOSAIC (http://www.mosaic-web.org), an NSF-funded effort to improve the teaching of statistics, calculus, sci- ence and computing in the undergraduate curriculum. In particular, we utilize the mosaic package, which was written to simplify the use of R for introductory statis- tics courses, and the mosaicData package which includes a number of data sets. A short summary of the R com- mands needed to teach introductory statistics can be found in the mosaic package vignette: https://cran. r-project.org/web/packages/mosaic.

Other related resources from Project MOSAIC may be helpful, including an annotated set of examples from the sixth edition of Moore, McCabe and Craig’s Introduction to the Practice of Statistics2 (see http://www.amherst.edu/ 2 D. S. Moore and G. P. McCabe.

Introduction to the Practice of Statistics. W.H.Freeman and Company, 6th edition, 2007

14 horton, kaplan, pruim

~nhorton/ips6e), the second and third editions of the Sta- tistical Sleuth3 (see http://www.amherst.edu/~nhorton/ 3 F. Ramsey and D. Schafer.

Statistical Sleuth: A Course in Methods of Data Analysis. Cen- gage, 2nd edition, 2002

sleuth), and Statistics: Unlocking the Power of Data by Lock et al (see https://github.com/rpruim/Lock5withR).

To use a package within R, it must be installed (one time), and loaded (each session). The mosaic package can be installed using the following commands:

> install.packages("mosaic") # note the quotation marks

The # character is a comment in R, and all text after that

RStudio features a simplified package installation tab (in the bottom right panel).

on the current line is ignored. Once the package is installed (one time only), it can be

loaded by running the command:

> require(mosaic)

The knitr/LATEX system allows experienced users to combine R and LATEX in the same docu- ment. The reward for learning this more complicated system is much finer control over the output format. But RMarkdown is much easier to learn and is adequate even for professional- level work.

Using Markdown or knitr/LATEX requires that the markdown package be installed.

The RMarkdown system provides a simple markup language and renders the results in PDF, Word, or HTML. This allows students to undertake their analyses using a workflow that facilitates “reproducibility” and avoids cut and paste errors.

We typically introduce students to RMarkdown very early, requiring students to use it for assignments and reports4. 4 B.S. Baumer, M. Çetinkaya

Rundel, A. Bray, L. Loi, and N. J. Horton. R Markdown: Integrating a reproducible analysis tool into introductory statistics. Technology Innovations in Statistics Education, 8(1):281– 283, 2014

2 Getting Started with RStudio

RStudio is an integrated development environment (IDE) for R that provides an alternative interface to R that has several advantages over other the default R interfaces:

A series of getting started videos are available at http://www.amherst.edu/ ~nhorton/rstudio.

• RStudio runs on Mac, PC, and Linux machines and pro- vides a simplified interface that looks and feels identical on all of them. The default interfaces for R are quite different on the various platforms. This is a distractor for students and adds an extra layer of support responsibility for the instructor.

• RStudio can run in a web browser. In addition to stand-alone desktop versions, RStudio can be set up as a server application that is accessed via the internet. The web interface is nearly identical to the desktop version. As with other web services, users login to Caution!

The desktop and server version of RStudio are so similar that if you run them both, you will have to pay careful attention to make sure you are working in the one you intend to be working in.

access their account. If students logout and login in again later, even on a different machine, their session is restored and they can resume their analysis right where they left off. With a little advanced set up, in- structors can save the history of their classroom R use and students can load those history files into their own environment. Note

Using RStudio in a browser is like Facebook for statistics. Each time the user returns, the previous session is restored and they can resume work where they left off. Users can login from any device with internet access.

• RStudio provides support for reproducible research. RStudio makes it easy to include text, statistical

analysis (R code and R output), and graphical displays all in the same document. The RMarkdown system provides a simple markup language and renders the results in HTML. The knitr/LATEX system allows users

16 horton, kaplan, pruim

to combine R and LATEX in the same document. The reward for learning this more complicated system is much finer control over the output format. Depending on the level of the course, students can use either of these for homework and projects. To use Markdown or

knitr/LATEX requires that the knitr package be installed on your system.• RStudio provides an integrated support for editing and

executing R code and documents.

• RStudio provides some useful functionality via a graph- ical user interface.

RStudio is not a GUI for R, but it does provide a GUI that simplifies things like installing and updating packages; monitoring, saving and loading environ- ments; importing and exporting data; browsing and exporting graphics; and browsing files and documenta- tion.

• RStudio provides access to the manipulate package. The manipulate package provides a way to create

simple interactive graphical applications quickly and easily.

While one can certainly use R without using RStudio, RStudio makes a number of things easier and we highly recommend using RStudio. Furthermore, since RStudio is in active development, we fully expect more useful fea- tures in the future.

We primarily use an online version of RStudio. RStudio is a innovative and powerful interface to R that runs in a web browser or on your local machine. Running in the browser has the advantage that you don’t have to install or configure anything. Just login and you are good to go. Furthermore, RStudio will “remember” what you were doing so that each time you login (even on a different machine) you can pick up right where you left off. This is “R in the cloud" and works a bit like GoogleDocs or Facebook for R.

R can also be obtained from http://cran.r-project. org/. Download and installation are pretty straightfor- ward for Mac, PC, or Linux machines. RStudio is available from http://www.rstudio.org/.

a student’s guide to r 17

2.1 Connecting to an RStudio server

RStudio servers have been set up at a number of schools to facilitate cloud-based computing.

RStudio servers have been in- stalled at many institutions. More details about (free) aca- demic licenses for RStudio Server Pro as well as setup instructions can be found at http://www.rstudio.com/ resources/faqs under the Academic tab.

Once you connect to the server, you should see a login screen:

The RStudio server doesn’t tend to work well with Internet Explorer.

Once you authenticate, you should see the RStudio interface:

Notice that RStudio divides its world into four panels. Several of the panels are further subdivided into multi-

18 horton, kaplan, pruim

ple tabs. Which tabs appear in which panels can be cus- tomized by the user.

R can do much more than a simple calculator, and we will introduce additional features in due time. But per- forming simple calculations in R is a good way to begin learning the features of RStudio.

Commands entered in the Console tab are immediately executed by R. A good way to familiarize yourself with the console is to do some simple calculator-like compu- tations. Most of this will work just like you would expect from a typical calculator. Try typing the following com- mands in the console panel.

> 5 + 3

[1] 8

> 15.3 * 23.4

[1] 358.02

> sqrt(16) # square root

[1] 4

This last example demonstrates how functions are called within R as well as the use of comments. Com- ments are prefaced with the # character. Comments can be very helpful when writing scripts with multiple com- mands or to annotate example code for your students.

You can save values to named variables for later reuse.

It’s probably best to settle on using one or the other of the right-to-left assignment opera- tors rather than to switch back and forth. We prefer the arrow operator because it represents visually what is happening in an assignment and because it makes a clear distinction be- tween the assignment operator, the use of = to provide values to arguments of functions, and the use of == to test for equality.

> product = 15.3 * 23.4 # save result

> product # display the result

[1] 358.02

> product <- 15.3 * 23.4 # <- can be used instead of =

> product

[1] 358.02

Once variables are defined, they can be referenced in other operations and functions.

a student’s guide to r 19

> 0.5 * product # half of the product

[1] 179.01

> log(product) # (natural) log of the product

[1] 5.880589

> log10(product) # base 10 log of the product

[1] 2.553907

> log2(product) # base 2 log of the product

[1] 8.483896

> log(product, base=2) # base 2 log of the product, another way

[1] 8.483896

The semi-colon can be used to place multiple com- mands on one line. One frequent use of this is to save and print a value all in one go:

> product <- 15.3 * 23.4; product # save result and show it

[1] 358.02

2.1.1 Version information

At times it may be useful to check what version of the mosaic package, R, and RStudioyou are using. Running sessionInfo() will display information about the version of R and packages that are loaded and RStudio.Version() will provide information about the version of RStudio.

> sessionInfo()

R version 3.2.2 (2015-08-14)

Platform: x86_64-apple-darwin13.4.0 (64-bit)

Running under: OS X 10.10.5 (Yosemite)

locale:

[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

20 horton, kaplan, pruim

attached base packages:

[1] grid stats graphics grDevices utils datasets

[7] methods base

other attached packages:

[1] mosaic_0.12.9003 mosaicData_0.9.9001 car_2.1-0

[4] ggplot2_1.0.1 dplyr_0.4.3 lattice_0.20-33

[7] knitr_1.11

loaded via a namespace (and not attached):

[1] Rcpp_0.12.1 magrittr_1.5 splines_3.2.2

[4] MASS_7.3-45 munsell_0.4.2 colorspace_1.2-6

[7] R6_2.1.1 ggdendro_0.1-17 minqa_1.2.4

[10] highr_0.5.1 stringr_1.0.0 plyr_1.8.3

[13] tools_3.2.2 nnet_7.3-11 parallel_3.2.2

[16] pbkrtest_0.4-2 nlme_3.1-122 gtable_0.1.2

[19] mgcv_1.8-9 quantreg_5.19 DBI_0.3.1

[22] MatrixModels_0.4-1 lme4_1.1-10 assertthat_0.1

[25] digest_0.6.8 Matrix_1.2-2 gridExtra_2.0.0

[28] nloptr_1.0.4 reshape2_1.4.1 formatR_1.2.1

[31] evaluate_0.8 stringi_1.0-1 scales_0.3.0

[34] SparseM_1.7 proto_0.3-10

2.2 Working with Files

2.2.1 Working with R Script Files

As an alternative, R commands can be stored in a file. RStudio provides an integrated editor for editing these files and facilitates executing some or all of the com- mands. To create a file, select File, then New File, then R Script from the RStudio menu. A file editor tab will open in the Source panel. R code can be entered here, and but- tons and menu items are provided to run all the code (called sourcing the file) or to run the code on a single line or in a selected section of the file.

2.2.2 Working with RMarkdown, and knitr/LATEX

A third alternative is to take advantage of RStudio’s sup- port for reproducible research. If you already know LATEX,

a student’s guide to r 21

you will want to investigate the knitr/LATEX capabili- ties. For those who do not already know LATEX, the sim- pler RMarkdown system provides an easy entry into the world of reproducible research methods. It also provides a good facility for students to create homework and re- ports that include text, R code, R output, and graphics.

To create a new RMarkdown file, select File, then New File, then RMarkdown. The file will be opened with a short template document that illustrates the mark up language.

The mosaic package includes two useful RMarkdown templates for getting started: fancy includes bells and whistles (and is intended to give an overview of features), while plain is useful as a starting point for a new analy- sis. These are accessed using the Template option when creating a new RMarkdown file.

22 horton, kaplan, pruim

Click on the Knit button to convert to an HTML, PDF, or Word file.

This will generate a formatted version of the docu- ment.

a student’s guide to r 23

There is a button (marked with a question mark) which provides a brief description of the supported markup commands. The RStudio web site includes more extensive tutorials on using RMarkdown. Caution!

RMarkdown, and knitr/LATEX files do not have access to the console environment, so the code in them must be self- contained.

It is important to remember that unlike R scripts, which are executed in the console and have access to the console environment, RMarkdown and knitr/LATEX files do not have access to the console environment This is a good feature because it forces the files to be self- contained, which makes them transferable and respects good reproducible research practices. But beginners, es- pecially if they adopt a strategy of trying things out in the console and copying and pasting successful code from the console to their file, will often create files that are incom- plete and therefore do not compile correctly.

2.3 The Other Panels and Tabs

2.3.1 The History Tab

As commands are entered in the console, they appear in the History tab. These histories can be saved and loaded, there is a search feature to locate previous commands, and individual lines or sections can be transferred back to the console. Keeping the History tab open will allow you to go back and see the previous several commands. This can be especially useful when commands produce a fair amount of output and so scroll off the screen rapidly.

24 horton, kaplan, pruim

2.3.2 Communication between tabs

RStudio provides several ways to move R code between tabs. Pressing the Run button in the editing panel for an R script or RMarkdown or other file will copy lines of code into the Console and run them.

2.3.3 The Files Tab

The Files tab provides a simple file manager. It can be navigated in familiar ways and used to open, move, re- name, and delete files. In the browser version of RStudio, the Files tab also provides a file upload utility for moving files from the local machine to the server. In RMarkdown and knitr files one can also run the code in a particular chunk or in all of the chunks in a file. Each of these fea- tures makes it easy to try out code “live” while creating a document that keeps a record of the code.

In the reverse direction, code from the history can be copied either back into the console to run them again (perhaps after editing) or into one of the file editing tabs for inclusion in a file.

2.3.4 The Help Tab

The Help tab is where RStudio displays R help files. These can be searched and navigated in the Help tab. You can also open a help file using the ? operator in the console. For example the following command will provide the help file for the logarithm function.

> ?log

2.3.5 The Environment Tab

The Environment tab shows the objects available to the console. These are subdivided into data, values (non- dataframe, non-function objects) and functions. The broom icon can be used to remove all objects from the en- vironment, and it is good to do this from time to time, es- pecially when running in RStudio server or if you choose

a student’s guide to r 25

to save the environment when shutting down RStudio since in these cases objects can stay in the environment essentially indefinitely.

2.3.6 The Plots Tab

Plots created in the console are displayed in the Plots tab. For example, the following commands display the num- ber of births in the United States for each day in 1978.

> # this will make lattice graphics available to the session

> # as well as the Births78 dataset

> require(mosaic)

> xyplot(births ~ dayofyear, data=Births78)

dayofyear

bi rt

hs

7000

8000

9000

10000

0 100 200 300

● ●

●● ● ●

●

●

●●● ●●

●

●

●● ●● ●

●

●

● ●

● ●

●

●

●

●●●●●

● ●

●●●● ●

● ●

●

●

● ● ●

●

●

● ●●● ●

●

●

●

●

● ●●

● ●

● ●

●●

●

● ●

●

● ●● ●

● ●

●

●

● ●●

● ●

● ●● ●

●

●

●

● ●

●● ●

●

●

● ●

●● ●

●

●

●

●

● ● ●

●

●

●●● ● ●

●

●

●

●

●

●

●

●

●

● ●●

●●

●

●

●● ● ●●

●

●

●●●●

●

● ●●

●● ● ●

●

●

●

● ● ● ●

●

●

●●●● ●

●

●

●●

● ● ●

●●

●

● ●●

●

● ●

●

●

●

●●

●

●

●● ● ●

●

●

●

●

●

●

●●

●

●

●

● ●● ●

● ●

● ● ● ●●

● ●

● ● ● ●●

●

●

●

●

● ●●

●

●

● ●

● ●●

●

●

●

● ● ● ●

● ●●

●

●

●●

●

●

● ● ●●●

●

●

●

● ●● ●

●●

●●● ● ●

●●

●

●

●●●

● ●

●

●

●●●

●

●

●

●

●● ●

●

●

● ● ● ●●

● ●

●●●●

●

● ●

● ● ●● ●

●

●

●

●

●● ●

●

●

●

●

●

●

●

● ●

●●

●● ●

●

●

●●

●●●

● ●

●

●

●●

●

●

●

●

●

●

●

●

●●●

●

● ● ●

●

●

From the Plots tab, you can navigate to previous plots and also export plots in various formats after interactively resizing them.

2.3.7 The Packages Tab

Much of the functionality of R is located in packages, many of which can be obtained from a central clearing house called CRAN (Comprehensive R Archive Network). The Packages tab facilitates installing and loading pack- ages. It will also allow you to search for packages that have been updated since you installed them.

3 One Quantitative Variable

3.1 Numerical summaries

R includes a number of commands to numerically sum- marize variables. These include the capability of calculat- ing the mean, standard deviation, variance, median, five number summary, interquartile range (IQR) as well as ar- bitrary quantiles. We will illustrate these using the CESD (Center for Epidemiologic Studies–Depression) measure of depressive symptoms (which takes on values between 0 and 60, with higher scores indicating more depressive symptoms).

To improve the legibility of output, we will also set the default number of digits to display to a more reasonable level (see ?options() for more configuration possibilities).

> require(mosaic)

> require(mosaicData)

> options(digits=3)

> mean(~ cesd, data=HELPrct)

[1] 32.8

Note that the mean() function in the mosaic package supports a formula interface common to lattice graphics and linear models (e.g., lm()). The mosaic package pro- vides many other functions that use the same notation, which we will be using throughout this document.

Digging Deeper

If you have not seen the for- mula notation before, the com- panion book, Start Teaching with R provides a detailed presen- tation. Start Modeling with R, another companion book, de- tails the relationship between the process of modeling and the formula notation.

The same output could be created using the following commands (though we will use the MOSAIC versions when available).

28 horton, kaplan, pruim

> with(HELPrct, mean(cesd))

[1] 32.8

> mean(HELPrct$cesd)

[1] 32.8

Similar functionality exists for other summary statistics.

> sd(~ cesd, data=HELPrct)

[1] 12.5

> sd(~ cesd, data=HELPrct)^2

[1] 157

> var(~ cesd, data=HELPrct)

[1] 157

It is also straightforward to calculate quantiles of the distribution.

> median(~ cesd, data=HELPrct)

[1] 34

By default, the quantile() function displays the quar- tiles, but can be given a vector of quantiles to display. Caution!

Not all commands have been upgraded to support the for- mula interface. For these functions, variables within dataframes must be accessed using with() or the $ operator.

> with(HELPrct, quantile(cesd))

0% 25% 50% 75% 100%

1 25 34 41 60

> with(HELPrct, quantile(cesd, c(.025, .975)))

2.5% 97.5%

6.3 55.0

a student’s guide to r 29

Finally, the favstats() function in the mosaic package provides a concise summary of many useful statistics.

> favstats(~ cesd, data=HELPrct)

min Q1 median Q3 max mean sd n missing

1 25 34 41 60 32.8 12.5 453 0

3.2 Graphical summaries

The histogram() function is used to create a histogram. Here we use the formula interface (as discussed in the Start Modeling with R book) to specify that we want a his- togram of the CESD scores.

> histogram(~ cesd, data=HELPrct)

cesd

D en

si ty

0.00

0.01

0.02

0.03

0 20 40 60

We can use the width() and center() options to con- trol the location of the bins.

> histogram(~ cesd, width=5, center=2.5, data=HELPrct)

30 horton, kaplan, pruim

cesd

D en

si ty

0.00

0.01

0.02

0.03

0 20 40 60

In the HELPrct dataset, approximately one quarter of the subjects are female.

> tally(~ sex, data=HELPrct)

female male

107 346

> tally(~ sex, format="percent", data=HELPrct)

female male

23.6 76.4

It is straightforward to restrict our attention to just the female subjects. If we are going to do many things with a subset of our data, it may be easiest to make a new dataframe containing only the cases we are interested in. The filter() function in the dplyr package can be used to generate a new dataframe containing just the women or just the men (see also section 13.5). Once this is created, the the stem() function is used to create a stem and leaf plot. Caution!

Note that the tests for equality use two equal signs> female <- filter(HELPrct, sex=='female')

> male <- filter(HELPrct, sex=='male')

> with(female, stem(cesd))

The decimal point is 1 digit(s) to the right of the |

a student’s guide to r 31

0 | 3

0 | 567

1 | 3

1 | 555589999

2 | 123344

2 | 66889999

3 | 0000233334444

3 | 5556666777888899999

4 | 00011112222334

4 | 555666777889

5 | 011122222333444

5 | 67788

6 | 0

Subsets can also be generated and used “on the fly" (this time including an overlaid normal density):

> histogram(~ cesd, fit="normal",

data=filter(HELPrct, sex=='female'))

cesd

D en

si ty

0.000

0.005

0.010

0.015

0.020

0.025

0.030

0 20 40 60

Alternatively, we can make side-by-side plots to com- pare multiple subsets.

> histogram(~ cesd | sex, data=HELPrct)

32 horton, kaplan, pruim

cesd

D en

si ty

0.00

0.01

0.02

0.03

0.04

0 20 40 60

female

0 20 40 60

male

The layout can be rearranged.

> histogram(~ cesd | sex, layout=c(1, 2), data=HELPrct)

cesd

D en

si ty

0.00 0.01 0.02 0.03 0.04

0 20 40 60

female 0.00 0.01 0.02 0.03 0.04

male

We can control the number of bins in a number of ways. These can be specified as the total number.

> histogram(~ cesd, nint=20, data=female)

a student’s guide to r 33

cesd

D en

si ty

0.00

0.01

0.02

0.03

0 10 20 30 40 50 60

The width of the bins can be specified.

> histogram(~ cesd, width=2, data=female)

cesd

D en

si ty

0.00

0.01

0.02

0.03

0.04

0 10 20 30 40 50 60

The dotPlot() function is used to create a dotplot for a smaller subset of subjects (homeless females). We also demonstrate how to change the x-axis label.

> dotPlot(~ cesd, xlab="CESD score",

data=filter(HELPrct, (sex=="female") & (homeless=="homeless")))

CESD score

C ou

nt

0

5

10

20 40 60

● ● ● ● ● ●

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

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

● ● ● ● ● ● ● ●

●

34 horton, kaplan, pruim

3.3 Density curves

Density plots are also sensi- tive to certain choices. If your density plot is too jagged or too smooth, try changing the adjust argument: larger than 1

for smoother plots, less than 1

for more jagged plots.

One disadvantage of histograms is that they can be sensi- tive to the choice of the number of bins. Another display to consider is a density curve.

Here we adorn a density plot with some additions to demonstrate how to build up a graphic for pedagogical purposes. We add some text, a superimposed normal density as well as a vertical line. A variety of line types and colors can be specified, as well as line widths. Digging Deeper

The plotFun() function can also be used to annotate plots (see section 10.2.1).> densityplot(~ cesd, data=female)

> ladd(grid.text(x=0.2, y=0.8, 'only females'))

> ladd(panel.mathdensity(args=list(mean=mean(cesd),

sd=sd(cesd)), col="red"), data=female)

> ladd(panel.abline(v=60, lty=2, lwd=2, col="grey"))

cesd

D en

si ty

0.000

0.005

0.010

0.015

0.020

0.025

0.030

0 20 40 60

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

only females

3.4 Frequency polygons

A third option is a frequency polygon, where the graph is created by joining the midpoints of the top of the bars of a histogram.

a student’s guide to r 35

> freqpolygon(~ cesd, data=female)

cesd

D en

si ty

0.01

0.02

0.03

10 20 30 40 50 60

3.5 Normal distributions

x is for eXtra.The most famous density curve is a normal distribution. The xpnorm() function displays the probability that a ran- dom variable is less than the first argument, for a normal distribution with mean given by the second argument and standard deviation by the third. More information about probability distributions can be found in section 11.

> xpnorm(1.96, mean=0, sd=1)

If X ~ N(0,1), then

P(X <= 1.96) = P(Z <= 1.96) = 0.975

P(X > 1.96) = P(Z > 1.96) = 0.025

[1] 0.975

de ns

ity

0.1

0.2

0.3

0.4

0.5

−2 0 2

1.96 (z=1.96)

0.975 0.025

36 horton, kaplan, pruim

3.6 Inference for a single sample

We can calculate a 95% confidence interval for the mean CESD score for females by using a t-test:

> t.test(~ cesd, data=female)

One Sample t-test

data: data$cesd

t = 30, df = 100, p-value <2e-16

alternative hypothesis: true mean is not equal to 0

95 percent confidence interval:

34.4 39.4

sample estimates:

mean of x

36.9

> confint(t.test(~ cesd, data=female))

mean of x lower upper level

1 36.9 34.4 39.4 0.95

Digging Deeper

More details and examples can be found in the mosaic package Resampling Vignette.

But it’s also straightforward to calculate this using a bootstrap. The statistic that we want to resample is the mean.

> mean(~ cesd, data=female)

[1] 36.9

One resampling trial can be carried out: Here we sample with re- placement from the original dataframe, creating a resam- pled dataframe with the same number of rows.

> mean(~ cesd, data=resample(female))

[1] 37.7

Even though a single trial is of little use, it’s smart having students do the calculation to show that they are (usually!) getting a different result than without resampling.

Another will yield different results:

> mean(~ cesd, data=resample(female))

[1] 34.9

a student’s guide to r 37

Now conduct 1000 resampling trials, saving the results in an object called trials:

> trials <- do(1000) * mean(~ cesd, data=resample(female))

> head(trials, 3)

mean

1 36.7

2 35.5

3 35.9

> qdata(~ mean, c(.025, .975), data=trials)

quantile p

2.5% 34.3 0.025

97.5% 39.4 0.975

4 One Categorical Variable

4.1 Numerical summaries

Digging Deeper

The Start Teaching with R com- panion book introduces the for- mula notation used throughout this book. See also Start Teach- ing with R for the connections to statistical modeling.

The tally() function can be used to calculate counts, percentages and proportions for a categorical variable.

> tally(~ homeless, data=HELPrct)

homeless housed

209 244

> tally(~ homeless, margins=TRUE, data=HELPrct)

homeless housed Total

209 244 453

> tally(~ homeless, format="percent", data=HELPrct)

homeless housed

46.1 53.9

> tally(~ homeless, format="proportion", data=HELPrct)

homeless housed

0.461 0.539

40 horton, kaplan, pruim

4.2 The binomial test

An exact confidence interval for a proportion (as well as a test of the null hypothesis that the population propor- tion is equal to a particular value [by default 0.5]) can be calculated using the binom.test() function. The standard binom.test() requires us to tabulate.

> binom.test(209, 209 + 244)

data: 209 out of 209 + 244

number of successes = 200, number of trials = 500, p-value =

0.1

alternative hypothesis: true probability of success is not equal to 0.5

95 percent confidence interval:

0.415 0.509

sample estimates:

probability of success

0.461

The mosaic package provides a formula interface that avoids the need to pre-tally the data.

> result <- binom.test(~ (homeless=="homeless"), data=HELPrct)

> result

data: HELPrct$(homeless == "homeless") [with success = TRUE]

number of successes = 200, number of trials = 500, p-value =

0.1

alternative hypothesis: true probability of success is not equal to 0.5

95 percent confidence interval:

0.415 0.509

sample estimates:

probability of success

0.461

As is generally the case with commands of this sort, there are a number of useful quantities available from the object returned by the function.

a student’s guide to r 41

> names(result)

[1] "statistic" "parameter" "p.value" "conf.int"

[5] "estimate" "null.value" "alternative" "data.name"

These can be extracted using the $ operator or an extrac- tor function. For example, the user can extract the confi- dence interval or p-value.

> result$statistic

number of successes

209

> confint(result)

probability of success lower upper level

1 0.461 0.415 0.509 0.95

> pval(result)

p.value

0.11

Digging Deeper

Most of the objects in R have a print() method. So when we get result, what we are seeing displayed in the console is print(result). There may be a good deal of additional information lurking inside the object itself.

In some situations, such as graphics, the object is returned invisibly, so nothing prints. That avoids your having to look at a long printout not intended for human consumption. You can still assign the returned object to a variable and pro- cess it later, even if nothing shows up on the screen. This is sometimes helpful for lattice graphics functions.

4.3 The proportion test

A similar interval and test can be calculated using the function prop.test(). Here is a count of the number of people at each of the two levels of homeless

> tally(~ homeless, data=HELPrct)

homeless housed

209 244

The prop.test() function will carry out the calcula- tions of the proportion test and report the result.

> prop.test(~ (homeless=="homeless"), correct=FALSE, data=HELPrct)

42 horton, kaplan, pruim

1-sample proportions test without continuity correction

data: HELPrct$(homeless == "homeless") [with success = TRUE]

X-squared = 3, df = 1, p-value = 0.1

alternative hypothesis: true p is not equal to 0.5

95 percent confidence interval:

0.416 0.507

sample estimates:

p

0.461

In this statement, prop.test is examining the homeless variable in the same way that tally() would. prop.test() More Info

We write homeless=="homeless" to de- fine unambiguously which proportion we are considering. We could also have written homeless=="housed".

can also work directly with numerical counts, the way binom.test() does.

prop.test() calculates a Chi- squared statistic. Most intro- ductory texts use a z-statistic. They are mathematically equiv- alent in terms of inferential statements, but you may need to address the discrepancy with your students.

> prop.test(209, 209 + 244, correct=FALSE)

1-sample proportions test without continuity correction

data: 209 out of 209 + 244

X-squared = 3, df = 1, p-value = 0.1

alternative hypothesis: true p is not equal to 0.5

95 percent confidence interval:

0.416 0.507

sample estimates:

p

0.461

4.4 Goodness of fit tests

A variety of goodness of fit tests can be calculated against a reference distribution. For the HELP data, we could test the null hypothesis that there is an equal proportion of subjects in each substance abuse group back in the original populations.

> tally(~ substance, format="percent", data=HELPrct)

alcohol cocaine heroin

39.1 33.6 27.4

a student’s guide to r 43

> observed <- tally(~ substance, data=HELPrct)

> observed

alcohol cocaine heroin

177 152 124 Caution! In addition to the format op- tion, there is an option margins to include marginal totals in the table. The default in tally() is margins=FALSE. Try it out!> p <- c(1/3, 1/3, 1/3) # equivalent to rep(1/3, 3)

> chisq.test(observed, p=p)

Chi-squared test for given probabilities

data: observed

X-squared = 9, df = 2, p-value = 0.01

> total <- sum(observed); total

[1] 453

> expected <- total*p; expected

[1] 151 151 151

We can also calculate the χ2 statistic manually, as a function of observed and expected values.

> chisq <- sum((observed - expected)^2/(expected)); chisq

[1] 9.31

> 1 - pchisq(chisq, df=2)

[1] 0.00951 The pchisq() function calcu- lates the probability that a χ2

random variable with df() de- grees is freedom is less than or equal to a given value. Here we calculate the complement to find the area to the right of the observed Chi-square statistic.

It may be helpful to consult a graph of the statistic, where the shaded area represents the value to the right of the observed value.

> plotDist("chisq", df=2, groups = x > 9.31, type="h")

44 horton, kaplan, pruim

0.0

0.1

0.2

0.3

0.4

0.5

0 5 10 15

Alternatively, the mosaic package provides a version of chisq.test() with more verbose output.

> xchisq.test(observed, p=p)

Chi-squared test for given probabilities

data: x

X-squared = 9, df = 2, p-value = 0.01

177 152 124

(151.00) (151.00) (151.00)

[4.4768] [0.0066] [4.8278]

< 2.116> < 0.081> <-2.197>

key:

observed

(expected)

[contribution to X-squared]

<residual> x in xchisq.test() stands for eXtra.

Objects in the workspace are listed in the Environment tab in RStudio. If you want to clean up that listing, remove objects that are no longer needed with rm().

> # clean up variables no longer needed

> rm(observed, p, total, chisq)

5 Two Quantitative Variables

5.1 Scatterplots

We always encourage students to start any analysis by graphing their data. Here we augment a scatterplot of the CESD (a measure of depressive symptoms, higher scores indicate more symptoms) and the MCS (mental component score from the SF-36, where higher scores indicate better functioning) for female subjects with a lowess (locally weighted scatterplot smoother) line, using a circle as the plotting character and slightly thicker line.

The lowess line can help to as- sess linearity of a relationship. This is added by specifying both points (using ‘p’) and a lowess smoother.

> females <- filter(HELPrct, female==1)

> xyplot(cesd ~ mcs, type=c("p","smooth"), pch=1,

cex=0.6, lwd=3, data=females)

mcs

ce sd

10

20

30

40

50

60

10 20 30 40 50 60

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

● ●

●

●

●

●

● ●

●

●

●

●

●

●

● ●

●

●

●●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

● ●●

●

●

●

●

●

●

●

●

●● ●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

46 horton, kaplan, pruim

Digging Deeper

The Start Modeling with R com- panion book will be helpful if you are unfamiliar with the modeling language. The Start Teaching with R also provides useful guidance in getting started.

It’s straightforward to plot something besides a char- acter in a scatterplot. In this example, the USArrests can be used to plot the association between murder and as- sault rates, with the state name displayed. This requires a panel function to be written.

> panel.labels <- function(x, y, labels='x',...) {

panel.text(x, y, labels, cex=0.4, ...)

}

> xyplot(Murder ~ Assault, panel=panel.labels,

labels=rownames(USArrests), data=USArrests)

Assault

M ur

de r

5

10

15

50 100 150 200 250 300 350

Alabama

Alaska

Arizona Arkansas California

Colorado

Connecticut

Delaware

Florida

Georgia

Hawaii

Idaho

Illinois

Indiana

Iowa

Kansas

Kentucky

Louisiana

Maine

Maryland

Massachusetts

Michigan

Minnesota

Mississippi

Missouri

Montana

Nebraska

Nevada

New Hampshire

New Jersey

New MexicoNew York

North Carolina

North Dakota

Ohio Oklahoma

Oregon

Pennsylvania

Rhode Island

South Carolina

South Dakota

Tennessee Texas

Utah Vermont

Virginia

Washington

West Virginia

Wisconsin

Wyoming

5.2 Correlation

Correlations can be calculated for a pair of variables, or for a matrix of variables.

> cor(cesd, mcs, data=females)

[1] -0.674

> smallHELP <- select(females, cesd, mcs, pcs)

> cor(smallHELP)

cesd mcs pcs

cesd 1.000 -0.674 -0.369

mcs -0.674 1.000 0.266

pcs -0.369 0.266 1.000

By default, Pearson correlations are provided. Other variants (e.g., Spearman) can be specified using the method

a student’s guide to r 47

option.

> cor(cesd, mcs, method="spearman", data=females)

[1] -0.666

5.3 Pairs plots

A pairs plot (scatterplot matrix) can be calculated for each pair of a set of variables.

The GGally package has sup- port for more elaborate pairs plots.

> splom(smallHELP, cex=0.3)

Scatter Plot Matrix

cesd 40 50 60

40 50 60

10 20 30

10 20 30 ●

●

● ●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

● ●

●

●

●

●

● ●

●

●

●

●

●

●

● ●

●

●

●● ●

●

●

●

●

●

●

●

●

● ●

●

●

● ●

●

●

●

●

●

●

●

● ●●

●

●

●

●

●

● ●

●

●● ●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

● ●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

●●

●

●

●

●

● ●

●

●

●

●

●

●

●●

●

●

● ● ●

●

●

●

●

●

●

●

●

●●

●

●

● ●

●

●

●

●

●

●

●

●● ●

●

●

●

●

●

● ●

●

●●●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

● ●

●

● ●

●

●

●

●

● ●

●

●

●

●

●

●

● ●

●

●

● ●

●

●

●

●

●

●

●

●

●

●●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

● mcs 40 50 60

40 50 60

10 20 30

10 20 30

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ● ●

●

● ●

●

● ●

●

●

●

●

● ●

●

●

●

●

●

●

● ●

●

●

● ●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

● ●

● ●

●

●

●

● ● ●

● ●

●

●

●

●

●

● ●

●

● ●

● ●

● ●

● ●

●●

●

●

●

●

●

●

●

●

●

●●

●

●

●

●

● ●

●

●

●

● ●

●

●

●

●

●

●

●

●

●● ●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

● ●

●

●

● ●

●

●

●

●

●●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●●

●●

●

●

●

● ● ●

● ●

●

●

●

●

●

●●

●

● ●

● ●

● ●

● ●

● ●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

● ●

●

●

●

● ●

●

●

●

●

●

●

●

●

●● ●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

● ●

●

●

● ●

●

●

●

●

● ●

●

●

●

●

●

pcs 50

60

70

50 60 70

20

30

40

20 30 40

5.4 Simple linear regression

We tend to introduce linear re- gression early in our courses, as a purely descriptive technique.

Linear regression models are described in detail in Start Modeling with R. These use the same formula interface introduced previously for numerical and graphical sum- maries to specify the outcome and predictors. Here we consider fitting the model cesd ∼ mcs.

> cesdmodel <- lm(cesd ~ mcs, data=females)

> coef(cesdmodel)

48 horton, kaplan, pruim

(Intercept) mcs

57.349 -0.707

It’s important to pick good names for modeling objects. Here the output of lm() is saved as cesdmodel, which denotes that it is a regression model of depressive symptom scores.

To simplify the output, we turn off the option to dis- play significance stars.

> options(show.signif.stars=FALSE)

> coef(cesdmodel)

(Intercept) mcs

57.349 -0.707

> msummary(cesdmodel)

Estimate Std. Error t value Pr(>|t|)

(Intercept) 57.3485 2.3806 24.09 < 2e-16

mcs -0.7070 0.0757 -9.34 1.8e-15

Residual standard error: 9.66 on 105 degrees of freedom

Multiple R-squared: 0.454,Adjusted R-squared: 0.449

F-statistic: 87.3 on 1 and 105 DF, p-value: 1.81e-15

> coef(summary(cesdmodel))

Estimate Std. Error t value Pr(>|t|)

(Intercept) 57.349 2.3806 24.09 1.42e-44

mcs -0.707 0.0757 -9.34 1.81e-15

> confint(cesdmodel)

2.5 % 97.5 %

(Intercept) 52.628 62.069

mcs -0.857 -0.557

> rsquared(cesdmodel)

[1] 0.454

> class(cesdmodel)

[1] "lm"

The return value from lm() is a linear model object. A number of functions can operate on these objects, as seen previously with coef(). The function residuals() re- turns a vector of the residuals.

The function residuals() can be abbreviated resid(). Another useful function is fitted(), which returns a vec- tor of predicted values.

a student’s guide to r 49

> histogram(~ residuals(cesdmodel), density=TRUE)

residuals(cesdmodel)

D en

si ty

0.00

0.01

0.02

0.03

0.04

−30 −20 −10 0 10 20

> qqmath(~ resid(cesdmodel))

qnorm

re si

d( ce

sd m

od el

)

−20

−10

0

10

20

−2 −1 0 1 2

● ● ●●

●●●●●● ●●●●

●●●●●● ●●●●●

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

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

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

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

●●●●●●● ●●●●

●●●● ●●●●●●●

● ●

> xyplot(resid(cesdmodel) ~ fitted(cesdmodel), type=c("p", "smooth", "r"),

alpha=0.5, cex=0.3, pch=20)

50 horton, kaplan, pruim

fitted(cesdmodel)

re si

d( ce

sd m

od el

)

−20

−10

0

10

20

20 30 40 50

The mplot() function can facilitate creating a variety of useful plots, including the same residuals vs. fitted scatterplots, by specifying the which=1 option.

> mplot(cesdmodel, which=1)

[[1]]

Residuals vs Fitted

Fitted Value

R es

id ua

l

−20

−10

0

10

20

20 30 40 50

●●

● ●

●

● ●

●

●

●

●

● ●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

●●

● ●

● ●

●

●

● ●

●

●

●

●

●●

●

●

●

●

●●

●

● ●

●

● ● ●

● ●

● ●

●●

●●

●

● ● ●

●

●

●

● ●

●●

●

● ●●

●

● ●

●

●

●

● ●

●●

●

●

●

●●

●

●

●

●

●

It can also generate a normal quantile-quantile plot (which=2),

> mplot(cesdmodel, which=2)

a student’s guide to r 51

Normal Q−Q

Theoretical Quantiles

S ta

nd ar

di ze

d re

si du

al s

−2

−1

0

1

2

−2 −1 0 1 2

● ● ●●

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

●●●●●● ●●●●●

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

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

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

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

●●●● ●●●●

●●●●●●● ●

●

scale vs. location,

> mplot(cesdmodel, which=3)

Scale−Location

Fitted Value

S ta

nd ar

di ze

d re

si du

al s

0.5

1.0

1.5

20 30 40 50

●● ● ●

●

●●

●

●

●

● ● ●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

● ●

●

● ●

●

●

● ●

●

●

●

●

●

●

●

●

●●

●●

●

●

●●

●

●

●

●

●● ●

●

● ●

●

●●

●●

●

●

● ●

●

●

● ●

●

● ●

●

●

●●

●

● ●

●

●

●

●

●

●● ●

●

●

● ●

●

●

●

●

●

Cook’s distance by observation number,

> mplot(cesdmodel, which=4)

52 horton, kaplan, pruim

Cook's Distance

Observation number

C oo

k' s

di st

an ce

0.00

0.05

0.10

0.15

0 20 40 60 80 100

●

●

●●

●

●●●

● ●

●●● ● ●●●●●

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

● ● ●●●●●

●● ● ●●●●●

● ●●

●

●●● ● ●●●

●

●● ●●●●●●●

●●●● ● ●●●

● ●●

●●●● ● ●●●●

●

● ●●●

● ●●

●

●

●

●

●●

●

●●●●

residuals vs. leverage,

> mplot(cesdmodel, which=5)

Residuals vs Leverage

Leverage

S ta

nd ar

di ze

d R

es id

ua ls

−2

−1

0

1

2

0.01 0.02 0.03 0.04 0.05 0.06 0.07

● ●

●●

●

● ●

●

●

●

●

●●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●●

●●

● ●

● ●

●

●

● ●

●

●

●

●

● ●

●

●

●

●

●●

●

● ● ●

● ●●

● ●

● ●

●●

●●

●

● ● ●

●

●

●

● ●

●●

●

● ●●

●

●●

●

●

●

● ●

● ●

●

●

●

● ●

●

●

●

●

●

and Cook’s distance vs. leverage.

> mplot(cesdmodel, which=6)

a student’s guide to r 53

Cook's dist vs Leverage

Leverage

C oo

k' s

di st

an ce

0.00

0.05

0.10

0.15

0.01 0.02 0.03 0.04 0.05 0.06 0.07

●

●

●●

●

●● ●

● ●

●●● ●

● ●● ●● ●

● ●

●● ●●

● ●

● ●●● ●

● ●●● ●●

● ● ●

● ●●● ●● ● ●

●

●● ● ●

●● ●

●

●● ● ●● ● ●●●

●● ●● ●

●● ● ●

● ●

●● ●● ●

●●●●

●

● ● ●●

● ●●

●

●

●

●

● ●

●

●● ●●

Prediction bands can be added to a plot using the panel.lmbands() function.

> xyplot(cesd ~ mcs, panel=panel.lmbands, cex=0.2, band.lwd=2, data=HELPrct)

mcs

ce sd

0

10

20

30

40

50

60

10 20 30 40 50 60

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

● ●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

● ●

● ●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●●

●

●

● ●

●

●

●

●

●

● ●

●

●

●

● ●

●

●

●

●

●

●

●

●

● ●

●

● ●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

● ●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

● ●

● ●

● ●●

●

●

● ●

●

●

●

●

●

● ●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

6 Two Categorical Variables

6.1 Cross classification tables

Cross classification (two-way or R by C) tables can be constructed for two (or more) categorical variables. Here we consider the contingency table for homeless status (homeless one or more nights in the past 6 months or housed) and sex.

> tally(~ homeless + sex, margins=FALSE, data=HELPrct)

sex

homeless female male

homeless 40 169

housed 67 177

We can also calculate column percentages:

> tally(~ sex | homeless, margins=TRUE, format="percent",

data=HELPrct)

homeless

sex homeless housed

female 19.1 27.5

male 80.9 72.5

Total 100.0 100.0

We can calculate the odds ratio directly from the table:

56 horton, kaplan, pruim

> OR <- (40/169)/(67/177); OR

[1] 0.625

The mosaic package has a function which will calcu- late odds ratios:

> oddsRatio(tally(~ (homeless=="housed") + sex, margins=FALSE,

data=HELPrct))

[1] 0.625

The CrossTable() function in the gmodels package also displays a cross classification table.

> require(gmodels)

> with(HELPrct, CrossTable(homeless, sex,

prop.r=FALSE, prop.chisq=FALSE, prop.t=FALSE))

Cell Contents

|-------------------------|

| N |

| N / Col Total |

|-------------------------|

Total Observations in Table: 453

| sex

homeless | female | male | Row Total |

-------------|-----------|-----------|-----------|

homeless | 40 | 169 | 209 |

| 0.374 | 0.488 | |

-------------|-----------|-----------|-----------|

housed | 67 | 177 | 244 |

| 0.626 | 0.512 | |

-------------|-----------|-----------|-----------|

Column Total | 107 | 346 | 453 |

| 0.236 | 0.764 | |

-------------|-----------|-----------|-----------|

a student’s guide to r 57

Graphical summaries of cross classification tables may be helpful in visualizing associations. Mosaic plots are one example, where the total area (all observations) is proportional to one. Here we see that males tend to be Caution!

The jury is still out regarding the utility of mosaic plots (also known as eikosograms), relative to the low data to ink ratio. We have found them to be helpful to reinforce understanding of a two way contingency table.

E. R. Tufte. The Visual Dis- play of Quantitative Information. Graphics Press, Cheshire, CT, 2nd edition, 2001

over-represented amongst the homeless subjects (as rep- resented by the horizontal line which is higher for the homeless rather than the housed).

The mosaic() function in the vcd package also makes mosaic plots.

> mytab <- tally(~ homeless + sex, margins=FALSE,

data=HELPrct)

> mosaicplot(mytab)

mytab

homeless

se x

homeless housed

fe m

al e

m al

e

6.2 Creating tables from summary statistics

Tables can be created from summary statistics using the do() function.

> HELPtable <- rbind(

do(40) * data.frame(sex="female", homeless="homeless"),

do(169) * data.frame(sex="male", homeless="homeless"),

do(67) * data.frame(sex="female", homeless="housed"),

58 horton, kaplan, pruim

do(177) * data.frame(sex="male", homeless="housed")

)

> tally(~ homeless + sex, data=HELPtable)

sex

homeless female male

homeless 40 169

housed 67 177

6.3 Chi-squared tests

> chisq.test(tally(~ homeless + sex, margins=FALSE,

data=HELPrct), correct=FALSE)

Pearson's Chi-squared test

data: tally(~homeless + sex, margins = FALSE, data = HELPrct)

X-squared = 4, df = 1, p-value = 0.04

There is a statistically significant association found: it is unlikely that we would observe an association this strong if homeless status and sex were independent in the population.

When a student finds a significant association, it’s im- portant for them to be able to interpret this in the context of the problem. The xchisq.test() function provides additional details (observed, expected, contribution to statistic, and residual) to help with this process.

x is for eXtra. > xchisq.test(tally(~homeless + sex, margins=FALSE,

data=HELPrct), correct=FALSE)

Pearson's Chi-squared test

data: x

X-squared = 4, df = 1, p-value = 0.04

40 169

a student’s guide to r 59

( 49.37) (159.63)

[1.78] [0.55]

<-1.33> < 0.74>

67 177

( 57.63) (186.37)

[1.52] [0.47]

< 1.23> <-0.69>

key:

observed

(expected)

[contribution to X-squared]

<residual>

We observe that there are fewer homeless women, and more homeless men that would be expected.

6.4 Fisher’s exact test

An exact test can also be calculated. This is computation- ally straightforward for 2 by 2 tables. Options to help constrain the size of the problem for larger tables exist (see ?fisher.test()). Digging Deeper

Note the different estimate of the odds ratio from that seen in section 6.1. The fisher.test() function uses a different estima- tor (and different interval based on the profile likelihood).

> fisher.test(tally(~homeless + sex, margins=FALSE,

data=HELPrct))

Fisher's Exact Test for Count Data

data: tally(~homeless + sex, margins = FALSE, data = HELPrct)

p-value = 0.05

alternative hypothesis: true odds ratio is not equal to 1

95 percent confidence interval:

0.389 0.997

sample estimates:

odds ratio

0.626

7 Quantitative Response, Categorical Predictor

7.1 A dichotomous predictor: numerical and graphical summaries

Here we will compare the distributions of CESD scores by sex. The mean() function can be used to calculate the mean CESD score separately for males and females.

> mean(cesd ~ sex, data=HELPrct)

female male

36.9 31.6

The favstats() function can provide more statistics by group.

> favstats(cesd ~ sex, data=HELPrct)

sex min Q1 median Q3 max mean sd n missing

1 female 3 29 38.0 46.5 60 36.9 13.0 107 0

2 male 1 24 32.5 40.0 58 31.6 12.1 346 0

Boxplots are a particularly helpful graphical display to compare distributions. The bwplot() function can be used to display the boxplots for the CESD scores sepa- rately by sex. We see from both the numerical and graph- ical summaries that women tend to have slightly higher CESD scores than men.

Although we usually put ex- planatory variables along the horizontal axis, page layout sometimes makes the other orientation preferable for these plots.

> bwplot(sex ~ cesd, data=HELPrct)

62 horton, kaplan, pruim

cesd

female

male

0 10 20 30 40 50 60

●

●

When sample sizes are small, there is no reason to summarize with a boxplot since xyplot() can handle categorical predictors. Even with 10–20 observations in a group, a scatter plot is often quite readable. Setting the alpha level helps detect multiple observations with the same value.

One of us once saw a biologist proudly present side-by-side boxplots. Thinking a major vic- tory had been won, he naively asked how many observations were in each group. “Four,” replied the biologist.

> xyplot(sex ~ length, alpha=.6, cex=1.4, data=KidsFeet)

length

se x

B

G

22 23 24 25 26 27

7.2 A dichotomous predictor: two-sample t

The Student’s two sample t-test can be run without (de- fault) or with an equal variance assumption.

> t.test(cesd ~ sex, var.equal=FALSE, data=HELPrct)

Welch Two Sample t-test

data: cesd by sex

t = 4, df = 200, p-value = 3e-04

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

2.49 8.09

a student’s guide to r 63

sample estimates:

mean in group female mean in group male

36.9 31.6

We see that there is a statistically significant difference between the two groups.

We can repeat using the equal variance assumption.

> t.test(cesd ~ sex, var.equal=TRUE, data=HELPrct)

Two Sample t-test

data: cesd by sex

t = 4, df = 500, p-value = 1e-04

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

2.61 7.97

sample estimates:

mean in group female mean in group male

36.9 31.6

The groups can also be compared using the lm() func- tion (also with an equal variance assumption). The mo- saic command msummary() provides a slightly terser ver- sion of the typical output from summary().

> msummary(lm(cesd ~ sex, data=HELPrct))

Estimate Std. Error t value Pr(>|t|)

(Intercept) 36.89 1.19 30.96 < 2e-16

sexmale -5.29 1.36 -3.88 0.00012

Residual standard error: 12.3 on 451 degrees of freedom

Multiple R-squared: 0.0323,Adjusted R-squared: 0.0302

F-statistic: 15.1 on 1 and 451 DF, p-value: 0.00012

The lm() function is part of a much more flexible modeling framework while t.test() is essentially a dead end. lm() uses of the equal variance as- sumption. See the companion book, Start Modeling in R for more details.

7.3 Non-parametric 2 group tests

The same conclusion is reached using a non-parametric (Wilcoxon rank sum) test.

64 horton, kaplan, pruim

> wilcox.test(cesd ~ sex, data=HELPrct)

Wilcoxon rank sum test with continuity correction

data: cesd by sex

W = 20000, p-value = 1e-04

alternative hypothesis: true location shift is not equal to 0

7.4 Permutation test

Here we extend the methods introduced in section 3.6 to undertake a two-sided test comparing the ages at baseline by gender. First we calculate the observed difference in means:

> mean(age ~ sex, data=HELPrct)

female male

36.3 35.5

> test.stat <- diffmean(age ~ sex, data=HELPrct)

> test.stat

diffmean

-0.784

We can calculate the same statistic after shuffling the group labels:

> do(1) * diffmean(age ~ shuffle(sex), data=HELPrct)

diffmean

1 0.268

> do(1) * diffmean(age ~ shuffle(sex), data=HELPrct)

diffmean

1 0.0357

> do(3) * diffmean(age ~ shuffle(sex), data=HELPrct)

diffmean

1 -0.515

2 -0.613

3 0.305

a student’s guide to r 65

Digging Deeper

More details and examples can be found in the mosaic package Resampling Vignette.

> rtest.stats <- do(500) * diffmean(age ~ shuffle(sex),

data=HELPrct)

> head(rtest.stats, 3)

diffmean

1 -0.1479

2 0.3171

3 0.0479

> favstats(~ diffmean, data=rtest.stats)

min Q1 median Q3 max mean sd n missing

-2.2 -0.527 0.0602 0.647 2.3 0.0502 0.871 500 0

> histogram(~ diffmean, n=40, xlim=c(-6, 6),

groups=diffmean >= test.stat, pch=16, cex=.8,

data=rtest.stats)

> ladd(panel.abline(v=test.stat, lwd=3, col="red"))

diffmean

D en

si ty

0.0

0.2

0.4

0.6

−4 −2 0 2 4

Here we don’t see much evidence to contradict the null hypothesis that men and women have the same mean age in the population.

7.5 One-way ANOVA

Earlier comparisons were between two groups. We can also consider testing differences between three or more groups using one-way ANOVA. Here we compare CESD scores by primary substance of abuse (heroin, cocaine, or

66 horton, kaplan, pruim

alcohol) with a line rather a dot to indicate the median.

> bwplot(cesd ~ substance, pch="|", data=HELPrct)

ce sd

0

10

20

30

40

50

60

alcohol cocaine heroin

●

> mean(cesd ~ substance, data=HELPrct)

alcohol cocaine heroin

34.4 29.4 34.9

> anovamod <- aov(cesd ~ substance, data=HELPrct)

> summary(anovamod)

Df Sum Sq Mean Sq F value Pr(>F)

substance 2 2704 1352 8.94 0.00016

Residuals 450 68084 151

While still high (scores of 16 or more are generally con- sidered to be “severe” symptoms), the cocaine-involved group tend to have lower scores than those whose pri- mary substances are alcohol or heroin.

> modintercept <- lm(cesd ~ 1, data=HELPrct)

> modsubstance <- lm(cesd ~ substance, data=HELPrct)

The anova() command can summarize models.

> anova(modsubstance)

Analysis of Variance Table

Response: cesd

Df Sum Sq Mean Sq F value Pr(>F)

substance 2 2704 1352 8.94 0.00016

Residuals 450 68084 151

a student’s guide to r 67

In this setting the results are identical (since there is only one predictor, with 2 df).

The anova() function can also be used to formally compare two (nested) models.

> anova(modintercept, modsubstance)

Analysis of Variance Table

Model 1: cesd ~ 1

Model 2: cesd ~ substance

Res.Df RSS Df Sum of Sq F Pr(>F)

1 452 70788

2 450 68084 2 2704 8.94 0.00016

7.6 Tukey’s Honest Significant Differences

There are a variety of multiple comparison procedures that can be used after fitting an ANOVA model. One of these is Tukey’s Honest Significant Differences (HSD). Other options are available within the multcomp package.

> favstats(cesd ~ substance, data=HELPrct)

substance min Q1 median Q3 max mean sd n missing

1 alcohol 4 26 36 42 58 34.4 12.1 177 0

2 cocaine 1 19 30 39 60 29.4 13.4 152 0

3 heroin 4 28 35 43 56 34.9 11.2 124 0

> HELPrct <- mutate(HELPrct, subgrp = factor(substance,

levels=c("alcohol", "cocaine", "heroin"),

labels=c("A", "C", "H")))

> mod <- lm(cesd ~ subgrp, data=HELPrct)

> HELPHSD <- TukeyHSD(mod, "subgrp")

> HELPHSD

Tukey multiple comparisons of means

95% family-wise confidence level

Fit: aov(formula = x)

68 horton, kaplan, pruim

$subgrp

diff lwr upr p adj

C-A -4.952 -8.15 -1.75 0.001

H-A 0.498 -2.89 3.89 0.936

H-C 5.450 1.95 8.95 0.001

> mplot(HELPHSD)

Tukey's Honest Significant Differences

difference in means

H−C

H−A

C−A

−5 0 5

●

●

●

subgrp

Again, we see that the cocaine group has significantly lower CESD scores than either of the other two groups.

8 Categorical Response, Quantitative Predictor

8.1 Logistic regression

Logistic regression is available using the glm() function, which supports a variety of link functions and distribu- tional forms for generalized linear models, including lo- gistic regression.

The glm() function has argu- ment family, which can take an option link. The logit link is the default link for the binomial family, so we don’t need to specify it here. The more verbose usage would be family=binomial(link=logit).

> logitmod <- glm(homeless ~ age + female, family=binomial,

data=HELPrct)

> msummary(logitmod)

Coefficients:

Estimate Std. Error z value Pr(>|z|)

(Intercept) 0.8926 0.4537 1.97 0.049

age -0.0239 0.0124 -1.92 0.055

female 0.4920 0.2282 2.16 0.031

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 625.28 on 452 degrees of freedom

Residual deviance: 617.19 on 450 degrees of freedom

AIC: 623.2

Number of Fisher Scoring iterations: 4

> exp(coef(logitmod))

(Intercept) age female

2.442 0.976 1.636

> exp(confint(logitmod))

Waiting for profiling to be done...

70 horton, kaplan, pruim

2.5 % 97.5 %

(Intercept) 1.008 5.99

age 0.953 1.00

female 1.050 2.57

We can compare two models (for multiple degree of freedom tests). For example, we might be interested in the association of homeless status and age for each of the three substance groups.

> mymodsubage <- glm((homeless=="homeless") ~ age + substance,

family=binomial, data=HELPrct)

> mymodage <- glm((homeless=="homeless") ~ age, family=binomial,

data=HELPrct)

> msummary(mymodsubage)

Coefficients:

Estimate Std. Error z value Pr(>|z|)

(Intercept) -0.0509 0.5164 -0.10 0.9215

age 0.0100 0.0129 0.77 0.4399

substancecocaine -0.7496 0.2303 -3.25 0.0011

substanceheroin -0.7780 0.2469 -3.15 0.0016

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 625.28 on 452 degrees of freedom

Residual deviance: 607.62 on 449 degrees of freedom

AIC: 615.6

Number of Fisher Scoring iterations: 4

> exp(coef(mymodsubage))

(Intercept) age substancecocaine substanceheroin

0.950 1.010 0.473 0.459

> anova(mymodage, mymodsubage, test="Chisq")

Analysis of Deviance Table

Model 1: (homeless == "homeless") ~ age

Model 2: (homeless == "homeless") ~ age + substance

Resid. Df Resid. Dev Df Deviance Pr(>Chi)

1 451 622

2 449 608 2 14.3 0.00078

a student’s guide to r 71

We observe that the cocaine and heroin groups are sig- nificantly less likely to be homeless than alcohol involved subjects, after controlling for age. (A similar result is seen when considering just homeless status and substance.)

> tally(~ homeless | substance, format="percent", margins=TRUE, data=HELPrct)

substance

homeless alcohol cocaine heroin

homeless 58.2 38.8 37.9

housed 41.8 61.2 62.1

Total 100.0 100.0 100.0

9 Survival Time Outcomes

Extensive support for survival (time to event) analysis is available within the survival package.

9.1 Kaplan-Meier plot

> require(survival)

> fit <- survfit(Surv(dayslink, linkstatus) ~ treat,

data=HELPrct)

> plot(fit, conf.int=FALSE, lty=1:2, lwd=2,

xlab="time (in days)", ylab="P(not linked)")

> legend(20, 0.4, legend=c("Control", "Treatment"),

lty=c(1,2), lwd=2)

> title("Product-Limit Survival Estimates (time to linkage)")

0 100 200 300 400

0. 0

0. 2

0. 4

0. 6

0. 8

1. 0

time (in days)

P (n

ot li

nk ed

)

Control Treatment

Product−Limit Survival Estimates (time to linkage)

74 horton, kaplan, pruim

We see that the subjects in the treatment (Health Eval- uation and Linkage to Primary Care clinic) were signif- icantly more likely to link to primary care (less likely to “survive”) than the control (usual care) group.

9.2 Cox proportional hazards model

> require(survival)

> summary(coxph(Surv(dayslink, linkstatus) ~ age + substance,

data=HELPrct))

Call:

coxph(formula = Surv(dayslink, linkstatus) ~ age + substance,

data = HELPrct)

n= 431, number of events= 163

(22 observations deleted due to missingness)

coef exp(coef) se(coef) z Pr(>|z|)

age 0.00893 1.00897 0.01026 0.87 0.38

substancecocaine 0.18045 1.19775 0.18100 1.00 0.32

substanceheroin -0.28970 0.74849 0.21725 -1.33 0.18

exp(coef) exp(-coef) lower .95 upper .95

age 1.009 0.991 0.989 1.03

substancecocaine 1.198 0.835 0.840 1.71

substanceheroin 0.748 1.336 0.489 1.15

Concordance= 0.55 (se = 0.023 )

Rsquare= 0.014 (max possible= 0.988 )

Likelihood ratio test= 6.11 on 3 df, p=0.106

Wald test = 5.84 on 3 df, p=0.12

Score (logrank) test = 5.91 on 3 df, p=0.116

Neither age nor substance group was significantly associated with linkage to primary care.

10 More than Two Variables

10.1 Two (or more) way ANOVA

We can fit a two (or more) way ANOVA model, without or with an interaction, using the same modeling syntax.

> median(cesd ~ substance | sex, data=HELPrct)

alcohol.female cocaine.female heroin.female alcohol.male

40.0 35.0 39.0 33.0

cocaine.male heroin.male female male

29.0 34.5 38.0 32.5

> bwplot(cesd ~ subgrp | sex, data=HELPrct)

ce sd

0

10

20

30

40

50

60

A C H

● ●

●

female

A C H

● ●

●

●

male

> summary(aov(cesd ~ substance + sex, data=HELPrct))

Df Sum Sq Mean Sq F value Pr(>F)

substance 2 2704 1352 9.27 0.00011

76 horton, kaplan, pruim

sex 1 2569 2569 17.61 3.3e-05

Residuals 449 65515 146

> summary(aov(cesd ~ substance * sex, data=HELPrct))

Df Sum Sq Mean Sq F value Pr(>F)

substance 2 2704 1352 9.25 0.00012

sex 1 2569 2569 17.57 3.3e-05

substance:sex 2 146 73 0.50 0.60752

Residuals 447 65369 146

There’s little evidence for the interaction, though there are statistically significant main effects terms for substance group and sex.

> xyplot(cesd ~ substance, groups=sex, lwd=2,

auto.key=list(columns=2, lines=TRUE, points=FALSE),

type='a', data=HELPrct)

substance

ce sd

0

10

20

30

40

50

60

alcohol cocaine heroin

female male

10.2 Multiple regression

Multiple regression is a logical extension of the prior commands, where additional predictors are added. This allows students to start to try to disentangle multivariate relationships.

We tend to introduce multiple linear regression early in our courses, as a purely descriptive technique, then return to it regularly. The motivation for this is described at length in the companion volume Start Modeling with R.

Here we consider a model (parallel slopes) for depres- sive symptoms as a function of Mental Component Score (MCS), age (in years) and sex of the subject.

a student’s guide to r 77

> lmnointeract <- lm(cesd ~ mcs + age + sex, data=HELPrct)

> msummary(lmnointeract)

Estimate Std. Error t value Pr(>|t|)

(Intercept) 53.8303 2.3617 22.79 <2e-16

mcs -0.6548 0.0336 -19.50 <2e-16

age 0.0553 0.0556 1.00 0.3200

sexmale -2.8993 1.0137 -2.86 0.0044

Residual standard error: 9.09 on 449 degrees of freedom

Multiple R-squared: 0.476,Adjusted R-squared: 0.473

F-statistic: 136 on 3 and 449 DF, p-value: <2e-16

We can also fit a model that includes an interaction be- tween MCS and sex.

> lminteract <- lm(cesd ~ mcs + age + sex + mcs:sex, data=HELPrct)

> msummary(lminteract)

Estimate Std. Error t value Pr(>|t|)

(Intercept) 55.3906 2.9903 18.52 <2e-16

mcs -0.7082 0.0712 -9.95 <2e-16

age 0.0549 0.0556 0.99 0.324

sexmale -4.9421 2.6055 -1.90 0.058

mcs:sexmale 0.0687 0.0807 0.85 0.395

Residual standard error: 9.09 on 448 degrees of freedom

Multiple R-squared: 0.477,Adjusted R-squared: 0.472

F-statistic: 102 on 4 and 448 DF, p-value: <2e-16

> anova(lminteract)

Analysis of Variance Table

Response: cesd

Df Sum Sq Mean Sq F value Pr(>F)

mcs 1 32918 32918 398.27 <2e-16

age 1 107 107 1.29 0.2563

sex 1 676 676 8.18 0.0044

mcs:sex 1 60 60 0.72 0.3952

Residuals 448 37028 83

> anova(lmnointeract, lminteract)

Analysis of Variance Table

78 horton, kaplan, pruim

Model 1: cesd ~ mcs + age + sex

Model 2: cesd ~ mcs + age + sex + mcs:sex

Res.Df RSS Df Sum of Sq F Pr(>F)

1 449 37088

2 448 37028 1 59.9 0.72 0.4

There is little evidence for an interaction effect, so we drop this from the model.

10.2.1 Visualizing the results from the regression

The makeFun() and plotFun() functions from the mosaic package can be used to display the predicted values from a regression model. For this example, we might display the predicted CESD values for a range of MCS (mental component score) values a hypothetical 36 year old male and female subject might have from the parallel slopes (no interaction) model.

> lmfunction <- makeFun(lmnointeract)

We can now plot the predicted values separately for male and female subjects over a range of MCS (mental component score) values, along with the observed data for all of the 36 year olds.

> xyplot(cesd ~ mcs, groups=sex, auto.key=TRUE,

data=filter(HELPrct, age==36))

> plotFun(lmfunction(mcs, age=36, sex="male") ~ mcs,

xlim=c(0, 60), lwd=2, ylab="predicted CESD", add=TRUE)

> plotFun(lmfunction(mcs, age=36, sex="female") ~ mcs,

xlim=c(0, 60), lty=2, lwd=3, add=TRUE)

a student’s guide to r 79

mcs

ce sd

10

20

30

40

50

20 30 40 50 60

●

●

●

female male

●

10.2.2 Coefficient plots

It is sometimes useful to display a plot of the coefficients for a multiple regression model (along with their associ- ated confidence intervals).

> mplot(lmnointeract, rows=-1, which=7)

95% confidence intervals

estimate

co ef

fic ie

nt

sexmale

age

mcs

−5 −4 −3 −2 −1 0

Darker dots indicate regression coefficients where the 95% confidence interval does not include the null hypothesis value of zero.

Caution! Be careful when fitting re- gression models with missing values (see also section 13.11).10.2.3 Residual diagnostics

It’s straightforward to undertake residual diagnostics for this model. We begin by adding the fitted values and residuals to the dataset.

The mplot() function can also be used to create these graphs.

Here we are adding two new variables into an existing dataset. It’s often a good practice to give the resulting dataframe a new name.

80 horton, kaplan, pruim

> HELPrct <- mutate(HELPrct,

residuals = resid(lmnointeract),

pred = fitted(lmnointeract))

> histogram(~ residuals, xlab="residuals", fit="normal",

data=HELPrct)

residuals

D en

si ty

0.00

0.01

0.02

0.03

0.04

−30 −20 −10 0 10 20

We can identify the subset of observations with ex- tremely large residuals.

> filter(HELPrct, abs(residuals) > 25)

age anysubstatus anysub cesd d1 daysanysub dayslink drugrisk e2b

1 43 0 no 16 15 191 414 0 NA

2 27 NA <NA> 40 1 NA 365 3 2

female sex g1b homeless i1 i2 id indtot linkstatus link mcs pcs

1 0 male no homeless 24 36 44 41 0 no 15.9 71.4

2 0 male no homeless 18 18 420 37 0 no 57.5 37.7

pss_fr racegrp satreat sexrisk substance treat subgrp residuals

1 3 white no 7 cocaine yes C -26.9

2 8 white yes 3 heroin no H 25.2

pred

1 42.9

2 14.8

> xyplot(residuals ~ pred, ylab="residuals", cex=0.3,

xlab="predicted values", main="predicted vs. residuals",

type=c("p", "r", "smooth"), data=HELPrct)

a student’s guide to r 81

predicted vs. residuals

predicted values

re si

du al

s

−20

−10

0

10

20

20 30 40 50

●

●

●

●

●

●

●

●

●

●

●

●● ●

●

●

●

●

●

●

●● ●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

● ●

● ●

●

● ●

●

●

●

●

●

●

● ●

●

● ●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●●

●●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

●●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

● ●

●

● ●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●●

●

●

●

●

●

●

●

●

●

● ●

● ●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●● ●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

● ●

●

●●

● ●●

●

● ●

●

● ●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●●

●

●

●

● ●

●

●

●

● ●

●

●

●

●

● ●

●

● ●

●

●

●

●

●

●

● ●

●

●

●

● ●

●

> xyplot(residuals ~ mcs, xlab="mental component score",

ylab="residuals", cex=0.3,

type=c("p", "r", "smooth"), data=HELPrct)

mental component score

re si

du al

s

−20

−10

0

10

20

10 20 30 40 50 60

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

● ●●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

● ●

● ●

●

● ●

●

●

●

●

●

●

● ●

●

● ●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

● ●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●●

●

●

● ●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ● ●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

● ● ●

●

● ●

●

●●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

● ●

●

●

●

● ●

●

●

●

●

●●

●

● ●

●

●

●

●

●

●

● ●

●

●

●

● ●

●

The assumptions of normality, linearity and homoscedas- ticity seem reasonable here.

11 Probability Distributions & Random Variables

R can calculate quantities related to probability distribu- tions of all types. It is straightforward to generate ran- dom samples from these distributions, which can be used for simulation and exploration.

> xpnorm(1.96, mean=0, sd=1) # P(Z < 1.96)

If X ~ N(0,1), then

P(X <= 1.96) = P(Z <= 1.96) = 0.975

P(X > 1.96) = P(Z > 1.96) = 0.025

[1] 0.975

de ns

ity

0.1

0.2

0.3

0.4

0.5

−2 0 2

1.96 (z=1.96)

0.975 0.025

84 horton, kaplan, pruim

> # value which satisfies P(Z < z) = 0.975

> qnorm(.975, mean=0, sd=1)

[1] 1.96

> integrate(dnorm, -Inf, 0) # P(Z < 0)

0.5 with absolute error < 4.7e-05

A similar display is available for the F distribution.

> xpf(3, df1=4, df2=20)

[1] 0.957

de ns

ity

0.2

0.4

0.6

0.8

2 4 6

0. 95

7 0.

04 3

The following table displays the basenames for prob- ability distributions available within base R. These func- tions can be prefixed by d to find the density function for the distribution, p to find the cumulative distribution function, q to find quantiles, and r to generate random draws. For example, to find the density function of an ex- ponential random variable, use the command dexp(). The qDIST() function is the inverse of the pDIST() function, for a given basename DIST.

a student’s guide to r 85

Distribution Basename Beta beta

binomial binom Cauchy cauchy

chi-square chisq exponential exp

F f gamma gamma

geometric geom hypergeometric hyper

logistic logis lognormal lnorm

negative binomial nbinom normal norm Poisson pois

Student’s t t Uniform unif Weibull weibull

Digging Deeper

The fitdistr() within the MASS package facilitates estimation of parameters for many distribu- tions.

The plotDist() can be used to display distributions in a variety of ways.

> plotDist('norm', mean=100, sd=10, kind='cdf')

0.0

0.2

0.4

0.6

0.8

1.0

60 80 100 120 140

> plotDist('exp', kind='histogram', xlab="x")

86 horton, kaplan, pruim

x

D en

si ty

0.0

0.2

0.4

0.6

0 2 4 6 8

> plotDist('binom', size=25, prob=0.25, xlim=c(-1,26))

0.00

0.05

0.10

0.15

0 5 10 15 20 25

● ●

●

●

●

●

●

●

●

●

●

● ● ● ● ●

Multiple distributions can be plotted on the same plot.

> plotDist("binom", size=100, prob=.3, col="black", lwd=3, pch=16)

> plotDist("norm", mean=30, sd=sqrt(100 * .3 * .7),

groups = abs(x - 30) > 6 , type="h", under=TRUE)

a student’s guide to r 87

0.00

0.02

0.04

0.06

0.08

20 30 40

●●●●●● ●

● ●

●

●

●

●

●

● ●●●

●

●

●

●

●

●

● ●

● ●●●●●● ●

The plotFun() function can be used to plot an arbi- trary function (in this case an exponential random vari- able).

> f <- makeFun(2*exp(-2*x) ~ x) # exponential with rate parameter 2

> x <- rexp(1000, rate=2)

> histogram(~ x, width=0.2, center=0.1)

> plotFun(f(x) ~ x, col="red", lwd=2, add=TRUE)

x

D en

si ty

0.0

0.5

1.0

1.5

0 1 2 3 4

12 Power Calculations

While not generally a major topic in introductory courses, power and sample size calculations help to reinforce key ideas in statistics. In this section, we will explore how R can be used to undertake power calculations using ana- lytic approaches. We consider a simple problem with two tests (t-test and sign test) of a one-sided comparison.

We will compare the power of the sign test and the power of the test based on normal theory (one sample one sided t-test) assuming that σ is known. Let X1, ..., X25 be i.i.d. N(0.3, 1) (this is the alternate that we wish to calculate power for). Consider testing the null hypothesis H0 : µ = 0 versus HA : µ > 0 at significance level α = .05.

12.1 Sign test

We start by calculating the Type I error rate for the sign test. Here we want to reject when the number of pos- itive values is large. Under the null hypothesis, this is distributed as a Binomial random variable with n=25 tri- als and p=0.5 probability of being a positive value. Let’s consider values between 15 and 19.

> xvals <- 15:19

> probs <- 1 - pbinom(xvals, size=25, prob=0.5)

> cbind(xvals, probs)

xvals probs

[1,] 15 0.11476

[2,] 16 0.05388

[3,] 17 0.02164

[4,] 18 0.00732

[5,] 19 0.00204

90 horton, kaplan, pruim

> qbinom(.95, size=25, prob=0.5)

[1] 17

So we see that if we decide to reject when the number of positive values is 17 or larger, we will have an α level of 0.054, which is near the nominal value in the problem.

We calculate the power of the sign test as follows. The probability that Xi > 0, given that HA is true is given by:

> 1 - pnorm(0, mean=0.3, sd=1)

[1] 0.618

We can view this graphically using the command:

> xpnorm(0, mean=0.3, sd=1, lower.tail=FALSE)

If X ~ N(0.3,1), then

P(X <= 0) = P(Z <= -0.3) = 0.3821

P(X > 0) = P(Z > -0.3) = 0.6179

[1] 0.618

de ns

ity

0.1

0.2

0.3

0.4

0.5

−2 0 2 4

0 (z=−0.3)

0.3821 0.6179

The power under the alternative is equal to the proba- bility of getting 17 or more positive values, given that p = 0.6179:

> 1 - pbinom(16, size=25, prob=0.6179)

[1] 0.338

The power is modest at best.

a student’s guide to r 91

12.2 T-test

We next calculate the power of the test based on normal theory. To keep the comparison fair, we will set our α level equal to 0.05388.

> alpha <- 1-pbinom(16, size=25, prob=0.5); alpha

[1] 0.0539

First we find the rejection region.

> n <- 25; sigma <- 1 # given

> stderr <- sigma/sqrt(n)

> zstar <- qnorm(1-alpha, mean=0, sd=1)

> zstar

[1] 1.61

> crit <- zstar*stderr

> crit

[1] 0.322

Therefore, we reject for observed means greater than 0.322.

To calculate the power of this one-sided test we find the probability under the alternative hypothesis to the right of this cutoff.

> power <- 1 - pnorm(crit, mean=0.3, sd=stderr)

> power

[1] 0.457

The power of the test based on normal theory is 0.457. To provide a check (or for future calculations of this sort) we can use the power.t.test() function.

> power.t.test(n=25, delta=.3, sd=1, sig.level=alpha, alternative="one.sided",

type="one.sample")$power

[1] 0.441

92 horton, kaplan, pruim

This analytic (formula-based approach) yields a simi- lar estimate to the value that we calculated directly.

Overall, we see that the t-test has higher power than the sign test, if the underlying data are truly normal.

Calculating power empirically demonstrates the power of simulations.

13 Data Management

The Start Teaching with R book features an extensive section on data management, including use of the read.file() function to load data into R and RStudio.

Data management is a key capacity to allow students (and instructors) to “compute with data” or as Diane Lambert of Google has stated, “think with data”. We tend to keep student data management to a minimum during the early part of an introductory statistics course, then gradually introduce topics as needed. For courses where students undertake substantive projects, data manage- ment is more important. This chapter describes some key data management tasks.

The dplyr and tidyr packages provide an elegant approach to data management and fa- cilitate the ability of students to compute with data. Hadley Wickham, author of the pack- ages, suggests that there are six key idioms (or verbs) imple- mented within these packages that allow a large set of tasks to be accomplished: filter (keep rows matching criteria), select (pick columns by name), ar- range (reorder rows), mutate (add new variables), summarise (reduce variables to values), and group by (collapse groups). See http://www.amherst.edu/ ~nhorton/precursors for more details and resources.

13.1 Inspecting dataframes

The inspect() function can be helpful in describing the variables in a dataframe (the name for a dataset in R).

> inspect(iris)

categorical variables:

name class levels n missing

1 Species factor 3 150 0

distribution

1 setosa (33.3%), versicolor (33.3%) ...

quantitative variables:

name class min Q1 median Q3 max mean sd n missing

1 Sepal.Length numeric 4.3 5.1 5.80 6.4 7.9 5.84 0.828 150 0

2 Sepal.Width numeric 2.0 2.8 3.00 3.3 4.4 3.06 0.436 150 0

3 Petal.Length numeric 1.0 1.6 4.35 5.1 6.9 3.76 1.765 150 0

4 Petal.Width numeric 0.1 0.3 1.30 1.8 2.5 1.20 0.762 150 0

The iris dataframe includes one categorical and four quantitative variables.

94 horton, kaplan, pruim

13.2 Adding new variables to a dataframe

We can add additional variables to an existing dataframe using mutate(). But first we create a smaller version of the iris dataframe.

> irisSmall <- select(iris, Species, Sepal.Length)

> # cut places data into bins

> irisSmall <- mutate(irisSmall,

Length = cut(Sepal.Length, breaks=4:8))

Multiple commands can be chained together using the %>% (pipe) operator:

> irisSmall <- iris %>%

select(Species, Sepal.Length) %>%

mutate(Length = cut(Sepal.Length, breaks=4:8))

Note that in this usage the first argument to select() is the first variable (as it inherits the data from the previous pipe).

The cut() function has an op- tion labels which can be used to specify more descriptive names for the groups.

> head(irisSmall)

Species Sepal.Length Length

1 setosa 5.1 (5,6]

2 setosa 4.9 (4,5]

3 setosa 4.7 (4,5]

4 setosa 4.6 (4,5]

5 setosa 5.0 (4,5]

6 setosa 5.4 (5,6]

The CPS85 dataframe contains data from a Current Population Survey (current in 1985, that is). Two of the variables in this dataframe are age and educ. We can estimate the number of years a worker has been in the workforce if we assume they have been in the workforce since completing their education and that their age at graduation is 6 more than the number of years of educa- tion obtained. We can add this as a new variable in the dataframe using mutate().

a student’s guide to r 95

> CPS85 <- mutate(CPS85, workforce.years = age - 6 - educ)

> favstats(~ workforce.years, data=CPS85)

min Q1 median Q3 max mean sd n missing

-4 8 15 26 55 17.8 12.4 534 0

In fact this is what was done for all but one of the cases to create the exper variable that is already in the CPS85 data.

> tally(~ (exper - workforce.years), data=CPS85)

0 4

533 1

13.3 Dropping variables

Since we already have the exper variable, there is no rea- son to keep our new variable. Let’s drop it. Notice the clever use of the minus sign.

> names(CPS85)

[1] "wage" "educ" "race"

[4] "sex" "hispanic" "south"

[7] "married" "exper" "union"

[10] "age" "sector" "workforce.years"

> CPS1 <- select(CPS85, select = -matches("workforce.years"))

> names(CPS1)

[1] "wage" "educ" "race" "sex" "hispanic" "south"

[7] "married" "exper" "union" "age" "sector"

Any number of variables can be dropped or kept in a similar manner.

> CPS1 <- select(CPS85, select = -matches("workforce.years|exper"))

96 horton, kaplan, pruim

13.4 Renaming variables

The column (variable) names for a dataframe can be changed using the rename() function in the dplyr pack- age.

> names(CPS85)

[1] "wage" "educ" "race"

[4] "sex" "hispanic" "south"

[7] "married" "exper" "union"

[10] "age" "sector" "workforce.years"

> CPSnew = rename(CPS85, workforce=workforce.years)

> names(CPSnew)

[1] "wage" "educ" "race" "sex" "hispanic"

[6] "south" "married" "exper" "union" "age"

[11] "sector" "workforce"

The row names of a dataframes can be changed by simple assignment using row.names().

The faithful data set (in the datasets package, which is always available) has very unfortunate names.

It’s a good idea to establish practices for choice of variable names from day one.

> names(faithful)

[1] "eruptions" "waiting"

The measurements are the duration of an eruption and the time until the subsequent eruption, so let’s give it some better names.

> faithful <- rename(faithful,

duration = eruptions,

time.til.next=waiting)

> names(faithful)

[1] "duration" "time.til.next"

a student’s guide to r 97

> xyplot(time.til.next ~ duration, alpha=0.5, data=faithful)

duration

tim e.

til .n

ex t

50

60

70

80

90

2 3 4 5

If the variable containing a dataframe is modified or used to store a different object, the original data from the pack- age can be recovered using data().

> data(faithful)

> head(faithful, 3)

eruptions waiting

1 3.60 79

2 1.80 54

3 3.33 74

13.5 Creating subsets of observations

We can also use filter() to reduce the size of a dataframe by selecting only certain rows.

> data(faithful)

> names(faithful) <- c('duration', 'time.til.next')

> # any logical can be used to create subsets

> faithfulLong <- filter(faithful, duration > 3)

> xyplot( time.til.next ~ duration, data=faithfulLong )

98 horton, kaplan, pruim

duration

tim e.

til .n

ex t

70

80

90

3.0 3.5 4.0 4.5 5.0

●

●

●

●

● ●●

●

● ●

●●

●

●

●

● ● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

● ●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

● ●

●

●

●

● ●

● ●

●

●

●

● ●

●

●

●

●

●

● ● ●

●●

●

●

●●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●●

●

● ● ●

●

●

●

●

●

● ●

●

●

●●

● ●

●

●

●

●

●

●

●

●

● ●● ●

●

●

●

●

●

● ●

●

●●

● ●

●

●

● ●

●

●

●

●

● ●

●●

●

●

●

●

13.6 Sorting dataframes

Data frames can be sorted using the arrange() function.

> head(faithful, 3)

duration time.til.next

1 3.60 79

2 1.80 54

3 3.33 74

> sorted <- arrange(faithful, duration)

> head(sorted, 3)

duration time.til.next

1 1.60 52

2 1.67 64

3 1.70 59

Caution! It is usually better to make new datasets rather than modifying the original.13.7 Merging datasets

The fusion1 dataframe in the fastR package contains genotype information for a SNP (single nucleotide poly- morphism) in the gene TCF7L2. The pheno dataframe contains phenotypes (including type 2 diabetes case/control status) for an intersecting set of individuals. We can join

a student’s guide to r 99

(or merge) these together to explore the association be- tween genotypes and phenotypes using merge().

> require(fastR)

> require(dplyr)

> fusion1 <- arrange(fusion1, id)

> head(fusion1, 3)

id marker markerID allele1 allele2 genotype Adose Cdose Gdose Tdose

1 1002 RS12255372 1 3 3 GG 0 0 2 0

2 1009 RS12255372 1 3 3 GG 0 0 2 0

3 1012 RS12255372 1 3 3 GG 0 0 2 0

> head(pheno, 3)

id t2d bmi sex age smoker chol waist weight height whr sbp dbp

1 1002 case 32.9 F 70.8 former 4.57 112.0 85.6 161 0.987 135 77

2 1009 case 27.4 F 53.9 never 7.32 93.5 77.4 168 0.940 158 88

3 1012 control 30.5 M 53.9 former 5.02 104.0 94.6 176 0.933 143 89

> require(tidyr)

> fusion1m <- inner_join(fusion1, pheno, by='id')

> head(fusion1m, 3)

id marker markerID allele1 allele2 genotype Adose Cdose Gdose Tdose t2d bmi

1 1002 RS12255372 1 3 3 GG 0 0 2 0 case 32.9

2 1009 RS12255372 1 3 3 GG 0 0 2 0 case 27.4

3 1012 RS12255372 1 3 3 GG 0 0 2 0 control 30.5

sex age smoker chol waist weight height whr sbp dbp

1 F 70.8 former 4.57 112.0 85.6 161 0.987 135 77

2 F 53.9 never 7.32 93.5 77.4 168 0.940 158 88

3 M 53.9 former 5.02 104.0 94.6 176 0.933 143 89

Now we are ready to begin our analysis.

> tally(~t2d + genotype, data=fusion1m)

genotype

t2d GG GT TT

case 737 375 48

control 835 309 27

100 horton, kaplan, pruim

13.8 Slicing and dicing

The tidyr package provides a flexible way to change the arrangement of data. It was designed for converting be- tween long and wide versions of time series data and its arguments are named with that in mind.

The vignettes that accompany the tidyr and dplyr packages feature a number of useful examples of common data manipulations.

A common situation is when we want to convert from a wide form to a long form because of a change in per- spective about what a unit of observation is. For example, in the traffic dataframe, each row is a year, and data for multiple states are provided.

> traffic

year cn.deaths ny cn ma ri

1 1951 265 13.9 13.0 10.2 8.0

2 1952 230 13.8 10.8 10.0 8.5

3 1953 275 14.4 12.8 11.0 8.5

4 1954 240 13.0 10.8 10.5 7.5

5 1955 325 13.5 14.0 11.8 10.0

6 1956 280 13.4 12.1 11.0 8.2

7 1957 273 13.3 11.9 10.2 9.4

8 1958 248 13.0 10.1 11.8 8.6

9 1959 245 12.9 10.0 11.0 9.0

We can reformat this so that each row contains a mea- surement for a single state in one year.

> longTraffic <- traffic %>%

gather(state, deathRate, ny:ri)

> head(longTraffic)

year cn.deaths state deathRate

1 1951 265 ny 13.9

2 1952 230 ny 13.8

3 1953 275 ny 14.4

4 1954 240 ny 13.0

5 1955 325 ny 13.5

6 1956 280 ny 13.4

We can also reformat the other way, this time having all data for a given state form a row in the dataframe.

a student’s guide to r 101

> stateTraffic <- longTraffic %>%

select(year, deathRate, state) %>%

mutate(year=paste("deathRate.", year, sep="")) %>%

spread(year, deathRate)

> stateTraffic

state deathRate.1951 deathRate.1952 deathRate.1953 deathRate.1954 deathRate.1955

1 ny 13.9 13.8 14.4 13.0 13.5

2 cn 13.0 10.8 12.8 10.8 14.0

3 ma 10.2 10.0 11.0 10.5 11.8

4 ri 8.0 8.5 8.5 7.5 10.0

deathRate.1956 deathRate.1957 deathRate.1958 deathRate.1959

1 13.4 13.3 13.0 12.9

2 12.1 11.9 10.1 10.0

3 11.0 10.2 11.8 11.0

4 8.2 9.4 8.6 9.0

13.9 Derived variable creation

A number of functions help facilitate the creation or re- coding of variables.

13.9.1 Creating categorical variable from a quantita- tive variable

Next we demonstrate how to create a three-level categor- ical variable with cuts at 20 and 40 for the CESD scale (which ranges from 0 to 60 points).

> favstats(~ cesd, data=HELPrct)

min Q1 median Q3 max mean sd n missing

1 25 34 41 60 32.8 12.5 453 0

> HELPrct <- mutate(HELPrct, cesdcut = cut(cesd,

breaks=c(0, 20, 40, 60), include.lowest=TRUE))

> bwplot(cesd ~ cesdcut, data=HELPrct)

102 horton, kaplan, pruim

ce sd

0

10

20

30

40

50

60

[0,20] (20,40] (40,60]

●

●

●

The ntiles() function can be used to automate creation of groups in this manner.

It might be preferable to give better labels.

> HELPrct <- mutate(HELPrct, cesdcut = cut(cesd,

labels=c("low", "medium", "high"),

breaks=c(0, 20, 40, 60), include.lowest=TRUE))

> bwplot(cesd ~ cesdcut, pch="|", data=HELPrct)

ce sd

0

10

20

30

40

50

60

low medium high

The derivedFactor() function is even more general and can also be used for this purpose.

> HELPrct <- mutate(HELPrct,

anothercut = derivedFactor(

low = cesd >= 0 & cesd <= 20,

medium = cesd > 20 & cesd <= 40,

high = cesd > 40))

a student’s guide to r 103

13.9.2 Reordering factors

By default R uses the first level in lexicographic order as the reference group for modeling. This can be overriden using the relevel() function (see also reorder()).

> tally(~ substance, data=HELPrct)

alcohol cocaine heroin

177 152 124

> coef(lm(cesd ~ substance, data=HELPrct))

(Intercept) substancecocaine substanceheroin

34.373 -4.952 0.498

> HELPrct <- mutate(HELPrct, subnew = relevel(substance,

ref="heroin"))

> coef(lm(cesd ~ subnew, data=HELPrct))

(Intercept) subnewalcohol subnewcocaine

34.871 -0.498 -5.450

13.10 Group-wise statistics

It can often be useful to calculate summary statistics by group, and add these into a dataset. The group_by() function in the dplyr package facilitates this process. Here we demonstrate how to add a variable containing the median age of subjects by substance group.

> favstats(age ~ substance, data=HELPrct)

substance min Q1 median Q3 max mean sd n missing

1 alcohol 20 33 38.0 43.0 58 38.2 7.65 177 0

2 cocaine 23 30 33.5 37.2 60 34.5 6.69 152 0

3 heroin 19 27 33.0 39.0 55 33.4 7.99 124 0

> ageGroup <- HELPrct %>%

group_by(substance) %>%

summarise(agebygroup = mean(age))

> ageGroup

104 horton, kaplan, pruim

Source: local data frame [3 x 2]

substance agebygroup

(fctr) (dbl)

1 alcohol 38.2

2 cocaine 34.5

3 heroin 33.4

> nrow(ageGroup)

[1] 3

> nrow(HELPrct)

[1] 453

> HELPmerged <- left_join(ageGroup, HELPrct, by="substance")

> favstats(agebygroup ~ substance, data=HELPmerged)

substance min Q1 median Q3 max mean sd n missing

1 alcohol 38.2 38.2 38.2 38.2 38.2 38.2 0 177 0

2 cocaine 34.5 34.5 34.5 34.5 34.5 34.5 0 152 0

3 heroin 33.4 33.4 33.4 33.4 33.4 33.4 0 124 0

> nrow(HELPmerged)

[1] 453

13.11 Accounting for missing data

Missing values arise in almost all real world investiga- tions. R uses the NA character as an indicator for missing data. The HELPmiss dataframe within the mosaicData package includes all n = 470 subjects enrolled at baseline (including the n = 17 subjects with some missing data who were not included in HELPrct).

> smaller <- select(HELPmiss, cesd, drugrisk, indtot, mcs, pcs,

substance)

> dim(smaller)

[1] 470 6

> summary(smaller)

a student’s guide to r 105

cesd drugrisk indtot mcs pcs

Min. : 1.0 Min. : 0.00 Min. : 4.0 Min. : 6.8 Min. :14.1

1st Qu.:25.0 1st Qu.: 0.00 1st Qu.:32.0 1st Qu.:21.7 1st Qu.:40.3

Median :34.0 Median : 0.00 Median :37.5 Median :28.6 Median :48.9

Mean :32.9 Mean : 1.87 Mean :35.7 Mean :31.5 Mean :48.1

3rd Qu.:41.0 3rd Qu.: 1.00 3rd Qu.:41.0 3rd Qu.:40.6 3rd Qu.:57.0

Max. :60.0 Max. :21.00 Max. :45.0 Max. :62.2 Max. :74.8

NA's :2 NA's :14 NA's :2 NA's :2

substance

alcohol:185

cocaine:156

heroin :128

missing: 1

Of the 470 subjects in the 6 variable dataframe, only the drugrisk, indtot, mcs, and pcs variables have missing values.

> favstats(~ mcs, data=smaller)

min Q1 median Q3 max mean sd n missing

6.76 21.7 28.6 40.6 62.2 31.5 12.8 468 2

> with(smaller, sum(is.na(mcs)))

[1] 2

> nomiss <- na.omit(smaller)

> dim(nomiss)

[1] 453 6

> nrow(nomiss)

[1] 453

> ncol(nomiss)

[1] 6

> favstats(~ mcs, data=nomiss)

min Q1 median Q3 max mean sd n missing

6.76 21.8 28.6 40.9 62.2 31.7 12.8 453 0

106 horton, kaplan, pruim

Alternatively, we could generate the same dataset us- ing logical conditions.

> nomiss <- filter(smaller,

(!is.na(mcs) & !is.na(indtot) & !is.na(drugrisk)))

> dim(nomiss)

[1] 453 6

14 Health Evaluation (HELP) Study

Many of the examples in this guide utilize data from the HELP study, a randomized clinical trial for adult inpa- tients recruited from a detoxification unit. Patients with no primary care physician were randomized to receive a multidisciplinary assessment and a brief motivational intervention or usual care, with the goal of linking them to primary medical care. Funding for the HELP study was provided by the National Institute on Alcohol Abuse and Alcoholism (R01-AA10870, Samet PI) and National Institute on Drug Abuse (R01-DA10019, Samet PI). The details of the randomized trial along with the results from a series of additional analyses have been published1. 1 J. H. Samet, M. J. Larson, N. J.

Horton, K. Doyle, M. Winter, and R. Saitz. Linking alcohol and drug dependent adults to primary medical care: A ran- domized controlled trial of a multidisciplinary health inter- vention in a detoxification unit. Addiction, 98(4):509–516, 2003; J. Liebschutz, J. B. Savetsky, R. Saitz, N. J. Horton, C. Lloyd- Travaglini, and J. H. Samet. The relationship between sexual and physical abuse and substance abuse consequences. Journal of Substance Abuse Treatment, 22(3):121–128, 2002; and S. G. Kertesz, N. J. Horton, P. D. Friedmann, R. Saitz, and J. H. Samet. Slowing the revolving door: stabilization programs reduce homeless persons’ sub- stance use after detoxification. Journal of Substance Abuse Treat- ment, 24(3):197–207, 2003

Eligible subjects were adults, who spoke Spanish or English, reported alcohol, heroin or cocaine as their first or second drug of choice, resided in proximity to the pri- mary care clinic to which they would be referred or were homeless. Patients with established primary care rela- tionships they planned to continue, significant dementia, specific plans to leave the Boston area that would prevent research participation, failure to provide contact informa- tion for tracking purposes, or pregnancy were excluded.

Subjects were interviewed at baseline during their detoxification stay and follow-up interviews were under- taken every 6 months for 2 years. A variety of continuous, count, discrete, and survival time predictors and out- comes were collected at each of these five occasions. The Institutional Review Board of Boston University Medical Center approved all aspects of the study, including the creation of the de-identified dataset. Additional privacy protection was secured by the issuance of a Certificate of Confidentiality by the Department of Health and Human

108 horton, kaplan, pruim

Services. The mosaicData package contains several forms of the

de-identified HELP dataset. We will focus on HELPrct, which contains 27 variables for the 453 subjects with min- imal missing data, primarily at baseline. Variables in- cluded in the HELP dataset are described in Table 14.1. More information can be found at: http://www.amherst. edu/~nhorton/r2. A copy of the study instruments can be found at: http://www.amherst.edu/~nhorton/help.

Table 14.1: Annotated description of variables in the HELPrct dataset

VARIABLE DESCRIPTION (VALUES) NOTE age age at baseline (in years) (range

19–60) anysub use of any substance post-detox see also

daysanysub

cesd Center for Epidemiologic Stud- ies Depression scale (range 0–60, higher scores indicate more depres- sive symptoms)

d1 how many times hospitalized for medical problems (lifetime) (range 0–100)

daysanysub time (in days) to first use of any substance post-detox (range 0–268)

see also anysubstatus

dayslink time (in days) to linkage to primary care (range 0–456)

see also linkstatus

drugrisk Risk-Assessment Battery (RAB) drug risk score (range 0–21)

see also sexrisk

e2b number of times in past 6 months entered a detox program (range 1–21)

female gender of respondent (0=male, 1=female)

g1b experienced serious thoughts of suicide (last 30 days, values 0=no, 1=yes)

homeless 1 or more nights on the street or shelter in past 6 months (0=no, 1=yes)

a student’s guide to r 109

i1 average number of drinks (standard units) consumed per day (in the past 30 days, range 0–142)

see also i2

i2 maximum number of drinks (stan- dard units) consumed per day (in the past 30 days range 0–184)

see also i1

id random subject identifier (range 1–470)

indtot Inventory of Drug Use Conse- quences (InDUC) total score (range 4–45)

linkstatus post-detox linkage to primary care (0=no, 1=yes)

see also dayslink

mcs SF-36 Mental Component Score (range 7-62, higher scores are bet- ter)

see also pcs

pcs SF-36 Physical Component Score (range 14-75, higher scores are better)

see also mcs

pss_fr perceived social supports (friends, range 0–14)

racegrp race/ethnicity (black, white, his- panic or other)

satreat any BSAS substance abuse treat- ment at baseline (0=no, 1=yes)

sex sex of respondent (male or female) sexrisk Risk-Assessment Battery (RAB) sex

risk score (range 0–21) see also drugrisk

substance primary substance of abuse (alco- hol, cocaine or heroin)

treat randomization group (randomize to HELP clinic, no or yes)

Notes: Observed range is provided (at baseline) for con- tinuous variables.

15 Exercises and Problems

The first part of the exercise number indicates which chapter it comes from.

3.1. Using the HELPrct dataset, create side-by-side his- tograms of the CESD scores by substance abuse group, just for the male subjects, with an overlaid normal den- sity.

5.1. Using the HELPrct dataset, fit a simple linear regres- sion model predicting the number of drinks per day as a function of the mental component score. This model can be specified using the formula: i1 ∼ mcs. Assess the distribution of the residuals for this model.

10.1. The RailTrail dataset within the mosaic package includes the counts of crossings of a rail trail in Northamp- ton, Massachusetts for 90 days in 2005. City officials are interested in understanding usage of the trail network, and how it changes as a function of temperature and day of the week. Describe the distribution of the variable avgtemp in terms of its center, spread and shape.

> favstats(~ avgtemp, data=RailTrail)

min Q1 median Q3 max mean sd n missing

33 48.6 55.2 64.5 84 57.4 11.3 90 0

> densityplot(~ avgtemp, xlab="Average daily temp (degrees F)",

data=RailTrail)

112 horton, kaplan, pruim

Average daily temp (degrees F)

D en

si ty

0.00

0.01

0.02

0.03

20 40 60 80 100

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

10.2. The RailTrail dataset also includes a variable called cloudcover. Describe the distribution of this variable in terms of its center, spread and shape.

10.3. The variable in the RailTrail dataset that provides the daily count of crossings is called volume. Describe the distribution of this variable in terms of its center, spread and shape.

10.4. The RailTrail dataset also contains an indicator of whether the day was a weekday (weekday==1) or a week- end/holiday (weekday==0). Use tally() to describe the distribution of this categorical variable. What percentage of the days are weekends/holidays?

10.5. Use side-by-side boxplots to compare the distri- bution of volume by day type in the RailTrail dataset. Hint: you’ll need to turn the numeric weekday variable into a factor variable using as.factor(). What do you conclude?

10.6. Use overlapping densityplots to compare the distri- bution of volume by day type in the RailTrail dataset. What do you conclude?

10.7. Create a scatterplot of volume as a function of avgtemp using the RailTrail dataset, along with a regression line and scatterplot smoother (lowess curve). What do you observe about the relationship?

10.8. Using the RailTrail dataset, fit a multiple regres- sion model for volume as a function of cloudcover, avgtemp,

a student’s guide to r 113

weekday and the interaction between day type and aver- age temperature. Is there evidence to retain the interac- tion term at the α = 0.05 level?

10.9. Use makeFun() to calculate the predicted number of crossings on a weekday with average temperature 60

degrees and no clouds. Verify this calculation using the coefficients from the model.

> coef(fm)

(Intercept) cloudcover avgtemp weekday1

378.83 -17.20 2.31 -321.12

avgtemp:weekday1

4.73

10.10. Use makeFun() and plotFun() to display predicted values for the number of crossings on weekdays and weekends/holidays for average temperatures between 30 and 80 degrees and a cloudy day (cloudcover=10).

10.11. Using the multiple regression model, generate a histogram (with overlaid normal density) to assess the normality of the residuals.

10.12. Using the same model generate a scatterplot of the residuals versus predicted values and comment on the linearity of the model and assumption of equal variance.

10.13. Using the same model generate a scatterplot of the residuals versus average temperature and comment on the linearity of the model and assumption of equal variance.

11.1. Generate a sample of 1000 exponential random vari- ables with rate parameter equal to 2, and calculate the mean of those variables.

11.2. Find the median of the random variable X, if it is exponentially distributed with rate parameter 10.

12.1. Find the power of a two-sided two-sample t-test where both distributions are approximately normally distributed with the same standard deviation, but the group differ by 50% of the standard deviation. Assume

114 horton, kaplan, pruim

that there are 25 observations per group and an alpha level of 0.054.

12.2. Find the sample size needed to have 90% power for a two group t-test where the true difference between means is 25% of the standard deviation in the groups (with α = 0.05).

13.1. Using faithful dataframe, make a scatter plot of eruption duration times vs. the time since the previous eruption.

13.2. The fusion2 data set in the fastR package contains genotypes for another SNP. Merge fusion1, fusion2, and pheno into a single data frame.

Note that fusion1 and fusion2 have the same columns.

> names(fusion1)

[1] "id" "marker" "markerID" "allele1" "allele2" "genotype" "Adose"

[8] "Cdose" "Gdose" "Tdose"

> names(fusion2)

[1] "id" "marker" "markerID" "allele1" "allele2" "genotype" "Adose"

[8] "Cdose" "Gdose" "Tdose"

You may want to use the suffixes argument to merge() or rename the variables after you are done merging to make the resulting dataframe easier to navigate.

Tidy up your dataframe by dropping any columns that are redundant or that you just don’t want to have in your final dataframe.

16 Bibliography

[BcRB+ 14] B.S. Baumer, M. Çetinkaya Rundel, A. Bray,

L. Loi, and N. J. Horton. R Markdown: In- tegrating a reproducible analysis tool into introductory statistics. Technology Innovations in Statistics Education, 8(1):281–283, 2014.

[HBW15] N.J. Horton, B.S. Baumer, and H. Wickham. Setting the stage for data science: integration of data management skills in introductory and second courses in statistics (http://arxiv. org/abs/1401.3269). CHANCE, 28(2):40–50, 2015.

[KHF+ 03] S. G. Kertesz, N. J. Horton, P. D. Friedmann,

R. Saitz, and J. H. Samet. Slowing the re- volving door: stabilization programs reduce homeless persons’ substance use after detox- ification. Journal of Substance Abuse Treatment, 24(3):197–207, 2003.

[LSS+02] J. Liebschutz, J. B. Savetsky, R. Saitz, N. J. Hor- ton, C. Lloyd-Travaglini, and J. H. Samet. The relationship between sexual and physi- cal abuse and substance abuse consequences. Journal of Substance Abuse Treatment, 22(3):121– 128, 2002.

[MM07] D. S. Moore and G. P. McCabe. Introduction to the Practice of Statistics. W.H.Freeman and Company, 6th edition, 2007.

116 horton, kaplan, pruim

[NT10] D. Nolan and D. Temple Lang. Computing in the statistics curriculum. The American Statistician, 64(2):97–107, 2010.

[RS02] F. Ramsey and D. Schafer. Statistical Sleuth: A Course in Methods of Data Analysis. Cengage, 2nd edition, 2002.

[SLH+ 03] J. H. Samet, M. J. Larson, N. J. Horton,

K. Doyle, M. Winter, and R. Saitz. Linking alcohol and drug dependent adults to primary medical care: A randomized controlled trial of a multidisciplinary health intervention in a detoxification unit. Addiction, 98(4):509–516, 2003.

[Tuf01] E. R. Tufte. The Visual Display of Quantitative Information. Graphics Press, Cheshire, CT, 2nd edition, 2001.

[Wor14] ASA Undergraduate Guidelines Workgroup. 2014 curriculum guidelines for undergradu- ate programs in statistical science. Techni- cal report, American Statistical Association, November 2014. http://www.amstat.org/ education/curriculumguidelines.cfm.

17 Index

%>%, 94

abs(), 79, 80

add option, 78, 87

adding variables, 94

all.x option, 99

alpha option, 49, 62

analysis of variance, 65

anova(), 77

anova(), 66, 70

aov(), 66, 75

arrange(), 98, 99

auto.key option, 76, 78

band.lwd option, 53

binom.test(), 40

binomial test, 40

bootstrapping, 36

breaks option, 101

bwplot(), 61, 66, 75

by.x option, 99

categorical variables, 39

center option, 87

cex option, 45, 62, 80

chisq.test(), 43, 58

class(), 48

coef(), 47, 48, 103

coefficient plots, 79

col option, 34

conf.int option, 73

confint(), 36, 41, 48

contingency tables, 39, 55

Cook’s distance, 51

cor(), 46

correct option, 41

correlation, 46

Cox proportional hazards model, 74

coxph(), 74

CPS85 dataset, 94

creating subsets, 97

cross classification tables, 55

CrossTable(), 56

cut(), 94, 101

data management, 93

data(), 97

dataframe, 94

dataframes inspecting, 93

merging, 98

reshaping, 100

sorting, 98

subsetting, 97

density option, 48

densityplot(), 34

derived variables, 101

derivedFactor(), 102

diffmean(), 64

dim(), 104

display first few rows, 94

dnorm(), 83

do(), 37

dotPlot(), 33

dplyr package, 30, 31

dropping variables, 95

exp(), 69

factor reordering, 103

factor(), 67

failure time analysis, 73

faithful dataset, 96

family option, 69

favstats(), 28, 61, 103, 105

filter(), 30, 95

Fisher’s exact test, 59

fisher.test(), 59

fit option, 31

fitted(), 79

format option, 30

freqpolygon(), 34

function(), 46

fusion1 dataset, 99

gather(), 100

glm(), 69

grid.text(), 34

group-wise statistics, 103

group_by(), 103

groups option, 65, 78, 86

head(), 94

Health Evaluation and Linkage to Primary Care study, 107

118 horton, kaplan, pruim

HELP study, 107

HELPmiss dataset, 104

HELPrct dataset, 27

histogram(), 29, 87

honest significant differences, 67

include.lowest option, 101

incomplete data, 104

inspect(), 93

inspecting dataframes, 93

install.packages(), 14

installing packages, 14

integrate(), 83

interactions, 77

iris dataset, 94

is.na(), 105

Kaplan-Meier plot, 73

knitr, 14

labels option, 67

ladd(), 34

layout option, 32

left_join(), 103

levels option, 67

leverage, 52

linear regression, 47

linearity, 45

lm(), 47, 67

loading packages, 14

logistic regression, 69

lowess, 45

lty option, 34

lwd option, 34, 45, 86

makeFun(), 78, 87

margins option, 39

markdown, 14

mean(), 27, 61

median(), 28, 75

merging dataframes, 98

missing data, 104

model comparison, 67

Modeling with R, 27

mosaic package, 27

mosaicplot(), 57

mplot(), 50, 68, 79

msummary(), 76

msummary(), 48, 63, 69

multiple comparisons, 67

multiple regression, 76

multivariate relationships, 76

mutate(), 67, 94, 101, 103

NA character, 104

na.omit(), 105

names(), 96

ncol(), 105

nint option, 32

nrow(), 103, 105

ntiles(), 102

oddsRatio(), 56

one-way ANOVA, 65

options(), 27

pairs plot, 47

panel.abline(), 34

panel.labels(), 46

panel.lmbands(), 53

panel.mathdensity(), 34

panel.text(), 46

paste(), 101

pbinom(), 90

pch option, 45, 86

pchisq(), 43

Pearson correlation, 46

permutation test, 64

pipe operator, 94

plotDist(), 43, 85

plotFun(), 78, 87

pnorm(), 83

polygons, 34

prediction bands, 53

print(), 41

prop.test(), 41

proportional hazards model, 74

pval(), 41

qdata(), 37

qnorm(), 83

qqmath(), 49

quantile(), 28

quantiles, 28

random variables, 83

read.file(), 93

regression, 47

regression diagnostics, 79

relevel(), 103

rename(), 96

renaming variables, 96

reordering factors, 103

reproducible analysis, 14

require(), 14, 27

resample(), 36

resampling, 36, 64

reshaping dataframes, 100

resid(), 79

residual diagnostics, 79

residuals(), 48

rexp(), 87

rnorm(), 83

row.names(), 96

rownames(), 46

rsquared(), 48

RStudio.Version(), 19

scale versus location, 51

scatterplot matrix, 47

scatterplots, 45

sd(), 28

select option, 95

select(), 94, 101, 103, 104

sessionInfo(), 19

shuffle(), 64

significance stars, 48

smoothers, 45

sorting dataframes, 98

Spearman correlation, 46

splom(), 47

spread(), 101

Start Modeling with R, 27

Start Teaching with R, 27

a student’s guide to r 119

stem(), 30

subsetting dataframes, 97

sum(), 43, 105

summarise(), 103

summary(), 63

Surv(), 73

survfit(), 73

survival analysis, 73

t.test(), 36, 62

tables, 39, 55

tally(), 30, 39, 55, 103

Teaching with R, 27

test option, 70

thinking with data, 93

tidyr package, 31, 99

time to event analysis, 73

transforming dataframes, 100

transposing dataframes, 100

Tukey’s HSD, 67

TukeyHSD(), 67

type option, 45, 80, 86

under option, 86

var(), 28

var.equal option, 62

vignettes, 13

which option, 50

width option, 33, 87

wilcox.test(), 63

with(), 27, 105

xchisq.test(), 44, 58

xlab option, 73

xlim option, 65

xpnorm(), 83

xyplot(), 45, 62, 78

ylab option, 78

Stat_205_hw_r_1.pdf