Install R

This page has 30 images

Download pdf of this page

Introduction

The first time installing R can seem intimidating. To start, be clear about the overall goal of the procedure: providing the student with an accessible environment for solving statistics problems.

In brief, this page explains how to get R set up on your computer. First, you need to download the R installer from the official CRAN website. When you run the installer,in general, accept the default choices. However, for Windows users, it’s important to right-click the file and choose “Run as administrator.” This step ensures that R has the proper permissions to install correctly and avoids problems with user access later on. Once installed, you can open R and test it by typing a simple command like `2 + 2` in the console to confirm everything is working.

The flow chart presented in Fig 1 suggests a way to orient the student to solving the task. First, understand which type of computing environment you have available. Second, if a macOS, then additional software (XQuartz) is required to provide full function to the R software. Third, download and install the base R software appropriate for the computer system.

Figure 1 suggests an R installation flow chart.

Figure 1. Suggested flow chart for R installation. 

It’s possible to follow the steps in Figure 1, accepting all default options presented along the way, to end up with a working R environment. As with many software processes, there are choices beyond defaults that can be made to improve the software use.

This page presents a detailed guide about how to install R onto your computer — this is referred to as building a local development environment or LDE. Additional install R help was provided in Chapter 1.1 – A quick look at R and R Commander.

Instructions for RStudio are also provided (optional for BI311 students). A guide to install R Commander is provided in Install R Commander.

Instructions for how to run R via a “cloud computing” (serverless) option — a remote development environment — are also provided, Use R in the Cloud.

For help upgrading installed packages after upgrading new R version, see R packages.

Run R on your computer (LDE install)

CoLab, skip this step: Instead, go to Use R in the cloud.

So why in this day in age should you install and build R on your own computer? The remote options to run R in the cloud are a wonderful option, convenient: you can access anywhere you have internet, from any device that connects to the internet. It’s easy to share and work together on projects, particularly those based on Jupyter Notebooks.

I think the main benefits to a local installation is it’s a more efficient environment to work in — you have control of everything and, provided your PC has power, a working R install on your computer will always be available to you. Since you can control the update cycle for your computer, you won’t run into times you cannot access the remote server to work on your project. Testing code is faster on a local install, feedback — think error messages — apply to your installed version. And, while remote R servers may come with low initial costs to students, any significant use will quickly require paid accounts. As a reminder, the good folks at the R-project continue to offer R as free software. All you need to do is work through the install process.

Start at the R-project homepage, r-project.org. To download software, first click on CRAN link,  located on left hand side of the screen (here, highlighted by green arrow, Fig 2).

Figure 2. Screenshot homepage for R-project.org. 

Figure 3 shows a screenshot of the CRAN mirror page. The idea is to select the mirror site closest to your location. In Hawaiʻi, that’s likely to be any of the sites in California. However, I recommend selecting the first in the list, 0-cloud, at cloud.r-project.org (highlighted by green arrow, Fig 3).

Figure 3. Screenshot of portion of R-Project CRAN mirror page.

After selecting the mirror site, the download page is presented (Fig 4). Click on the link that corresponds to your computer system (Linux, macOS, or Windows).

Figure 4. Screenshot of portion of base R download page.

Once the installation file is located onto your computer, proceed to install base R.

Detailed instructions

For screenshots of installation steps on WinPC, see Win10/11 setup, screenshots

1. Windows PCs, download the base application from selected CRAN Mirror site, select Download R for Windows, and install the R software as you would any other software. All of you are likely to have the 64-bit version of Windows 11, so install the 64-bit version of R. Follow the instructions as they are presented. Screenshots of the install process are available at the end of this page (click here or scroll down to Win11 setup, Screenshots).
   • Current versions of Microsoft Windows come in several flavors, the simplest distinction is between home and pro. R runs perfectly well on both.
     • Windows 10 is reaching end of life cycle.
     • Some inexpensive Microsoft Windows PCs are built on ARM64, not Intel or AMD64 CPU. Thus, installing R and or RStudio may prove problematic.
     • Also note: Windows in S mode only run applications from the Microsoft store. To install R, you first may need to switch out of S mode (see Microsoft FAQ about S mode).
   • You should install R with Administrator privileges. Highlight the install file, right-click the file, and select “Run as administrator” from the popup menu.
   • When you first try to run R you may get a popup screen “Windows protected your PC,” locate and click on the “More info” link and select “Run anyway.”
     • This in no way will harm your computer — provided you have downloaded from official sites. R is a verified program. Microsoft has taken an aggressive line on developers and favors apps that are part of their app store.
     • It is advisable to confirm for yourself: check the md5sum against the fingerprint on the CRAN server
   • When prompted, I recommend that you change the install directory to root folder, e.g., C:\R\R-4.5.1. This will allow for installation of packages to the common library as opposed to a personal library.
     • I recommend this change because of how Windows assigns home folders. During initial setup Windows 10 prompted you to choose a username and whether you wanted your work stored locally or in your OneDrive folder. A worse case scenario? You select a username with spaces, e.g.,”Mike Dohm,” and you selected OneDrive. Both will cause challenges later for running and or installing packages for R.
       • If you install R anywhere but the default Program files folder on your Win10/11 PC, chances are you will need to add the folder containing the executable, r.exe, to you path.
         • Search: “env”
         • Open “Edit the system environment variables” in the Control Panel
         • Click Advanced tab, then click on Environment Variables… button (lower right of panel)
         • Under System variables, scroll to and select Path.
         • Click Edit… button, then click New button.
         • Type in the path to the folder containing R.exe. That’s likely to be C:\R\R-4.5.1, assuming R-4.5.1 is the latest version of R installed on your computer.
     • I made a video for you. Video is about 26 minutes long; at 22 minute mark, video includes how to install R Commander (instructions provided Install R Commander).
       • Apologies for the production quality — videos are not my thing.

https://youtu.be/upjmBieh3bM

For screenshots of installation steps on WinPC, see MacOS setup, screenshots

1. macOS PCs, first you must download and install XQuartz from https://www.xquartz.org. Best to restart your mac after installing XQuartz then proceed to install R.
   After installing XQuartz, then return to https://cran.r-project.org, select Download for Mac(OS) X, and run the installer. Screenshots of the install process are available at the end of this page (click here or scroll down to Macos setup, Screenshots).
   • As of August 2025, be advised that there are two distinct R versions for your MacBook or iMac.
     • For MacBook or iMac with Apple’s M1 or M2 ARM chip sets, download and install R-4.5.1-arm64.pkg.
       • If you recently purchased a new MacBook or iMac (2020 to present), then you probably have the M2 or M3 chipset (check by clicking the Apple icon, then selecting About this Mac or System Information (/Applications/Utilities/System Information.app)).
       • XQuartz version 2.8.5 works on macs with either the M1 or Intel chipsets.
     • For older MacBook or iMacs with Intel processors, download R-4.5.1-x86_64.pkg.
     • Depreciated 8/4/2021: Be advised that these instructions are for Intel-based macs. At the time of writing these instructions (April 2021), the installation of XQuartz and R should work on new M1-based macs. At the time of this writing (April 2021), however, R will not run natively on your M1 mac. It will run using Rosetta 2, an emulator that is included with your M1 mac. The R folks are busy working on a version that will run natively, which may be ready within a few months.
   • At the completion of the install process, don’t forget to drag the R app to your Applications folder.
   • When you first try to run R, you may get a popup screen which provides no option to start the app, and perhaps even a rather ominous option to move the app to Trash. Just close the warning message and right-click on the R app. A new screen pops up, which looks very much like the previous warning, but now you will see and option to open the app. Click on open to start R.
   • Like the message to Windows PC users, bypassing Apple’s Gatekeeper to run R in no way will harm your computer — provided you have downloaded from official sites — R is a verified program. Apple has taken an aggressive line on developers and favors apps that are part of their app store.
2. LINUX distros. If your PC platform is Linux, then you should be comfortable with installation and updating of software. R base is already included in Debian distributions (e.g., Mint, Ubuntu). See https://cloud.r-project.org/ for additional instructions.
   • For Chromebook users, if you can install a Linux subsystem, then you can also install and run R. For instructions to install R see Levi’s excellent writeup at levente.littvay.hu/chromebook/.

Test R

For both macOS and Windows PCs, successful installation of R on your computer installs base R programming language and a simple graphical user interface. Test your install by running code in the terminal (one line at a time) or via script:

1. Rgui.exe (Windows PC)
2. File → New script

Enter code in script editor, e.g.,

myX <- c(1,2,3,4)
myY <- c(5,10,15,20)
plot(myY,myX)

3. Run code: Ctrl+R

Figure 5 shows screenshot of R running on WinPC. RGui.exe (1), script editor (2), and results of plot() (3) on WinPC.

Figure 5. Screenshot of RGui.exe (1), script editor (2), and results of plot() (3) on WinPC.

1. R.app (macOS): run code in the terminal or via script
2. File → New Document

Enter code in script editor, e.g.,

myX <- c(1,2,3,4)
myY <- c(5,10,15,20)
plot(myY,myX)

3. Run code: Cmd+Enter

Figure 6 shows screenshot of R running on macOS. R.app (1), script editor (2), and results of plot() (3).

Figure 6. Screenshot of R.app (1), script editor (2), and results of plot() (3) on macOS.

Many of you would like a video. Do a little search and you’ll find plenty, although most are also showing how to install RStudio in addition to base R.

For BI311, we also use R Commander package

R Commander is a package that adds function to R; it provides a familiar point-and-click interface to R, which allows the user to access functions via a drop-down menu system (Fox 2017).

Go to Install R Commander guide.

Personalize R startup

At startup, R app looks to a series of environment files to instruct how it will provide a working session environment. One important file to know about is the .Rprofile file, which will reside in your home directory on your computer. This file generally won’t be present after the first installation, but can be easily added.

First, verify the home directory where R will first look for an .Rprofile file.

Identify home directory

At the R prompt, verify your home directory, type and run

path.expand("~")

On a macOS, this should return something like

"~Users/USERNAME"

On a Windows 11 machine, this may return

"C:/Users/USERNAME/OneDrive/Documents"

if default instructions were accepted during initial Windows setup.

How to modify .Rprofile

Good practice suggests setting a project folder to contain work done for a project. For biostatistics students, minimally, this might be the course name, BIOSTATS. Within this folder we’d expect to find all course data files, homework, R coding scripts and related for the class. While it is easy enough to point R to a data file (see Part 07. Working with your own data, Mike’s Workbook for Biostatistics), it’s often easier to tell R in advance where to expect to find these files by setting a working folder at the beginning of a session using the setwd() command.

Alternatively, you can edit .Rprofile so that R always starts in the folder you want.

To open (or create if the first time) .Rprofile, at the R prompt type and submit

file.edit("~/.Rprofile")

options(test_success = "Success! .Rprofile is loaded.")

getOption("test_success")

Windows 11 Setup

• Open your local home documents folder (usually C:\Users\YOUR_NAME\Documents).
• Create or open a file named exactly .Rprofile (Ensure it does not end in .txt).

Copy this code for Windows:

.First <- function() {
# 1. Grab the literal, local hard drive user folder
local_root <- Sys.getenv("USERPROFILE")

# 2. Open a native Windows browse window for folder selection
message("Please select your local project folder...")
flush.console()

chosen_path <- choose.dir(
default = local_root,
caption = "Select your local project folder"
)

# 3. Set the directory to what they clicked on
if (!is.na(chosen_path)) {
setwd(chosen_path)
message("Working directory successfully set to: ", chosen_path)
} else {
warning("No folder selected. Working directory unchanged.")
}
}

macOS Setup

• Ensure R.app has been granted Full Disk Access in your Mac’s System Settings > Privacy & Security.
• Open Terminal and run open -a TextEdit ~/.Rprofile to open your profile configuration.

Copy this code for Mac:

.First <- function() {
# 1. Open a macOS Finder folder selection window
message("Please select your local project folder...")
flush.console()

# Calls native macOS AppleScript to open a secure folder picker
chosen_path <- tryCatch({
system("osascript -e 'POSIX path of (choose folder with prompt \"Select your local project folder:\")'", intern = TRUE)
}, error = function(e) {
return(NA)
})

# 2. Set the directory to what they clicked on
if (!is.na(chosen_path) && length(chosen_path) > 0) {
setwd(chosen_path)
message("Working directory successfully set to: ", chosen_path)
} else {
warning("No folder selected. Working directory unchanged. Keep in mind R.app needs 'Full Disk Access' in macOS System Settings.")
}
}

Set your CRAN mirror site

To add function to R, you’ll download many R packages. These, too, are stored at the mirror sites. Each time you run install.packages(), R will ask for the mirror site. The R app on Windows pops up a long selection box; macOS R.app is a bit more friendly, but regardless, it’s best to choose a mirror site at the beginning of a session. At the R prompt, type and submit the command

options(repos = c(CRAN = "https://cloud.r-project.org"))

this will remain in effect for the current R session.

To make the change permanent, add to .Rprofile in your personal home directory

local({
r <- getOption("repos")
r["CRAN"] <- "https://cloud.r-project.org"
options(repos = r)
})

If you use textedit or notepad, make sure to return after the }); the file requires an empty new line to run.

Restart R, then enter at the prompt

file.exists("~/.Rprofile")

If you’ve created the Rprofile file correctly, R should repeat.

Run R from the terminal

Whether yours is a mac or win11 pc, you have a powerful computing environment lurking beneath the glossy graphic user interface. It’s called the terminal. The terminal is a place where text-based instructions may be written and submitted to command your computer to do something.

On win11 pc, the modern terminal is the Power Shell. Search “PowerShell,” then open Windows PowerShell. Alternatively, shortcut + X, then I key. This opens terminal in your home folder.

• You’ll want to navigate to your working folder, e.g., BI311 on the Desktop. Setting your working folder can be done in R of course, but at the terminal type the command “cd \users\default\desktop\bi311“.

On macOS, Spotlight (search) “Terminal.app“. Alternatively, right-click on your working folder icon and select New terminal at Folder from the popup menu.

On Ubuntu Linux, shortcut keys + alt, then t key.

Open the included script editor

Use RStudio

Rstudio is a very popular data science IDE, specifically written for R programming. Therefore, if your goal is to become adept at R programming, go with RStudio. Like R Commander, Rstudio provides an environment to write an manage code, generate markdown reports, and importantly, manage files in R projects. Figure 7 shows a screenshot of RStudio: Clockwise, four panes: Source, Environment, Output, Console on WinPC. BI311 students, please note that I recommend use of R Commander and Google CoLab.

Figure 7. Screenshot of RStudio IDE.

Installation of RStudio desktop is straightforward.

Go to https://posit.co/download/rstudio-desktop/

Install R first. Since at this point you’ve likely already installed R, go directly to 2 (green box, Fig 8).

Figure 8. Screenshot RStudio desktop download page, current as of September 2025.

After downloading the installation software, double click the file to begin. Follow on screen instructions for your computer system.

Run R in the “Cloud”

If you do not wish to install R, or, if you have a Chromebook and, therefore cannot gracefully install R, then there are alternatives; Run R in the Cloud. I’ll list three ways to run R in the cloud for free. Go to Use R in the Cloud guide.

MacOS setup, Screenshots

Download R install package from R-project.org, then select the R install package from your Download folder (Fig 9).

Figure 9. Screenshot — Find the R install file in your download folder.

First screen, R install for macOS (Fig 10). Select continue.

Figure 10. Screenshot: first instructions pop-up window.

Second screen, R install on macOS (Fig 11). Select continue.

Figure 11. Screenshot: second instructions pop-up window.

Third screen, R install on macOS (Fig 12). Select continue.

Figure 12. Screenshot: third instructions pop-up window.

Fourth screen, R install on macOS (Fig 13). Agree to continue.

Figure 13. Screenshot: fourth instructions pop-up window.

Fifth screen, R install on macOS (Fig 14). Select Install.

Figure 14. Screenshot: fifth instructions pop-up window.

Sixth screen, R install on macOS (Fig 15). Enter your username and password for your computer, then select Install Software.

Figure 15. Screenshot: sixth instructions pop-up window.

Seventh screen, R install on macOS (Fig 16). Several screens will popup, reporting progress. Be patient!

Figure 16. Screenshot: seventh instructions pop-up window.

Eighth and final screen, R install on macOS (Fig 17). Select close.

Figure 17. Screenshot: eighth and final instructions pop-up window.

Optional — Keep or discard the install file (Fig 18). I keep and then do manual delete after I’ve confirmed the installation.

Figure 18. Screenshot: MacOS will prompt with option to delete the installation file. This has no effect on the installation of R.

From Applications folder, start R.app. You should see the R Console (Fig 19).

Figure 19. Screenshot: R Console in the R.app on macOS.

Wind10/11 setup, screenshots

Download from R-project.org, then right-click the R install package from your Download folder. Run as administrator (Fig 20).

Figure 20. WinPC  screenshot: Pop-up menu after right-click on the R installation file in windows Explorer.

First screen, select language (Fig 21). Select OK to continue

Figure 21. WinPC screenshot: first instructions pop-up window.

Second screen (Fig 22), click Next to continue

Figure 22. WinPC screenshot: second instructions pop-up window.

Third screen (Fig 23). Change the default location (show in the screenshot) to root folder, e.g., C:\R\R-4.5.1 (current version as of August 2025).

Figure 23. WinPC screenshot: third instructions pop-up window. Note — R-4.5.1 (current version as of August 2025).

Fourth screen (Fig 24). Change startup options. Select Yes (customized startup) to continue.

Figure 24. WinPC screenshot: fourth instructions pop-up window.

Fifth screen (Fig 25), select SDI — single document interface — not the default MDI — multiple document interface, then Next to continue.

Figure 25. WinPC screenshot: fifth instructions pop-up window. Recommend setting to SDI.

Sixth screen (Fig 26), select HTML help, then Next to continue.

Figure 26. WinPC screenshot: sixth instructions pop-up window.

Seventh screen (Fig 27), leave start menu folder as is (R), then Next to continue.

Figure 27. WinPC screenshot: seventh instructions pop-up window.

Eighth screen (Fig 28), check all boxes, then Next to continue.

Figure 28. WinPC screenshot: eighth instructions pop-up window.

Ninth screen (Fig 29), a series of status updates during the installation.

Figure 29. WinPC screenshot: ninth instructions pop-up window.

Final screen (Fig 30), successful install.

Figure 30. WinPC screenshot: final instructions pop-up window — successful installation.

/MD

20.11 – Plot a Newick tree

Introduction

The phrase “paradigm shift”, attributed to Kuhn (1962, see Wikipedia), may be well-worn and even abused today (Naughton 2012), but the shift in thinking from essential types and group thinking (essentialism) to viewing species as varying individuals in populations (populating thinking) revolutionized biology (O’Hara 1997, Sandvik 2008). Tree thinking is the manifestation of Charles Darwin’s “descent with modification” metaphor (Gregory 2008). Thus, every biology student should have ability to work with, and interpret, phylogenetic trees (tree thinking). The subject of creating and working with phylogenetic graphs is complicated with an extensive library. A good review is available from Holder and Lewis (2003) and readers should know Felsenstein’s book (2004).

Here, I include a modest, incomplete primer on working with trees in R.

• Loading the tree file
• Change tip names
• Write tip names to a text file
• Plot the tree as phylogram or cladogram
• Get node labels
• Re-root the tree
• Write a tree to a file

I assume that the student already has a set of species or other taxa, gathered sequences (DNA or protein), aligned the sequences, estimated gene or phylogeny tree, and wishes to view and manipulate the tree in R. While these kinds of analyses can be done with R and R packages (see Task view: Phylogenetics), other software may be better choice for the student just beginning with phylogenetic tree building (see Unipro UGENE and MEGA, for examples). If the goal is just to view a tree file, or add annotations, then I recommend the iTOL tools.

Data formats

Phylogeny and gene trees are special cases of network graphs. Newick format (Wikipedia) is a common but limited representation of the tree which uses parentheses (groupings) and commas (branching). Other formats permit additional information; examples are Nexus file (Wikipedia) and the extension of Nexus to XML, NeXML (Wikipedia), and phyloXML (Wikipedia) formats. Our example uses Newick format.

Data set

I’ll use a “time tree” for an example. Tree from timetree.org, list of species (copy/paste list to a text file, load the text file Load list of of Species, then save the tree as a Newick file).

Alligator mississippiensis
Felis catus
Bos taurus
Gallus gallus
Pan troglodytes
Canis lupus
Homo sapiens
Anolis carolinensis
Macaca mulatta
Mus musculus
Didelphis virginiana
Sus scrofa
Oryctolagus cuniculus
Rattus norvegicus

R code

Requires the ape package. Phylotools and Phytools packages provide additional handy functions. References for these packages are listed at the end of this page.

require(ape)
require(phytools)
require(phylotools)
#If tree file, then
read.tree(file="tree14.nwk")
or
tree14 <- read.tree(file.choose())
#If no tree file saved, copy the Newick data use text="", replace example tree with your Newick tree
tree14 <-read.tree(text="((Anolis_carolinensis:279.65697667,(Gallus_gallus:236.50266286,Alligator_mississippiensis:236.50266286)'14':43.15431381)'13':32.24694470,(Didelphis_virginiana:158.59758758,(((Felis_catus:54.32144118,Canis_lupus:54.32144118)'11':23.43351523,(Bos_taurus:61.96598852,Sus_scrofa:61.96598852)'10':15.78896789)'19':18.70743276,((Oryctolagus_cuniculus:82.14079889,(Rattus_norvegicus:20.88741740,Mus_musculus:20.88741740)'9':61.25338149)'22':7.68238853,(Macaca_mulatta:29.44154682,(Pan_troglodytes:6.65090500,Homo_sapiens:6.65090500)'8':22.79064182)'6':60.38164060)'30':6.63920175)'29':62.13519841)'27':153.30633379);")
#return information about the object
tree14

Output returned by R

Phylogenetic tree with 14 tips and 13 internal nodes.
Tip labels:
Anolis_carolinensis, Gallus_gallus, Alligator_mississippiensis, Didelphis_virginiana, Felis_catus, Canis_lupus, ...
Node labels:
, 13, 14, 27, 29, 19, ...
Rooted; includes branch lengths.

Change the tip names. Create a data frame with the tip labels and new tip names.

require(phylotools)
timeTreeTips <- tree14$tip.label
replaceTips <- c("Alligator", "Cat", "Chicken", "Chimpanzee", "Cow", "Dog", "Human", "Lizard", "Macaque", "Mouse", "Opossum", "Pig", "Rabbit", "Rat")
myDat <- data.frame(timeTreeTips,replaceTips)
ntree14<- sub.taxa.label(tree14,myDat)

Collect and write the tip names to a text file

#Extract tips from newick file, write to text file
require(ape)
my.tips <- sort(tree14$tip.label)
#option 1
cat(my.tips,file="outfile.txt",sep="\n")
#option 2
my_conn = file("outfile.txt")
writeLines(my.tips,my_conn)
close(my_conn)

Next, make the plot.

plot(ntree14)

Result, a simple phylogram, i.e., a tree diagram with branching patterns and branch lengths proportional to amount of character change.

Figure 1. Phylogram plot of 14 taxa

Or, change from default “phylogram” to “cladogram” view.

plot(tree14, type="cladogram")

Figure 2.  Cladogram view, same 14 taxa.

Note that while the tree is rooted, it’s a midpoint rooting, the default setting in Newick files. For true root based on outgroup(s), identify the nodes, then select root.

Add node labels; plot() must be run first.

nodelabels()

Figure 3. Plot of tree with labeled nodes.

The outgroup(s) were the reptiles (Alligator, Chicken, Lizard), so reroot at node 16.

rrTree<- root(tree14, node=16)
plot(rrTree, type="cladogram")

Figure 4. Re-rooted tree.

Write tree to a file

require(ape)

To export tree to Newick format

write.tree(tree14, file = "filename.nwk")

for Nexus format

write.nexus(tree14, file = "filename.nex")

Star phylogeny

Collapse the tree to a star phylogeny, an unlikely evolutionary model in which the species resulted from “… a single explosive adaptive radiation” (Felsenstein 1985). Star phylogeny is an extreme tree shape, or multifurcation (polytomy), where all tips derive from the same node (Colijn and Plazzotta 2018). This type of phylogeny can be viewed as a null model for inference (but see Bayesian “star phylogeny paradox,” cf. Kolaczkowski and Thornton 2006).

require(phytools)
ctree14 <-collapse.to.star(tree14, 15)
plot(ctree14, type="cladogram")

Figure 5. Star phylogeny

Under a star phylogeny model, all taxa are assumed independent of each other, in contrast to the nested hierarchical model of evolution (e.g., Fig. 4), which shows a lack of independence among the taxa. More succinctly, a fitted to uncorrected taxa comparisons may violated the assumption that errors are independent and identically distributed. Phylogenetically correct methods attempt to address the lack of independence among taxa for comparative analysis (Felsenstein 1985, Uyeda et al 2016). Biologists should know about Felsenstein’s 1985 paper. Felsenstein’s paper created a paradigm shift in how to analyze comparative datasets and has been cited more than ten thousand times (1 August 2023, Google Scholar). To put that number in context, the 1986 paper by Kary Mullis et al., which announced invention of PCR with thermally stable polymerase that has revolutionized molecular biology, has been cited 6721 times over that same period.

Additional packages of note

The R package tanggle works with the package ggtree and advantage of the ggplot2 environment. Contains many functions to work with phylogeny graphs including re-rooting and swapping nodes. The package is available from Bioconductor,

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("tanggle")

ggtree is also a Bioconductor package, not available at CRAN.

Online viewers

Many browser-based tree viewers online, including icytree.org and iTOL tools. Additional tree viewers listed at Wikipedia.

References and suggested readings

Colijn C, Plazzotta G. A Metric on Phylogenetic Tree Shapes. Syst Biol. 2018 Jan 1;67(1):113-126.

Felsenstein, J (1985) Phylogenies and the comparative method. American Naturalist 125(1):1-15.

Felsenstein, J. (2004). Inferring phylogenies. Sunderland, MA: Sinauer associates.

Gregory, T. R. (2008). Understanding evolutionary trees. Evolution: Education and Outreach, 1(2), 121-137.

Holder, M., & Lewis, P. O. (2003). Phylogeny estimation: traditional and Bayesian approaches. Nature reviews genetics, 4(4), 275-284.

Kuhn, T. S. (1962). The structure of scientific revolutions. University of Chicago Press: Chicago.

Mullis, K., Faloona, F., Scharf, S., Saiki, R. K., Horn, G. T., & Erlich, H. (1986, January). Specific enzymatic amplification of DNA in vitro: the polymerase chain reaction. In Cold Spring Harbor symposia on quantitative biology (Vol. 51, pp. 263-273). Cold Spring Harbor Laboratory Press.

Naughton, J. (2012, August 18). Thomas Kuhn: The man who changed the way the world looked at science. The Guardian. https://www.theguardian.com/science/2012/aug/19/thomas-kuhn-structure-scientific-revolutions

O’Hara, R. J. (1997). Population thinking and tree thinking in systematics. Zoologica scripta, 26(4), 323-329.

Paradis, E. (2012) Analysis of Phylogenetics and Evolution with R (Second Edition). New York: Springer.

Paradis, E. and Schliep, K. (2019) ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics, 35, 526–528.

Revell, L. J. (2012) phytools: An R package for phylogenetic comparative biology (and other things). Methods in Ecology and Evolution, 3, 217-223.

Sandvik, H. (2008). Tree thinking cannot taken for granted: challenges for teaching phylogenetics. Theory in Biosciences, 127(1), 45-51.

Uyeda, J. C., Zenil-Ferguson, R., & Pennell, M. W. (2018). Rethinking phylogenetic comparative methods. Systematic Biology, 67(6), 1091-1109.

Zhang, J., Pei, N., & Mi, X. (2012). phylotools: Phylogenetic tools for Eco-phylogenetics. R package version 0.1, 2.

Chapter 20 contents

• Additional topics
• Area under the curve
• Peak detection
• Baseline correction
• Surveys
• Time series
• Dimensional analysis
• Estimating population size
• Diversity indexes
• Survival analysis
• Growth equations and dose response calculations
• Plot a Newick tree
• Phylogenetically independent contrasts
• How to get the distances from a distance tree
• Binary classification
• Meta-analysis

20.10 – Growth equations and dose response calculations

Introduction

In biology, growth may refer to increase in cell number, change in size of an individual across development, or increase of number of individuals in a population over time. Nonlinear, time-series, several models proposed to fit growth data, including the Gompertz, logistic, and the von Bertalanffy. These models fit many S-shaped growth curves. These models are special cases of generalized linear models, also called Richard curves.

Growth example

This page describes how to use R to analyze growth curve data sets.

Key to parameter estimates: y0 is the lag, mumax is the growth rate, and K is the asymptotic stationary growth phase. The spline function does not return an estimate for K.

R code

require(growthrates)
#Enter the data. Replace these example values with your own 
#time variable (Hours)
Hours <- c(0.000, 0.274, 0.384, 0.658, 0.986, 1.260, 1.479, 1.699, 1.973, 2.137, 2.466, 2.795, 3.123, 3.671, 4.110, 4.548, 5.151, 5.534, 5.863, 6.521, 7.068, 7.671, 8.164, 8.877, 9.644, 10.466, 11.233, 11.890, 12.493, 13.205, 14.082)
#absorbance or concentration variable (Abs)
Abs <- c(0.002207, 0.010443, 0.033688, 0.063257, 0.111848, 0.249240, 0.416236, 0.515578, 0.572632, 0.589528, 0.619091, 0.608486, 0.621136, 0.616850, 0.614689, 0.614643, 0.612465, 0.606082, 0.603933, 0.595407, 0.589006, 0.578372, 0.567749, 0.559217, 0.546451, 0.537907, 0.537826, 0.529300, 0.516551, 0.505905, 0.491013)

#Make a dataframe and check the data; If error, then check that variables have equal numbers of entries
Yeast <- data.frame(Hours,Abs); Yeast

#Obtain growth parameters from fit of a parametric growth model
#First, try some reasonable starting values
p <- c(y0 = 0.001, mumax = 0.5, K = 0.6)
model.par <- fit_growthmodel(FUN = grow_logistic, p = p, Hours, Abs, method=c("L-BFGS-B"))
summary(model.par)
coef(model.par)

#Obtain growth parameters from fit of a nonparametric smoothing spline
model.npar <- fit_spline(Hours,Abs)
summary(model.npar)
coef(spline.md)

#Make plots
par(mfrow = c(2, 1))
plot(Yeast, ylim=c(0,1), cex=1.5,pch=16, main="Parametric Nonlinear Growth Model", xlab="Hours", ylab="Abs")
lines(model.par, col="blue", lwd=2)
plot(model.npar, ylim=c(0,1), lwd=2, main="Nonparametric Spline Fit", xlab="Hours", ylab="Abs")

Results from example code

Parameters:
      Estimate Std. Error t value Pr(>|t|)    
y0    0.003765   0.001049    3.59  0.00125 ** 
mumax 3.954458   0.238810   16.56 5.37e-16 ***
K     0.579354   0.008092   71.60  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.03833 on 28 degrees of freedom

Parameter correlation:
            y0   mumax        K
y0     1.00000 -0.8755  0.05847
mumax -0.87550  1.0000 -0.14623
K      0.05847 -0.1462  1.00000
         y0       mumax           K 
0.003764518 3.954457995 0.579354378 
Fitted smoothing spline:
Call:
smooth.spline(x = time, y = ylog)

Smoothing Parameter  spar= 0.4409789  lambda= 1.221762e-05 (12 iterations)
Equivalent Degrees of Freedom (Df): 14.70929
Penalized Criterion (RSS): 0.5024203
GCV: 0.05868796

Parameter values of exponential growth curve:
Maximum growth at x= 0.0004255438 , y= 0.003072228 
y0 = 0.003065586 
mumax = 5.085845 

r2 of log transformed data= 0.9906254 
         y0       mumax 
0.003065586 5.085845036

Figure 1. Top: Parametric Nonlinear Growth Model; Bottom: Nonparametric Spline Fit

LD₅₀

In toxicology, the dose of a pathogen, radiation, or toxin required to kill half the members of a tested population of animals or cells is called the lethal dose, 50%, or LD₅₀. This measure is also known as the lethal concentration, LC₅₀, or properly after a specified test duration, the LCt₅₀ indicating the lethal concentration and time of exposure. LD₅₀ figures are frequently used as a general indicator of a substance’s acute toxicity. A lower LD₅₀ is indicative of increased toxicity.

The point at which 50% response of studied organisms to range of doses of a substance (e.g., agonist, antagonist, inhibitor, etc.) to any response, from change in behavior or life history characteristics up to and including death can be described by the methods described in this chapter. The procedures outlined below assume that there is but one inflection point, i.e., an “s-shaped” curve, either up or down; if there are more than one inflection points, then the logistic equations described will not fit the data well and other choices need to be made (see Di Veroli et al 2015). We will use the drc package (Ritz et al 2015).

Example

First we’ll work through use of R. We’ll follow up with how to use Solver in Microsoft Excel.

After starting R, load the drc library.

library(drc)

Consider some hypothetical 24-hour survival data for yeast exposed to salt solutions. Let resp equal the variable for frequency of survival (e.g., calculated from OD₆₆₀ readings) and NaCl equal the millimolar (mm) salt concentrations or doses.

At the R prompt type

resp <- c(1,1,1,.9,.7,.3,.4,.2,0,0,0)
NaCl=seq(0,1000,100) 
#Confirm sequence was correctly created; alternatively, enter the values.
NaCl 
[1] 0 100 200 300 400 500 600 700 800 900 1000 
#Make a plot
plot(NaCl,resp,pch=19,cex=1.2,col="blue",xlab="NaCl [mm]",ylab="Survival frequency")

And here is the plot (Fig. 2).

Figure 2. Hypothetical data set, survival of yeast in different salt concentrations.

Note the sigmoidal “S” shape — we’ll need an logistic equation to describe the relationship between survival of yeast and NaCl doses.

The equation for the four parameter logistic curve, also called the Hill-Slope model, is

where c is the parameter for the lower limit of the response, d is the parameter for the upper limit of the response, e is the relative EC₅₀, the dose fitted half-way between the limits c and d, and b is the relative slope around the EC₅₀. The slope, b, is also known as the Hill slope. Because this experiment included a dose of zero, a three parameter logistic curve would be appropriate. The equation simplifies to

EC₅₀ from 4 parameter model

Let’s first make a data frame

dose <- data.frame(NaCl, resp)

Then call up a function, drm, from the drc library and specify the model as the four parameter logistic equation, specified as LL.4(). We follow with a call to the summary command to retrieve output from the drm function. Note that the four-parameter logistic equation

model.dose1 = drm(dose,fct=LL.4())
summary(model.dose1)

And here is the R output.

Model fitted: Log-logistic (ED50 as parameter) (4 parms) 
Parameter estimates: 
                   Estimate  Std. Error    t-value  p-value 
b:(Intercept)      3.753415    1.074050   3.494636   0.0101 
c:(Intercept)     -0.084487    0.127962  -0.660251   0.5302 
d:(Intercept)      1.017592    0.052460  19.397441   0.0000 
e:(Intercept)    492.645128   47.679765  10.332373   0.0000 

Residual standard error: 
0.0845254 (7 degrees of freedom)

The EC₅₀, technically the LD₅₀ because the data were for survival, is the value of e: 492.65 mM NaCl.

You should always plot the predicted line from your model against the real data and inspect the fit.

At the R prompt type

plot(model.dose1, log="",pch=19,cex=1.2,col="blue",xlab="NaCl [mm]",ylab="Survival frequency")

As long as the plot you made in earlier steps is still available, R will add the line specified in the lines command. Here is the plot with the predicted logistic line displayed (Fig. 3).

Figure 3. Logistic curve added to Figure 1 plot.

While there are additional steps we can take to decide is the fit of the logistic curve was good to the data, visual inspection suggests that indeed the curve fits the data reasonably well.

More work to do

Because the EC₅₀ calculations are estimates, we should also obtain confidence intervals. The drc library provides a function called ED which will accomplish this. We can also ask what the survival was at 10% and 90% in addition to 50%, along with the confidence intervals for each.

At the R prompt type

ED(model.dose1,c(10,50,90), interval="delta")

And the output is shown below.

Estimated effective doses
(Delta method-based confidence interval(s))
     Estimate Std. Error   Lower   Upper
1:10  274.348     38.291 183.803  364.89
1:50  492.645     47.680 379.900  605.39
1:90  884.642    208.171 392.395 1376.89

Thus, the 95% confidence interval for the EC₅₀ calculated from the four parameter logistic curve was between the lower limit of 379.9 and an upper limit of 605.39 mm NaCl.

EC₅₀ from 3 parameter model

Looking at the summary output from the four parameter logistic function we see that the value for c was -0.085 and p-value was 0.53, which suggests that the lower limit was not statistically different from zero. We would expect this given that the experiment had included a control of zero mm added salt. Thus, we can explore by how much the EC₅₀ estimate changes when the additional parameter c is no longer estimated by calculating a three parameter model with LL.3().

model.dose2 = drm(dose,fct=LL.3())
summary(model.dose2)

R output follows.

Model fitted: Log-logistic (ED50 as parameter) with lower limit at 0 (3 parms)
Parameter estimates:
               Estimate Std. Error   t-value p-value
b:(Intercept)   4.46194    0.76880   5.80378   4e-04
d:(Intercept)   1.00982    0.04866  20.75272   0e+00
e:(Intercept) 467.87842   25.24633  18.53253   0e+00
Residual standard error:
 0.08267671 (8 degrees of freedom)

The EC₅₀ is the value of e: 467.88 mM NaCl.

How do the four and three parameter models compare? We can rephrase this as as statistical test of fit; which model fits the data better, a three parameter or a four parameter model?

At the R prompt type

anova(model.dose1, model.dose2)

The R output follows.

1st mode
fct:      LL.3()
2nd model
fct:      LL.4()
ANOVA table
          ModelDf      RSS Df F value p value
2nd model       8 0.054684                   
1st model       7 0.050012  1  0.6539  0.4453

Because the p-value is much greater than 5% we may conclude that the fit of the four parameter model was not significantly better than the fit of the three parameter model. Thus, based on our criteria we established in discussions of model fit in Chapter  16 – 18, we would conclude that the three parameter model is the preferred model.

The plot below now includes the fit of the four parameter model (red line) and the three parameter model (green line) to the data (Fig. 4).

Figure 4. Four parameter (red) and three parameter (green) logistic models fitted to data.

The R command to make this addition to our active plot was

lines(dose,predict(model.dose2, data.frame(x=dose)),col="green")

We continue with our analysis of the three parameter model and produce the confidence intervals for the EC₅₀ (modify the ED() statement above for model.dose2 in place of model.dose1).

Estimated effective doses
(Delta method-based confidence interval(s))
     Estimate Std. Error   Lower  Upper
1:10  285.937     33.154 209.483 362.39
1:50  467.878     25.246 409.660 526.10
1:90  765.589     63.026 620.251 910.93

Thus, the 95% confidence interval for the EC₅₀ calculated from the three parameter logistic curve was between the lower limit of 409.7 and an upper limit of 526.1 mm NaCl. The difference between upper and lower limits was 116.4 mm NaCl, a smaller difference than the interval calculated for the 95% confidence intervals from the four parameter model (225.5 mm NaCl). This demonstrates the estimation trade-off: more parameters to estimate reduces the confidence in any one parameter estimate.

Additional notes of EC₅₀ calculations

Care must be taken that the model fits the data well. What if we did not have observations throughout the range of the sigmoidal shape? We can explore this by taking a subset of the data.

dd = dose[1:6,]

Here, all values after dose 500 were dropped

dd
  resp dose
1  1.0    0
2  1.0  100
3  1.0  200
4  0.9  300
5  0.7  400
6  0.3  500

and the plot does not show an obvious sigmoidal shape (Fig. 5).

Figure 5. Plot of reduced data set.

We run the three parameter model again, this time on the subset of the data.

model.dosedd = drm(dd,fct=LL.3())
summary(model.dosedd)

Output from the results are

Model fitted: Log-logistic (ED50 as parameter) with lower limit at 0 (3 parms)
Parameter estimates:
                Estimate Std. Error    t-value p-value
b:(Intercept)   6.989842   0.760112   9.195801  0.0027
d:(Intercept)   0.993391   0.014793  67.153883  0.0000
e:(Intercept) 446.882542   5.905728  75.669344  0.0000
Residual standard error:
 0.02574154 (3 degrees of freedom)

Conclusion? The estimate is different, but only just so, 447 vs. 468 mm NaCl.

Thus, within reason, the drc function performs well for the calculation of EC₅₀. Not all tools available to the student will do as well.

NLopt and nloptr

draft

Free open source library for nonlinear optimization.

Steven G. Johnson, The NLopt nonlinear-optimization package, http://github.com/stevengj/nlopt

https://cran.r-project.org/web/packages/nloptr/vignettes/nloptr.html

Alternatives to R

What about online tools? There are several online sites that will allow students to perform these kinds of calculations. Students familiar with MatLab know that it can be used to solve nonlinear equation problems. An open source alternative to MatLab is called GNU Octave, which can be installed on a personal computer or run online at http://octave-online.net. Students also may be aware of other sites, e.g., mycurvefit.com and IC50.tk.

Both of these free services performed well on the full dataset (results not shown), but fared poorly on the reduced subset: mycurvefit returned a value of 747.64 and IC50.tk returned an EC₅₀ estimate of 103.6 (Michael Dohm, pers. obs.).

Simple inspection of the plotted values shows that these values are unreasonable.

EC₅₀ calculations with Microsoft Excel

Most versions of Microsoft Excel include an add-in called Solver, which will permit mathematical modeling. The add-in is not installed as part of the default installation, but can be installed via the Options tab in the File menu for a local installation or via Insert for the Microsoft Office online version of Excel (you need a free Microsoft account). The following instructions are for a local installation of Microsoft Office 365 and were modified from information provided by sharpstatistics.co.uk.

After opening Excel, set up worksheet as follows. Note that in order to use my formulas your spreadsheet values need to be set up exactly as I describe.

• Dose values in column A, beginning with row 2
• Response values in column B, beginning with row 2
• In column F type b in row 2, c in row 3, d in row 4, and e in row 5.
• In cells G2 – G5, enter the starting values. For this example, I used b = 1, c = 0, d = 1, e = 400.
  • For your own data you will have to explore use of different starting values.
• Enter headers in row 1.
  • Column A Dose
  • Column B Response
  • Column C Predicted
  • Column D Squared difference
  • Column F Constants
• In cell C14 enter sum squares
• In cell D14 enter =sum(D2:D12)

Here is an image of the worksheet (Fig. 6), with equations entered and values updated, but before Solver run.

Figure 6. Screenshot Microsoft Excel worksheet containing our data set (col A & B), with formulas added and calculated. Starting values for constants in column G, rows 2 – 4.

Next, enter the functions.

• In cell C2 enter the four parameter logistic formula — type everything between the double quotes: “=$G$3+(($G$4-$G$3)/(1+(A2/$G$5)^$G$2))“. Next, copy or drag the formula to cell C12 to complete the predictions.
  • Note: for a three parameter model, replace the above formula with “=(($G$4)/(1+(A2/$G$5)^$G$2))“.
• In cell D2 type everything between the double quotes: “=(B2-C2)^2“.
• Next, copy or drag the formula to cell D12 to complete the predictions.
• Now you are ready to run solver to estimate the values for b, c, d, and e.
  • Reminder: starting values must be in cells G2:G5
• Select cell D14 and start solver by clicking on Data and looking to the right in the ribbon.

Solver is not normally installed as part of the default installation of Microsoft Excel 365. If the add-in has been installed you will see Solver in the Analyze box (Fig. 7).

Figure 7. Screenshot Microsoft Excel, Solver add-in available.

If Solver is not installed, go to File → Options → Add-ins (Fig. 8).

Figure 8. Screenshot Microsoft Excel, Solver add-in available and ready for use.

Go to Microsoft support for assistance with add-ins.

With the spreadsheet completed and Solver available, click on Solver in the ribbon (Fig. 7) to begin. The screen shot from the first screen of Solver is shown below (Fig. 9).

Figure 9. Screenshot Microsoft Excel solver menu.

Setup Solver

• Set Objective, enter the absolute cell reference to the sum squares value, \$D$14
• Set To: Min.
• For By Changing Variable Cells, enter the range of cells for the four parameters, $G$2:$G$5
• Uncheck the box by “Make Unconstrained Variables Non-Negative.”
  • Where constraints refers to any system of equalities or inequalities equations imposed on the algorithm.
• Select a Solving method, choose “GRG Nonlinear,” the nonlinear programming solver option (not shown in Fig. 7, select by clicking on the down arrow).
• Click Solve button to proceed.

If all goes well, this next screen (Fig. 10) will appear which shows the message, “Solver has converged to the current solution. All constraints are satisfied.”

Figure 10. Screenshot solver completed run.

Click OK and note that the values of the parameters have been updated (Table 1).

Table 1. Four parameter logistic model, results from Solver 

Note the values obtained by Solver are virtually identical to the values obtained in the drc R package. The differences are probably because of the solver algorithm.

More interestingly, how did Solver do on the subset data set? Here are the results from a three parameter logistic model (Table 2).

Table 2. Three parameter logistic model, results from Solver

The results are again very close to results from the drc R package.

Thus, we would conclude that Solver and Microsoft Excel would be a reasonable choice for EC₅₀ calculations, and much better than IC50.tk and mycurvefit.com. The advantages of R over Microsoft Excel is that the model building is more straightforward than entering formulas in the cell reference format required by Excel.

Questions

[pending]

References and suggested readings

Beck B, Chen YF, Dere W, et al. Assay Operations for SAR Support. 2012 May 1 [Updated 2012 Oct 1]. In: Sittampalam GS, Coussens NP, Nelson H, et al., editors. Assay Guidance Manual [Internet]. Bethesda (MD): Eli Lilly & Company and the National Center for Advancing Translational Sciences; 2004-. Available from: www.ncbi.nlm.nih.gov/books/NBK91994/

Di Veroli G. Y., Fornari C., Goldlust I., Mills G., Koh S. B., Bramhall J. L., Richards, F. M., Jodrell D. I. (2015) An automated fitting procedure and software for dose-response curves with multiphasic features. Scientific Reports 5: 14701.

Ritz, C., Baty, F., Streibig, J. C., Gerhard, D. (2015) Dose-Response Analysis Using R PLOS ONE, 10(12), e0146021

Tjørve, K. M., & Tjørve, E. (2017). The use of Gompertz models in growth analyses, and new Gompertz-model approach: An addition to the Unified-Richards family. PloS one, 12(6), e0178691.

Chapter 20 contents

• Additional topics
• Area under the curve
• Peak detection
• Baseline correction
• Surveys
• Time series
• Dimensional analysis
• Estimating population size
• Diversity indexes
• Survival analysis
• Growth equations and dose response calculations
• Plot a Newick tree
• Phylogenetically independent contrasts
• How to get the distances from a distance tree
• Binary classification
• Meta-analysis

/MD

20.9 – Survival analysis

draft

Introduction

As the name suggests, survival analysis is a branch of statistics used to account for death of organisms, or more generally, failure in a system.  In general, this kind of analysis models time to event, where event would be death or failure.  The basics of the method is defined by the survival function

where is time, is a variable that represents time of death or other end point, and is probability of an event occurring later than at time . Two methods are common: the descriptive Kaplan-Meier estimator and the the semi-parametric Cox proportional hazards model.

Excellent resources available, series of articles in volume 89 of British Journal of Cancer: Clark et al (2003a), Bradburn et al (2003a), Bradburn et al (2003b), and Clark et al (2003b).

Hazard function

Defined as the event rate at time t based on survival for time times equal to or greater than t.

Censoring

Censoring is a a missing data problem typical of survival analysis.  Distinguish right-censored and left-censored.

Kaplan-Meier plot

Kaplan-Meier (KM) estimator of survival function. Other survival function estimators Fleming-Harrington. The KM estimator, is

where d_i is the number of events that occurred at time t_i, n_i is the number of individuals known to have survived or not been censored. Because it’s an estimate, a statistic, need an estimate of the error variance. Several options, the default in R is the Greenwood estimator.

The KM plot, censoring times noted with plus.

R code

Download and install the RcmdrPlugin.survival package.

Example

data(heart, package="survival")
attach(heart)
#Get help with the data set
help("heart", package="survival")

head(heart)

 start stop event        age      year surgery transplant id
1    0   50     1 -17.155373 0.1232033       0          0  1
2    0    6     1   3.835729 0.2546201       0          0  2
3    0    1     0   6.297057 0.2655715       0          0  3
4    1   16     1   6.297057 0.2655715       0          1  3
5    0   36     0  -7.737166 0.4900753       0          0  4
6   36   39     1  -7.737166 0.4900753       0          1  4

Run basic survival analysis. After installing the RcmdrPlugin.survival, from Rcmdr select estimate survival function.

Figure 1. Screenshot of menu call for survival analysis in Rcmdr

Get survival estimator and KM plot (Figure 2)

R output

.Survfit <- survfit(Surv(start, event) ~ 1, conf.type="log", conf.int=0.95, 
Rcmdr+ type="kaplan-meier", error="greenwood", data=heart)

 .Survfit
Call: survfit(formula = Surv(start, event) ~ 1, data = heart, error = "greenwood", 
conf.type = "log", conf.int = 0.95, type = "kaplan-meier")

n events median 0.95LCL 0.95UCL 
172 75 26 17 37

plot(.Survfit, mark.time=TRUE) 

quantile(.Survfit, quantiles=c(.25,.5,.75))
$quantile
25 50 75 
3 26 67

$lower
25 50 75 
1 17 46

$upper
25 50 75 
12 37 NA

#by default, Rcmdr removes the object
remove(.Survfit)

Figure 2. Kaplan-Meier plot of heart data. Dashed lines are upper and lower confidence intervals about the survival function.

Interpreting the plot is straightforward — every time an event occurs, the curve drops.

ylim = c(0,1), ylab="Survival probability", xlab="Days", lwd=2, col="blue"

the data set includes age and whether or not subjects had heart surgery before transplant. Compare.

Variable surgery is recorded 0,1, so need to create a factor

fSurgery <- as.factor(surgery)

Now, to compare

Rcmdr: Statistics → Survival analysis → Compare survival functions…

R output

Rcmdr> survdiff(Surv(start,event) ~ fSurgery, rho=0, data=heart)
Call:
survdiff(formula = Surv(start, event) ~ fSurgery, data = heart, 
rho = 0)

              N Observed Expected (O-E)^2/E (O-E)^2/V
fSurgery=No 143       66     58.7     0.902      4.56
fSurgery=Yes 29        9     16.3     3.255      4.56

Chisq= 4.6 on 1 degrees of freedom, p= 0.03

Get the KM estimator and make a KM plot

mySurvfit <- survfit(Surv(start, event) ~ surgery, conf.type="log", conf.int=0.95,
type="kaplan-meier", error="greenwood", data=heart)

plot(mySurvfit, mark.time=TRUE, ylim=c(0,1),lwd=2, col=c("blue","red"), xlab="Number of days", ylab="Survival probability")
legend("topright", legend = paste(c("Surgery - No", "Surgery - Yes")), col = c("blue", "red"), pch = 19, bty = "n")

Figure 3. KM plot comparing survival of patients who received surgery (red) and those who did not have the surgery (blue).

The comparison plot can be made in Rcmdr by selecting our Surgery factor in Strata setting (Fig. 4). Recall that strata refers to subgroups of a population.

Figure 4. Screenshot of Survival estimator menu in Rcmdr.

Cox proportional hazards model

coxph()

Questions

[pending]

Chapter 20 contents

• Additional topics
• Area under the curve
• Peak detection
• Baseline correction
• Surveys
• Time series
• Dimensional analysis
• Estimating population size
• Diversity indexes
• Survival analysis
• Growth equations and dose response calculations
• Plot a Newick tree
• Phylogenetically independent contrasts
• How to get the distances from a distance tree
• Binary classification
• Meta-analysis

/MD

20.1 – Area under the curve

Introduction.

The classical physiological area under the curve, AUC, represents the total change in given change in . For example, if is time, and is oxygen consumption — eg, respirometry: recorded as an analog voltage output (millivolts or volts) from a gas analyzer to a data acquisition (DAQ) board — an AUC would be appropriate to quantify the total oxygen consumption following strenuous exercise (Excess post-exercise oxygen consumption, EPOC) or following a large meal (Specific Dynamic Action, SDA).  The area represents the total, cumulative volume of oxygen consumed over that specific timeframe, eg, calculating an individual’s total EPOC.

In machince learning, the axes are probabilities and error rates. The area tells you how accurate a decision-making model is. We plot True Positive Rate (Sensitivity) against False Positive Rate (1 – Specificity) across every possible mathematical decision threshold (from 0% to 100% confidence). The curve now represents the shifting balance between catching true positives and committing false positives and the area summarizes how cleanly the model separates two distinct classes, eg, Sick vs. healthy. An area of 1.0 means perfect separation with zero overlap.

Area under the relative (receiver) operating carrier, AUROC, shows characteristics of a diagnostic model, a graphic used to show trade off between sensitivity and specificity. Classifier performance. Used to find the appropriate cut-off. Plot true positive rates against false positive rates as cumulative functions, shows the relationship between sensitivity and specificity for every possible cut off value. Can then calculate AUC to get a measure of the intervention’s ability to discriminate between true and false positive rates.

Related, area under precision-recall curve, AUPRC, is a single scalar metric that measures a binary classification model’s ability to balance precision (how many predicted positives are actually positive) and recall (how many actual positives are caught) across all possible decision thresholds. Here, precision and recall refer to both the accuracy of the model’s positive predictions and its ability to find all the positive cases.

From Chapter 7 we defined

where is true positives and is false positives.

We now introduce recall, which is also refered to as sensitivity

Precision measures the quality of your model’s positive guesses; Recall measures how many actual positive cases in the dataset were identified by the model.

AUC evaluates a model’s ability to distinguish between both and equally, while AUPRC focuses strictly on the model’s performance regarding the positive (rare) class.

AUC would be preferred when your dataset is balanced (equal positives and negatives) OR when you care equally about performance on both classes. For example, dataset of patients with flu or Covid from a clinic during winter months where there is roughly balanced 50/50 split of patients presenting with either Flu or Covid. You want a model that accurately classifies both which favors AUC approach.

AUPRC would be the choice when your dataset is highly imbalanced, ie, the positive class is rare. In contrast, AUPRC would be favored for scenarios where a researcher is screening 100,000 asymptomatic patients for a rare disease that only 10 people actually have.

Physiology Area under the curve

Download and install R package MESS; requires geepack, geeM, and Matrix packages

R code

x <- seq(1:10) 
y <- c(1,4,5,2,11,22,9,7,5,1) 
#length(x)==length(y)
#smooth the data 
loxy <- loess(y~x)
#Make a plot (Fig. 1)
plot(x,y, pch=19, cex=2, col="blue") 
lines(predict(loxy), type="l", col="red")

where == is an R comparison operator.

And R output

Figure 1. Area under the curve example.

library(MESS) 
auc(x,y,from=0,rule=2) 
auc(x,loxy$fitted,from=0,rule=2)

And R output

#area under curve for raw data
[1] 67
#area under curve for smoothed data
[1] 66.77616

Machine Learning area under the curve.

# Install and load the ROCR package if you haven't already
# install.packages("ROCR")
# modified from https://rviews.rstudio.com/2019/03/01/some-r-packages-for-roc-curves/

library(ROCR)
df <- data.frame(ROCR.simple)
pred <- prediction(df$predictions, df$labels)
perf <- performance(pred,"tpr","fpr")

# 1. Calculate the actual numeric AUC value
auc_perf <- performance(pred, "auc")
auc_value <- round(unlist(auc_perf@y.values), 4)

# 2. Plot the ROC curve with colorized thresholds and custom labels
plot(perf, 
colorize = TRUE, 
main = "ROC Curve",
xlab = "False Positive Rate (1 - Specificity)",
ylab = "True Positive Rate (Sensitivity)")

# 3. Add the AUC value as a subtitle
mtext(paste("AUC =", auc_value), side = 3, line = 0.5)

# 4. Add a diagonal 50/50 reference line
abline(a = 0, b = 1, lty = 2, col = "darkgray")

End R code.

R output.

The calculated AUC was . The right-hand axes are color codes by AUC values (Fig 2): good tests AUC between 0.8 and 0.9, very good tests greater than 0.9.

Figure 2. Example ROC curve.

Area under the precision recall curve

AUPRC measures the overall performance of a binary classification model by summarizing how well it balances precision and recall across all possible decision thresholds. A score closer to 1.0 means the model is very good at finding real positive cases without making many mistakes, which is especially useful when your dataset has very few positive examples. A score closer to 1.0 means the model is very good at finding real positive cases without making many mistakes.

R code.

library(ROCR)

# 1. Prepare data and prediction object
df <- data.frame(ROCR.simple)
pred <- prediction(df$predictions, df$labels)

# 2. Extract Precision and Recall for the plot
perf_pr <- performance(pred, measure = "prec", x.measure = "rec")

# 3. Calculate the numeric AUPRC value
# In ROCR, AUPRC is calculated by finding the Area Under the Precision/Recall curve ("aucpr")
perf_auprc <- performance(pred, measure = "aucpr")
auprc_value <- round(unlist(perf_auprc@y.values), 4)

# Print the final AUPRC score to the console
cat("The calculated AUPRC is:", auprc_value, "\n")

# 4. Plot the Precision-Recall Curve using base R graphics
plot(perf_pr, 
colorize = TRUE, 
main = paste("Precision-Recall Curve (AUPRC =", auprc_value, ")"),
xlab = "Recall (Sensitivity)", 
ylab = "Precision (Positive Predictive Value)")

End R code.

R output.

The calculated AUPRC was .

Figure 3. Precision recall curve.

The differences: AUC evaluates how well the model separates both classes equally, while AUPRC evaluates how reliably the model handles the positive class alone.

The Interpretive baselines:

• Figure 2: AUC Baseline is always 0.5. A completely random model yields a diagonal line from (0,0) to (1,1), resulting in an AUC of 0.5. This baseline never changes, regardless of how rare your disease is.
• Figure 3: AUPRC Baseline depends on prevalence. A random model yields a flat horizontal line across the plot. The height of this line equals the proportion of positive cases in your dataset (). If only 1% of the patients have a rare condition, a random model’s AUPRC baseline is 0.01.

What a shifting threshold means:

• On an AUC Curve: As you lower your decision threshold to catch more sick people, you track how fast your True Positive Rate rises compared to your False Positive Rate. You are looking at the overall trade-off between the entire healthy population and the entire sick population.
• On an AUPRC Curve: As you lower your threshold to catch more sick people, you track Recall against Precision. This evaluates the quality of the “Positive” alerts. It directly measures how quickly the pool of patients you flagged as “high risk” becomes diluted with false alarms.

References and suggested readings

Çorbacıoğlu, Ş. K., & Aksel, G. (2023). Receiver operating characteristic curve analysis in diagnostic accuracy studies: A guide to interpreting the area under the curve value. Turkish Journal of Emergency Medicine, 23(4), 195–198.

Lobo, J. M., Jiménez-Valverde, A., & Real, R. (2008). AUC: A misleading measure of the performance of predictive distribution models. Global Ecology and Biogeography, 17(2), 145–151.

Martin, B., Bennett, T. D., DeWitt, P. E., Russell, S., & Sanchez-Pinto, L. N. (2025). Use of the Area Under the Precision-Recall Curve to Evaluate Prediction Models of Rare Critical Illness Events. Pediatric Critical Care Medicine, 26(6), e855–e859.

Ozenne, B., Subtil, F., & Maucort-Boulch, D. (2015). The precision–recall curve overcame the optimism of the receiver operating characteristic curve in rare diseases. Journal of Clinical Epidemiology, 68(8), 855–859. https://doi.org/10.1016/j.jclinepi.2015.02.010

Steyerberg, E. W., & Vergouwe, Y. (2014). Towards better clinical prediction models: Seven steps for development and an ABCD for validation. European Heart Journal, 35(29), 1925–1931.

Questions

1. Write up three learning outcomes for this page. Hint: Point your favorite generative AI to this page and ask for help

Chapter 20 contents

• Additional topics
• Area under the curve
• Peak detection
• Baseline correction
• Surveys
• Time series
• Dimensional analysis
• Estimating population size
• Diversity indexes
• Survival analysis
• Growth equations and dose response calculations
• Plot a Newick tree
• Phylogenetically independent contrasts
• How to get the distances from a distance tree
• Binary classification
• Meta-analysis

/MD

6.6 – Continuous distributions

Large numbers and the Central Limit Theorem

Imagine we’ve collected (sampled) data from a population and now want to summarize the data sample. How do we proceed? A good starting point is to plot the data in a histogram and note the shape of the sample distribution. Not to get too far ahead of ourselves here, but much of the classical inferential statistics demands that we are able to assume that the sampled values come from a certain kind of distribution called the normal, or Gaussian distribution.

Consider a random sample drawn from a normally distributed population of the following series of graphs, Figures 1 – 4:

Figure 1. Sample size = 20, drawn from population with known μ = 0 and σ = 1.

Figure 2. Sample size = 100, also drawn from population with known μ = 0 and σ = 1.

Figure 3. Sample size n = 1000, once again drawn from population with known μ = 0 and σ = 1.

Figure 4. And lastly, sample size n = 1 million also drawn from population with known μ = 0 and σ = 1.

These graphs illustrate a fundamental point in statistics: for many kinds of measurements in biology, the more data you sample, the more likely the data will approach a normal distribution. This series of simulations was a quick and dirty “proof” of the Central Limit Theorem, which is one of the two fundamental theorems of probability, the other being that of Law of Large Numbers (ie, large-sample statistics). Basically the CLT says that for a large number of random samples, the sample mean will approach the population mean,  μ, and the sample variance will approach the population variance σ²; the distribution of the large sample will converge on the normal distribution.

As the sample size gets bigger and bigger, the resulting sample means and standard deviations get closer and closer to the true value (remember — I TOLD the program to grab numbers from the Z distribution with a mean of zero and standard deviation = 1), obeying the Law of Large Numbers.

Simulation

I used the computer to generate sample data from a population. This process is called a simulation. R can make new data sets by sampling from known populations with specified distribution properties that we determine in advance — a very powerful tool — a technique used for many kinds of statistics (eg, Monte Carlo methods, bootstrapping, etc., see Chapter 19).

1. The first graph is for a sample of 20 points;
2. the second for 100;
3. the third for 1,000;
4. and lastly, 1 million points.

To generate a sample from a normal population, the Rcmdr call the menu by selecting:

Rcmdr: Distributions → Continuous distributions → Normal distribution → Sample from normal distribution…

Figure 5. Screenshot Rcmdr menu, sample from a normal distribution.

The menu pops up. I entered Mean (μ) = 0 and Standard deviation (σ) = 1 number of samples = 10, unchecked all boxes under the “add to data set.” I left the object name as “NormalSamples,” you can, of course, change it as needed. R code derived from these requests were

normalityTest(~obs, test="shapiro.test", data=NormalSamples)
NormalSamples <- as.data.frame(matrix(rnorm(10*1, mean=0, sd=1), ncol=1))
rownames(NormalSamples) <- paste("sample", 1:10, sep="")
colnames(NormalSamples) <- "obs"

This results in a new data.frame called NormalSamples with a single variable called obs.

for (i in 1:3){
print(rnorm(5))
}
[1] -0.4221672 -1.4317800 -1.8310352 0.4181184 -1.1596058
[1] -0.2034944 1.1809083 1.5925296 -2.0763677 1.6982357
[1] -1.0967218 -0.3205041 -1.7513838 -0.3335311 -1.8808454

for (i in 1:3){
set.seed(1) 
print(rnorm(5))
}
[1] -0.6264538 0.1836433 -0.8356286 1.5952808 0.3295078
[1] -0.6264538 0.1836433 -0.8356286 1.5952808 0.3295078
[1] -0.6264538 0.1836433 -0.8356286 1.5952808 0.3295078

After creating a sample of numbers drawn from the normal distribution, make a histogram, Rcmdr: Graphs → Histogram… (see Chapter 4.2).

The normal distribution

More on the importance of normal curves in a moment. Sometimes people call these “bell-shaped” curves. You may also see such distributions referred to as Gaussian Distributions, but the normal curve is but one of many Gaussian -type distributions. Moreover, not all “bell-shaped” curves are NORMAL. For a distribution to be “normally distributed” it must follow a specific formula.

This formula has a lot of parts, but once we break down the parts, we’ll see that there is just two things to know. First, let’s identify the parts of the equation:

 is the height of the curve (normal density)
(pie) is a constant = 3.14159… (R, use pi)
µ is the population mean
σ² is the population variance
σ is the square-root of the variance or the population standard deviation
e is the natural logarithm (R, use exp())
 is the individual’s value

Why the Normal distribution is so important in classical statistics

With these distinctions out of the way, the first important message about the normal curve is that it permits us to say how likely (ie, how probable) is a particular value if the observation comes from a population with mean µ and standard deviation σ and the population from which the sample was drawn came from a normal distribution.

The second message: all we need to know to recreate the normal distribution for a set of data is the mean and the variance (or the standard deviation) for the population!! With just these two parameters, we can then determine the expected proportion of observations expected for each value of X. Note — we generally do not know these two because they are population parameters: we must estimate them from samples, using our sample statistics, and that’s where the first big assumption in conducting statistical analyses comes into play!!

An example for calculating the normal distribution when knowing the mean and variance

µ = 5, σ² = 10 (thus, the standard deviation, σ = 3.16)

The formula becomes

Now, plug in different values of X (for example, what’s the probability that a value of X could be 0, 1, 2, …, 10? If we really do have a normal curve with mean = 5, and variance = 10).

The normal equation returns the proportion of observations in a normal population for each X value:

When i = 5 Y₅ = 0.12616 This is the proportion of all data points that have an X = 5 When i = 1 Y₁ = 0.019995 This is the proportion of all data points that have an X = 1

We can keep going, getting the proportion expected for each value of X, then make the plot (Fig 6).

Figure 6. Frequency expected for a few points (X: 0 – 10) drawn from a normal distribution, calculated using the formula and example values. 

Here’s the R code for the plot

X = seq(0,10, by=1)
Y = (0.398947/3.16)*exp((-1*(X-5)^2)/20)
plot(Y~X, ylab="Frequency (Y)", cex=1.5, pch=19,col="blue")
lines(X,Y, col="blue", lwd=2)

Next up, more about the normal distribution, Chapter 6.7.

Questions

1. For a mean of 0 and standard deviation of 1, apply the equation for the normal curve for X = (-4, -3.5, -3, -2.5, -2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4). Plot your results.
2. Sample from a normal distribution with different sample size, means, and standard deviations, each time make a histogram and compare the shape of the histograms.
3. For  a subset of

Quiz Chapter 6.6

Continuous distributions

Start Quiz

Question

✕

✓

Your answer:

Correct answer:

Next

Please wait..

↺

You got {{SCORE_CORRECT}} out of {{SCORE_TOTAL}}

Your Answers

Chapter 6 contents

• Introduction
• Some preliminaries
• Ratios and proportions
• Combinations and permutations
• Types of probability
• Discrete probability distributions
• Continuous distributions
• Normal distribution and the normal deviate
• Moments
• Chi-square distribution
• t distribution
• F distribution
• References and suggested readings

6.4 – Types of probability

Introduction

By probability, we mean a number that quantifies the uncertainty associated with a particular event or outcome. If an outcome is certain to occur, the probability is equal to 1; if an outcome cannot happen, the probability is equal to zero. Most events we talk about in biology have probability somewhere in between.

A basic understanding, or at least an appreciation for probability is important to your education in biostatistics. Simply put, there is no certainty in biology beyond simple expectations like Benjamin Franklin’s famous quip about “death and taxes…”

Discrete probability

You are probably already familiar with discrete probability. For example, what is the probability that a single toss of a fair coin will result in a “heads?” The outcomes of the coin toss are discrete, categorical, either “heads” or “tails.”

And for ten tosses of the coin? And more cogently, what is the probability that you will toss a coin ten times and get all ten “heads?” While different in tone, these are discrete outcomes. An important concept is independence. Are the multiple events independent? In other words, does the toss of coin on the first attempt affect the toss of the coin on the second attempt and so on up to the tenth toss? At least in principle the repeated tosses are independent so to find the probability you just multiply each event’s probability to get the total.

In contrast, if one or more events are not independent but somehow influence the behavior of the next event, then you add the probabilities for each dependent event. We can do better than simply multiply or add events one at a time; depending on the number of discrete outcomes it is very likely that someone has already calculated all possible outcomes and come up with an equation. In the case of the tossing of two coins this is a binomial equation problem and repeat tosses can be modeled by use of the Bernoulli distribution.

Now, try this on for size. What is the probability that the next child born at Kapiolani Medical Center in Honolulu will be assigned female?

We just described a discrete random variable, which can only take on discrete or “countable” numbers.  This distribution of values is the probability mass function. The probability of a one fatal airline accident in a year exactly on 20.1 is practically zero (the area under a point along the curve is zero), so we can get the probability of a range of values around the point as our answer.

Continuous probability

Many events in biology are of degree, not kind. It is kind of awkward to think about it, but for a sample of adult house mice drawn from a population, what is the probability of obtaining a mouse that is 20.0000 grams (g) in weight? Each possible value of body mass for a mouse is considered an event, just like in our example of tossing a coin. But clearly, we don’t expect to get a lot of mice that are exactly 20.0000 g in weight. For variables like body mass, the type of data we collect is continuous, and the probability values need to be rethought along a continuum of possible values and, in turn, how likely each value is for a mouse. Although it is theoretically possible that a mouse could weigh ten pounds, we know by experience that this is impossible. Adult mice weigh between 15 and 50 g or there about.

We just described a continuous random variable, which can take on any value between a specific interval of values. This distribution of values is the probability density function (PDF). The probability of a mouse’s weight falling exactly on 20.1 is practically zero (the area under a point along the curve is zero), so we can get the probability of a range of values around the point as our answer.

In statistical inference, following our measurements of the variables from our sample drawn from a population, we make conclusions with the following kind of caveat: “the mean body mass for this strain of mouse is 20 g.” That is our best estimate of the mean (middle) for the population of mice, more specifically, for the body mass of the mice. Here, the variable body mass is more formally termed a random variable. This implies that there is in fact a true population mean body mass for the mice and that any deviations from that mean are due to chance.

In statistics we don’t settle for a single point estimate of the population mean. You will find most reporting of estimates of random variables is accompanies by a statement like, the mean was 20g with a 95% probability that the true population mean is between 18.9 and 21 g. This is called the 95% confidence interval for the mean and it takes into account how good of an estimate our sample is likely to be relative to the true population value. Not only are we saying that we think the population mean is 20g, but we’re willing to say that we are 95% certain that the true value must be between a lower limit (18.9 g) and and upper limit (21 g). In order to make this kind of statement we have to assume a distribution that will describe the probability of mouse weights. For many reasons we usually assume a normal distribution. Once we make this assumption we can calculate how probable a particular weight is for a mouse.

We introduced how to calculate Confidence Intervals in Chapter 3.4 and will extend this in any chapter that introduces an estimation.

Types of probability

To begin refining our concept of probability, it is sometimes useful to distinguish among kinds of probabilities:

• between theoretical and empirical;
• between subjective and objective.

In most cases, including your statistics book, we would begin our discussion of probability by talking of some probabilities for events we’re familiar with.

1. The theoretical probability of heads appearing the next time you flip a fair coin is 1/2 or 50%. As long as we’re talking about a fair coin, the probability of a heads appearing each time you flip the coin remains 50%. We can check this by conducting and experiment: out of 10 tosses, how many heads appear? The answer would be an empirical probability, and we understand the chance in an objective manner (no interpretation needed).
2. The theoretical probability that a “5” will appear on the face of a fair dice after a toss is 1/6 or 16.667%. Again, as long as we’re talking about a fair dice, the probability of a “5” appearing each time you roll the dice remains 16.667%.
3. The probability that at birth, a human baby’s sex will be male about 1/2 or 50%. This is an empirical probability based on millions of observations. Changes in technology, and ethical standards notwithstanding, the probability will remain the same.
4. The probability of the birth of a Downs syndrome baby is 1/800, but increases with age until by age 45, the chance is 1/12. Again, these are empirical and objective.
5. The probability of winning the Publisher’s Clearing House Sweepstakes is about 1 in 100 million. This probability is theoretical, it is also objective; however, by adding lots of twists to the game, by having multiple opportunities and by giving the appearance that a person must purchase a magazine, some players perceive their chances as increasing or decreasing by their efforts (=subjective).

R and distributions

R offers access to many named distributions. R generally provides four functions available for each.

density (d):the probability density function (PDF) for continuous distributions or the probability mass function (PMF) for discrete distributions.

probability (p): the cumulative distribution function (CDF), which is the chance (probability) of a random variable being less than or equal to a given value

quantile (q): the quantile function (inverse CDF), which returns the value at a given probability

r: draw random numbers from the specified distribution.

Examples of functions for the normal distribution then, are: dnorm(), pnorm(), qnorm(), and rnorm().

Coding example: Find the PDF for x = 0 in a standard normal distribution. # Find the PDF value for x = 0 in a standard normal distribution (mean=0, sd=1)

dnorm(x = 0, mean = 0, sd = 1)

R output

0.398942280401433

Thus, we conclude the chance that x = 0 is about 40%.

Coding example: Find the probability that for x is equal to -1.2 or less in a standard normal distribution.

pnorm(q = -1.2, mean = 0, sd = 1, lower.tail = TRUE)

R output

0.115069670221708

Thus, we conclude the chance that x is equal to or less than -1.2 is about 11.5%.

Coding example: Assuming normal distribution of test scores, what score will 50% of students score? Find the probability that for x is equal to -1.2 or less in a standard normal distribution.

qnorm(0.5, mean = 75, sd = 10)

R output

75

As expected, at 50% it returns the mean. Now, confirm that the lower 25% and upper 75% are 68.26 and 81.75, respectively.

Last coding example, draw random sample of three scores from the same exam.

rnorm(3, mean = 75, sd = 10)

R output

1. 66.0384856946112

2. 69.1887923903669

3. 77.6342297053815

R Commander: Distributions menus give four options

• Quantiles
• Probabilities
• Plot
• Sampling

Questions

1. Define and distinguish, with examples

1. discrete and continuous probability
2. theoretical and empirical probability
3. subjective and objective probability

Quiz Chapter 6.4

Types of probability

Start Quiz

Question

✕

✓

Your answer:

Correct answer:

Next

Please wait..

↺

You got {{SCORE_CORRECT}} out of {{SCORE_TOTAL}}

Your Answers

Chapter 6 contents

• Introduction
• Some preliminaries
• Ratios and proportions
• Combinations and permutations
• Types of probability
• Discrete probability distributions
• Continuous distributions
• Normal distribution and the normal deviate
• Moments
• Chi-square distribution
• t distribution
• F distribution
• References and suggested readings

4.6 – Adding a second Y axis

Introduction

Scatter plots are used to show association between two continuous variables. However, it is not uncommon to have a third variable for which association between the same X variable is expected. Thus, a common scatter plot includes a second Y axis. These graphs are a bit more involved to make regardless of which application used. The purpose of this short section is to provide a way to create a plot with two Y axes against a common X axis. Additional plotting options are also provided and explained.

The data set

The first draft of this post was written in September 4, 2019, peak time for hurricanes. Hurricane Dorian passed the Bahamas as a category 4 storm heading for Florida. For context I include a screenshot of imagery from NOAA; you can see Dorian along the Florida coastline as well as four additional storms lined up from the coast of Africa across the Atlantic Ocean (Fig 1).

Figure 1. Screenshot from NOAA GOES-East – Sector view: Tropical Atlantic – GeoColor, 4 September 2019. Click image to view full sized image.

The data set here is the number of Atlantic Ocean Category 1, 2, 3, 4, or 5 storms since 1900, tabulated by decade. Storms are categorized by wind speed according to the Saffir-Simpson Hurricane Wind Scale, or SSHWS.

Additional data, average levels of carbon dioxide (CO₂), a greenhouse gas, measured at Mauna Loa and average global temperature index from NOAA were also acquired for the same period. (The temperature index is called anomaly data as it is the difference of temperature between the average temperatures between 1950 and 1980). 

#Create the time data, 12 decades, with 2010 based on 8 years only -- up to 2018
decade <- seq(1900,2010, by = 10)
#Saffir-Simpson Hurricane Wind Scale; compiled from Wikipedia, checked against NOAA
cat01 <- c(2,0,1,3,0,9,4,11,18,11,27,46)
cat02 <- c(2,0,1,3,10,4,5,4,3,13,9,8)
cat03 <- c(7,2,3,1,4,10,9,5,4,12,11,7)
cat04 <- c(1,6,5,8,9,10,9,5,5,10,15,8)
cat05 <- c(0,0,2,6,0,2,4,3,3,2,9,4)
#land ocean anomaly temperature index, NOAA
tempIndex <- c(-0.317,-0.329,-0.241,-0.123,0.042,-0.048,-0.028,0.034,0.247,0.387,0.591,0.791)
#mauna kea CO2, mean by decade, record starts 1958
co2 <- c(NA,NA,NA,NA,NA,316.0,320.3,330.9,345.5,360.5,378.6,399.0)

Combine all of the variables into a single data frame.

storms <- data.frame(decade,tempIndex,cat01,cat02,cat03,cat04,cat05,co2)
head(storms)

Output from R should look like

  decade tempIndex cat01 cat02 cat03 cat04 cat05 co2
1   1900    -0.317     2     2     7     1     0  NA
2   1910    -0.329     0     0     2     6     0  NA
3   1920    -0.241     1     1     3     5     2  NA
4   1930    -0.123     3     3     1     8     6  NA
5   1940     0.042     0    10     4     9     0  NA
6   1950    -0.048     9     4    10    10     2 316

Create a new variable with all categories of storms summed by decade. 

#Sum all of the storms
allHur <- cat01+cat02+cat03+cat04+cat05
#Add new variable to the data frame
storms$allHur <- allHur
#Attach the data frame so you don't have to keep referring to the data frame when you call variables.
attach(storms)

Now, create a simple scatter plot, number of hurricanes by decade (Fig 2).

plot(decade,allHur,pch=16)

Figure 2. Plot of hurricanes from 1900 to present by decade.

A note on data management

This note is not needed for the plots. However, it’s as good a place as any to comment about how to work with your data in R. For the hurricane data, I set the data frame storms as an unstacked (wide) worksheet. The code that follows shows how to use the reshape2 package, part of the tidyverse set of packages for data manipulation, to convert the data frame from unstacked (wide) to stacked (long).

#code modified from examples presented by Thomas Eisfeld
library(reshape2)
library(dplyr)
tryMe <- melt(storms, id.vars=c("decade", "tempIndex", "co2"), variable.name = "SSHWS", value.name = "hurricanes")
head(tryMe)

You can also take a stacked (long) worksheet and convert it to unstacked (wide) with dcast().

tryAgain <- dcast(tryMe, decade + tempIndex + co2 ~ SSHWS)

Make the plots

At last, here’s the code for adding a second Y axis. Code modified from https://www.r-bloggers.com/r-single-plot-with-two-different-y-axes/. The work flow begins by setting parameters of the plotting window, creating the first plot, followed up by adding the second plot (par=TRUE), setting graphing elements, then adding labels and a legend.

First plot example: Total number of hurricanes by decades, with Temperature Index by decades. Number of hurricanes represented on first (left) axis and Temperature Index represented on second (right) axis (Fig 3).

par(mar = c(5,5,2,5))
#Create the first plot
plot(tryAgain$decade, tryAgain$allHur, pch=16, cex=1.5, col="black", xlab="Decades", ylab="Hurricanes", main="Atlantic Hurricanes by Decade")

the syntax par(mar = c(bottom, left, top, right))

#Add the second plot
par(new=T)
plot(tryAgain$decade, tryAgain$tempIndex, type="l", lty="dashed", lwd=2, col="red3", axes=F, xlab=NA, ylab=NA)
axis(side=4)
mtext(side=4,line=3,"Temperature Index")

#Add a legend
legend("topleft",legend=c("Hurricanes", "Temperature Index"),lty=c(0,2), pch=c(16,NA), col=c("black", "red3"))

Figure 3. Total number of hurricanes by decades, with Temperature Index by decades. Number of hurricanes represented on first (left) axis and Temperature Index represented on second (right) axis.

Next plot example: Total number of hurricanes by decades, with Atmospheric CO₂ measured at Mauna Loa by decades. Number of hurricanes represented on first (left) axis and Atmospheric CO₂ represented on second (right) axis (Fig 4).

par(mar = c(5,5,2,5))
plot(tryAgain$decade, tryAgain$allHur, pch=16, cex=1.5, col="black", xlab="Decades", ylab="Hurricanes", main="Atlantic Hurricanes by Decade")

Create space for the second plot in the same frame (use T or TRUE, both work).

par(new=T)

Code for the second plot. Note that axes=F (“F” or “FALSE” works) is used to suppress printing axes — Recall that lines were already created by the first plot; if you do not add this, then additional lines are added to the plot.”NA” is used to suppress adding labels; if you do not add this code, then labels will be printed over the existing labels created for the first plot. The code axis=4 sets the right-hand Y axis to active. The next code, mtext() is used to place the axis label for the second Y axis.

plot(tryAgain$decade, tryAgain$co2, type="l", lty="dashed", lwd=2, col="red3", axes=F, xlab=NA, ylab=NA)
axis(side=4)
mtext(side=4,line=3,"Mean CO2, ppm")

Add the legend.

legend("topleft",legend=c("Hurricanes", "CO2"),lty=c(0,2), pch=c(16,NA), col=c("black", "blue"))

Figure 4. Total number of hurricanes by decades, with Atmospheric CO₂ measured at Mauna Kea by decades. Number of hurricanes represented on first (left) axis and Atmospheric CO₂ represented on second (right) axis.

Go ahead and try these plots. Change the settings (e.g., lty, pch, cex, col), and note how the graphs look.

Bubble (chart) plot

From the basic scatter plot a variety of variants have been developed to convey additional information about ratio scale variables. The bubble chart is one common example and can be used as an alternative to adding a second Y-axis in some cases. A bubble plot uses the area of a bubble to represent a third variable, while a second y-axis is used to display a second set of data with its own scale on the same chart.

In brief, a third variable is added, mapped to circle size. Spreadsheet apps also include options to make bubble charts and the wonderful plot.ly library provides a vehicle to produce sophisticated graphs including bubble charts online or as a package in R — and it works well with the ggplot2 package. A nice bubble chart introduction is available at Atlassian. 

Examples of bubble charts with R code are provided in Ch16.3. That said, what’s one more?  Here’s example R code:

ggplot(storms, aes(x = decade, y = tempIndex, size = total_storms)) +
  geom_point(alpha = 0.7, color = "lightgreen") + # Added color = "lightgreen"
  scale_size_continuous(range = c(3, 15)) +
  labs(title = "Temperature Index vs. Decade with Total Storms",
       x = "Decade",
       y = "Land Ocean Anomaly Temperature Index",
       size = "Total Storms") +
  theme_classic()

and the plot, Figure 4.

Figure 4. Bubble plot of association between total hurricanes and Land-Ocean Temperature Anomalies by Decade.

If you are interested in trying yourself, a bubble chart yourself, CO₂ datas with the hurricane data set is a good candidate for a bubble plot.

Questions

1. Using the plot() command, make the following graphs
   • One scatter plot for each category of storm by decade
   • Explore the kinds of graph elements available by changing pch values. Create your own 
   • Change point size by changing valued for cex
2. Create a new plot with Decade on X axis, Temperature Index on first Y axis, and CO₂ on the second Y axis.
   • Include a legend for the plot.
3. What is a Q-Q plot used for in statistics?
4. Looking at the plot in Figure 5, explain why the confidence lines get further and further away from the straight line.
5. What are some advantages and disadvantages of heat map for communicating data.
6. What color pallet is considered “color-friendly” for accessible visualization?
7. Create a table with three columns.

• In column 1, list the name for each graph type discussed in Chapter 5.
• In column 2, state the kind of data appropriate for each graph.
• In column 3, enter the name of a function in the base R package you can use to make the graph
• In column 4, enter the name of a function from the ggplot2 package you can use to make the graph

Quiz Chapter 4.6

Adding a second Y axis

Start Quiz

Question

Confirm

✕

✓

Your answer:

Correct answer:

Next

Please wait..

↺

You got {{SCORE_CORRECT}} out of {{SCORE_TOTAL}}

Your Answers

Chapter 4 contents

• Graphs and tables (How to report statistics)
• Bar (column) charts
• Histograms
• Box plot
• Mosaic plots
• Scatter plots
• Add a second Y axis
• Q-Q plot
• Ternary plots
• Heat maps
• Graph software
• References and suggested readingsh

4.2 – Histograms

Introduction

Displaying interval or continuously scaled data. A histogram (frequency or density distribution), is a useful graph to summarize patterns in data, and is commonly used to judge whether or not the sample distribution approximates a normal distribution. Three kind of histograms, depending on how the data are grouped and counted. Lump the data into sequence of adjacent intervals or bins (aka classes), then count how many individuals have values that fall into one of the bins — the display is referred to as a frequency histogram. Sum up all of the frequencies or counts in the histogram and they add to the sample size. Convert from counts to percentages, then the heights of the bars are equal to the relative frequency (percentage) — the display is referred to as a percentage histogram (aka relative frequency histogram). Sum up all of the bin frequencies and they equal one (100%).

Figure 1 shows two frequency histograms of the distribution of ages for female (left panel) and male (right panel) runners at the 2013 Jamba Juice Banana 5K race in Honolulu, Hawaii (link to data set).

Figure 1. Histograms of age distribution of runners who completed the 2103 Jamba Juice 5K race in Honolulu

The graphs in Fig 1 were produced using R package ggplot2.

The third kind of histogram is referred to as a probability density histogram. The height of the bars are the probability densities, generally expressed as a decimal. The probability density is the bin probability divided by the bin width (size). The area of the bar gives the bin probability and the total area under the curve sums to one. 

Which to choose? Both relative frequency  histogram and density histogram convey similar messages because both “sum to one (100%), i.e., bin width is the same across all intervals. Frequency histograms may have different bin widths, with more numerous observations, the bin width is larger than cases with fewer observations. 

Purpose of the histogram plot

The purpose of displaying the data is to give you or your readers a quick impression of the general distribution of the data. Thus, from our histogram one can see the range of the data and get a qualitative impression of the variability and the central tendency of the data. 

Kernel density estimation

Kernel density estimation (KDE) is a non-parametric approach to estimate the probability distribution function. The “kernel” is a window function, where an interval of points is specified and another function is applied only to the points contained in the window. The function applied to the window is called the bandwidth. The kernel smoothing function then is applied to all of the data, resulting in something that looks like a histogram, but without the discreteness of the histogram.

The chief advantage of kernel smoothing over use of histograms is that histogram plots are sensitive to bin size, whereas KDE plot shape are more consistent across different kernel algorithms and bandwidth choices.

Today, statisticians use kernel smoothing functions instead of histograms; these reduce the impact that binning has on histograms, although kernel smoothing still involves choices (smoothing function? default is Gaussian. Widths or bandwidths for smoothing? Varies, but the default is from the variance of the observations). Figure 2 shows smoothed plot of the 752 age observations.

Figure 2. KDE plot of age distribution of female runners who completed the 2103 Jamba Juice 5K race in Honolulu.

R commands typed at the R prompt were, in order

d <- density(w) #w is a vector of the ages of the 752 females
plot(d, main="Gaussian smoothing")
polygon(d, col="black", border="black") #col and border are settings which allows you to color in the polygon shape, here, black with black border.

Conclusion? A histogram is fine for most of our data. Based on comparing the histogram and the kernel smoothing graph I would reach the same conclusion about the data set. The data are right skewed, maybe kurtotic (peaked), and not normally distributed (see Ch 6.7).

Design criteria for a histogram

The X axis (horizontal axis) displays the units of the variable (eg, age). The goal is to create a graph that displays the sample distribution. The short answer here is that there is no single choice you can make to always get a good histogram — in fact, statisticians now advice you to use a kernel function in place of histograms if the goal is to judge the distribution of samples.

For continuously distributed data the X-axis is divided into several intervals or bins:

1. The number of bins (intervals) depends (somewhat) on the sample size and (somewhat) on the range of values. Thus, the shape of the histogram is dependent on your choice of intervals: too many bins and the plot flattens and stretches to either end (over smoothing); too few bins and the plot stacks up and the spread of points is restricted (under smoothing). For both you lose the
2. A general rule of thumb: try to have 10 to 15 different intervals. This number of intervals will generally give enough information.
3. For large sample size (N=1000 or more) you can use more intervals.

The intervals on the X-axis should be of equal size on the scale of measurement.

1. They will not necessarily have the same number of observations in each interval.
2. If you do this by hand you need to first determine the range of the data and then divide this number by the number of categories you want. This will give you the size of each category (eg, range is 25; 25 / 10 = 2.5; each category would be 2.5 unit).

For any given X category the Y-axis then is the Number or Frequency of individuals that are found within that particular X category.

Software

Rcmdr: Graphs → Histogram…

Accepting the defaults to a subset of the 5K data set we get (Fig 3)

Figure 3. Histogram of 752 observations, Sturge’s rule applied, default histogram.

The subset consisted of all females or n = 752 that entered and finished the 5K race with an official time.

R Commander plugin KMggplot2

Install the RcmdrPlugin.KMggplot2as you would any package in R. Start or restart Rcmdr and load the plugin by selecting Rcmdr: Tools → Load Rcmdr plug-in(s)… Once the plugin is installed select Rcmdr: KMggplot2 → Histogram… The popup menu provides a number of options to set to format the image. Settings for the next graph were No. of bins “Scott,” font family “Bookman,” Colour pattern “Set 1,” Theme “theme_wsj2.”

Figure 4. Histogram of 752 observations, Scott’s rule applied, ggplot2 histogram

Selecting the correct bin number

You may be saying to yourself, wow, am I confused. Why can’t I just get a graph by clicking on some buttons? The simple explanation is that the software returns defaults, not finished products. It is your responsibility to know how to present the data. Now, the perfect graph is in the eye of the beholder, but as you gain experience, you will find that the default intervals in R bar graphs have too many categories (recall that histograms are constructed by lumping the data into a few categories, or bins, or intervals, and counting the number of individuals per category => “density” or frequency). How many categories (intervals) should we use?

Figure 5. Default histogram with different bin size.

R’s default number for the intervals seems too much to me for this data set; too many categories with small frequencies. A better choice may be around 5 or 6. Change number of intervals to 5 (click Options, change from automatic to number of intervals = 5). Why 5? Experience is a guide; we can guess and look at the histograms produced.

Improving estimation of bin number

I gave you a rough rule of thumb. As you can imagine there have been many attempts over the years to come up with a rational approach to selecting the intervals used to bin observations for histograms. Histogram function in Microsoft’s Excel (Data Analysis plug-in installed) uses the square root of the sample size as the default bin number. Sturge’s rule is commonly used, and the default choice in some statistical application software (e.g., Minitab, Prism, SPSS). Scott’s approach (Scott 2009), a modification to Sturge’s rule, is the default in the ggplot() function in the R graphics package (the Rcmdr plugin is RcmdrPlugin.KMggplot2). And still another choice, which uses interquartile range (IQR), was offered by Freedman and Diacones (1981).  Scargle et al (1998) developed a method, Bayesian blocks, to obtain optimum binning for histograms of large data sets. 

What is the correct number of intervals (bins) for histograms?

• Use the square root of the sample size, e.g., at R prompt sqrt∗n∗ = 4.5, round to 5). Follow Sturges’ rule (to get the suggested number of intervals for a histogram, let k = the number of intervals, and k = 1 + 3.322(log₁₀n),

• where n is the sample size. I got k = 5.32, round to nearest whole number = 5).
• Another option was suggested by Freedman and Diacones (1981): find IQR for the set of observations and then the solution to the bin size is k = 2IQRn^_-1/3, where n is the sample size.

Select Histogram Options, then enter intervals. Here’s the new graph using Sturge’s rule…

Figure 6. Default histogram, bin size set by Sturge’s rule.

OK, it doesn’t look much better. And of course, you’ll just have to trust me on this — it is important to try and make the bin size appropriate given the range of values you have in order that the reader/viewer can judge the graphic correctly.

Two histograms on same plot

In many cases you may wish to compare sample distributions among two or more groups. 

Figure 7. Two histograms on same plot with ggplot2.

R code

Stacked data set, data from Kaggle.

myColors=c(female="blue", male="red")
ggplot(sexChol, aes(chol, fill = sex, color=sex)) + 
  geom_histogram(binwidth = 10, alpha=0.75) + 
  labs(x = "Cholesterol") + 
  theme_bw() + 
  theme(legend.title= element_blank(),legend.text=element_text(size=14), 
     axis.text=element_text(size=14),axis.title=element_text(size=18)) +
  scale_color_manual(values=myColors) + 
  scale_fill_manual(values=myColors)

Base graphics version

Figure 8. Same data as Fig 7, but using base hist() and plot() functions.

Not an aesthetically pleasing graph (Fig 8), but the idea was to recreate Figure 7 and show regions of overlap. We do this by setting transparent colors (h/t DataAnalytics.org.uk).

First, pick you colors. I like blue and red. Next, we identify how the colors are declared in the three rgb channels. Use the col2rgb() function.

col2rgb("blue")
[,1]
red 0
green 0
blue 255

Create my transparent blue with the alpha setting (125 = 50%). Note the command asks for names, but we won’t use it further.

myBlue <- rgb(0, 0, 255, max = 255, alpha = 125, names = "blue50")

Repeat for 50% transparent red.

col2rgb("red")
[,1]
red 255
green 0
blue 0
myRed <- rgb(255, 0, 0, max = 255, alpha = 125, names = "red50")

Now, we make the plots. Data unstacked, a column with cholesterol values for females (fChol), and another for male cholesterol values mChol). Bin size is different from Figure 7. 

hist1 <- hist(fChol, breaks=6, plot=FALSE)
hist2 <- hist(mChol, add=T, breaks=6, plot=FALSE)

Open a Graphics Device so that we can save image to png format. I’ve set width and height in numbers of pixels.

png("base_histograms_two_chol.png", width = 800, height = 600)

Make the graph

plot(hist1, col = myBlue, main="", xlab="Cholesterol", ylab="Counts")
plot(hist2, col = myRed, add = TRUE)

Close the device, the image file is saved to working folder.

dev.off()

Draw region on a histogram

From time to time it may be necessary to highlight a region on a histogram (eg., see Chapters 6.7 – Normal distribution and the normal deviate (Z) and 8 – Inferential statistics). Various solutions to this, but here’s example code using only base R graph options.

First, generate some data from a normal distribution with mean 50 and standard deviation of 10. We show two examples: highlight regions between 40 and 60, and highlight regions 

set.seed(123)
data <- rnorm(1000, mean = 50, sd = 10)

As always, a good idea to view at least a portion of the data, stored in object data.

head(data)

R returns

[1] 44.39524 47.69823 65.58708 50.70508 51.29288 67.15065

Next, create the histogram. It’s a good idea to set plot = FALSE, which prevents immediate plotting.

h <- hist(data, breaks = 20, plot = FALSE) 
plot(h, main = "Histogram with Highlighted Regions", 
     xlab = "Value", ylab = "Frequency")

And last, as the goal was to highlight a region from 40 to 60, we set color to red and impose 30% transparency, no border. The plot is shown in Figure 9.

rect(xleft = 40, ybottom = 0, xright = 60, ytop = max(h$counts),
col = rgb(1, 0, 0, 0.3), border = NA)

Figure 9. Using only base R graphics, a lot with region between 40 and 60 highlighted in red.

The second example, we highlight another region from 70 to 80 in blue, and this time 50% transparency, no border.

rect(xleft = 70, ybottom = 0, xright = 80, ytop = max(h$counts),
col = rgb(0, 0, 1, 0.5), border = NA) 

Other options

The base graphs are OK for teaching purposes (and homework answers!), but not examples of data visualization excellence. As part of the tidyverse approach to R programming, ggplot2 includes options to accomplish highlighting regions on plots, including histograms.To draw regions on a histogram plot using ggplot2 in R, you can use geom_rect() or geom_ribbon(). geom_rect() is a good choice to highlight a specific rectangular areas;  geom_ribbon() is useful for shading areas under a curve, such as a density plot overlaid on a histogram.

Figure 10. Using geom_ribbon() and ggplot2 package, a plot with region between 40 and 60 highlighted in red.

R code with # comments was:

set.seed(123)
data_normal <- data.frame(rnorm(1000, mean = 50, sd = 10))

# View the dataframe
head(data_normal)

# Output shows need to rename the column.
colnames(data_normal) <- "value"

Get the histogram. first, create the plot object, but don’t draw it yet.

myPlot <- ggplot(data_normal, aes(x = value)) +
geom_histogram(bins = 30)

Next, use ggplot_build to extract the bin data.

bin_data <- ggplot_build(myPlot)$data[[1]]

Create a new data frame for the geom_ribbon layer. The ribbon needs x, ymin, and ymax values.

ribbon_data <- bin_data

Set ymin to 0 for the shaded area.

ribbon_data$ymin <- 0

Next, set ymax to the height of the bars (the ‘count’ or ‘y’ column from bin_data).

ribbon_data$ymax <- ribbon_data$count

Filter the ribbon data so that only the region between 40 and 60 is included.

ribbon_data <- ribbon_data[ribbon_data$x >= 40 & ribbon_data$x <= 60, ]

Plot with both layers. Draw the main histogram first (optional: fill to distinguish). Overlay the shaded region using geom_ribbon with the filtered data. Add axes labels. I selected the classic theme; try other themes on your own.

ggplot(data_normal, aes(x = value)) + 
     geom_histogram(aes(y = after_stat(count)), 
     bins = 30, fill = "lightblue", color = "white") + 
     geom_ribbon(data = ribbon_data, aes(x = x, 
     ymin = ymin, ymax = ymax), fill = "red", 
     alpha = 0.5) +
     labs(x = "Value",y = "Count") +
     theme_classic()

A package tigerstats, now removed from CRAN, was a nice set of tools to accomplish adding regions to a histogram plot with base R graphics. I leave an example which was used to draw regions on a standard normal curve for Chapter 8.

library(tigerstats)
pnormGC(1.645, region="above", mean=0, sd=1,graph=TRUE)
pnormGC(-1.645, region="below", mean=0, sd=1,graph=TRUE)

Questions

Example data set, comet tails and tea

Figure 11. Examples of comet assay results.

The Comet assay, also called single cell gel electrophoresis (SCGE) assay, is a sensitive technique to quantify DNA damage from single cells exposed to potentially mutagenic agents. Undamaged DNA will remain in the nucleus, while damaged DNA will migrate out of the nucleus (Figure 11). The basics of the method involve loading exposed cells immersed in low melting agarose as a thin layer onto a microscope slide, then imposing an electric field across the slide. By adding a DNA selective agent like Sybr Green, DNA can be visualized by fluorescent imaging techniques. A “tail” can be viewed: the greater the damage to DNA, the longer the tail. Several measures can be made, including the length of the tail, the percent of DNA in the tail, and a calculated measure referred to as the Olive Moment, which incorporates amount of DNA in the tail and tail length (Kumaravel et al 2009). 

The data presented in Table 1 comes from an experiment in my lab; students grew rat lung cells (ATCC CCL-149), which were derived from type-2 like alveolar cells. The cells were then exposed to dilute copper solutions, extracts of hazel tea, or combinations of hazel tea and copper solution. Copper exposure leads to DNA damage; hazel tea is reported to have antioxidant properties (Thring et al 2011). 

Copy the table into a data frame.

1. Create histograms for tail, tail percent, and olive moment
   • Change bin size
2. Repeat, but with a kernel function.
3. Looking at the results from question 1 and 2, how “normal” (i.e., equally distributed around the middle) do the distributions look to you?
4. Plot means to compare Tail, Tail percent, and olive moment. Do you see any evidence to conclude that one of the teas protects against DNA damage induced by copper exposure?

Quiz Chapter 4.2

Histograms

Start Quiz

Question

✕

✓

Your answer:

Correct answer:

Next

Please wait..

↺

You got {{SCORE_CORRECT}} out of {{SCORE_TOTAL}}

Your Answers

Chapter 4 contents

• Graphs and tables (How to report statistics)
• Bar (column) charts
• Histograms
• Box plot
• Mosaic plots
• Scatter plots
• Add a second Y axis
• Q-Q plot
• Ternary plots
• Heat maps
• Graph software
• References and suggested readings

2.2 – Why do we use R Software?

If this page is TL;DR, and you’re in a rush, then instructions to install R onto your macOS or Windows PC are provided in Install R.  To run in the cloud only, see Use R in the cloud.

Why R? The case for using R in biostatistics education.

.

Why do we use R Software? Or put another way: DrD, Why are you making me use R?

Truth? You can use just about any acceptable statistical application to get the work done and achieve the learning objectives we have for beginning biostatistics. However, we will use the R statistical language as our primary statistical software in this course. Part of the justification is that all statistical software applications come with a learning curve, so you’d start at zero regardless of which application I used for the course.

In selecting software for statistics I have several criteria. The software should be:

• if not exactly easy, the software should have a reasonable learning curve
• widely accessible and compatible with all most personal computers
• well-respected and widely used by professionals
• free software
• open source
• well-supported for the purposes of data analysis and data processing
• really good for making graphics, from the basics to advanced
• capable to handle diverse kinds of statistical tests

R meets all of these criteria. R history began back in 1993 and has always been available as free software under the terms of the Free Software Foundation’s GNU General Public License in source code form. R compiles and runs on a wide variety of UNIX platforms and similar systems, including GNU/LINUX, FreeBSD, and various Linux distros like the popular Ubuntu^®, in addition to their more famous Microsoft Windows^® and Apple macOS^® distributions. To facilitate access to the software, numerous mirror sites are available from sites around the world, with cloud.r-project.org supported by RStudio perhaps the most widely used. From January 2024 to December 2024, more than 8.5 million downloads of base R were made from the RStudio CRAN mirror site (CRAN stands for Comprehensive R Archive Network.

install.packages("cranlogs")
library(cranlogs)
# How many downloads of base R around start of Fall semester?
 out <- cran_downloads("R", from = "2025-08-01", to = "2025-08-31")
sum(out$count)

[1] 582524

R is straightforward to use once you learn how to work with the language, but has a steep learning curve; after all, it’s a programming language. I recommend installing the program onto your computer, but ready to go R is also available at multiple cloud services, including Google CoLab and Posit Cloud. The remainder of this page is directed at users who choose to install R to their computer, although much of the discussion holds regardless of the R environment choice, local or Cloud.

For local installs of R, I also recommend use of the GUI R Commander.  Use of R Commander helps with the learning curve, and eventually, your use of code will become second nature. We don’t use RStudio, but as your skills improve, I recommend transitioning to the IDE RStudio; after the initial growing pains are behind you, RStudio likely will be a better solution over R Commander.

This discussion implies learning R programming. However, while we need statistical software to do statistics, students in my BI311 course must keep in mind that learning objectives for a biostatistics course is about the concepts and interpretation of statistics, not just use of the software. In other words, learning how to use R is not the focus of BI311 nor will you likely achieve R programming competency by the end of the semester. I certainly encourage students to strive for competency and I give frequent bonus opportunities to demonstrate coding skills during the semester.

Why R and not some other statistics software?

Thus, you might ask if the purpose of the course isn’t to learn R, why work with R instead of a more familiar app or software, eg, Microsoft Excel^® (hereafter simply referred to as Excel), or Google Sheets, or even my favorite open-source office alternative, LibreOffice Calc? Or, perhaps even just one of the many online calculators, if the course learning objective is to “just” learn about statistics?

First, I believe that real data derived from real biology or biomedical problems are essential elements to a first course in biostatistics. That’s not a particularly unique perspective although I don’t have survey results of other statistics instructors to back up the claim. Real problems involve observations on multiple subjects, many variables — large data sets; this alone precludes use of hand calculations and calculators. As a corollary, we will not spend a great deal of time learning the in’s and out’s of the algorithms that form particular statistical tests. Now, do understand that there is a tremendous benefit to understanding statistics by working through the equations, by looking at the algorithms, and there’s no escaping the need for understanding that probability provides the foundation of statistics inference (Chapter 8). Thus, for most of us, the statistical software available to us provides an appropriate framework for applying correct statistical tests to our projects. Therefore, the decision is about which statistical package we should use.

Second, R is perhaps the choice in academia for statistical software. A PUBMED search found more than 1500 citations of R. Visit Robert A. Muenchen’s web page (The popularity of data analysis software, r4stats.com) to see updated statistics on statistical software use. Those of you continuing on to graduate school or to professional schools will find that many of your statistically literate colleagues use R and not one of the commercial programs. While there are many excellent commercial packages (Table 1), and in some cases you can make spreadsheet programs do statistics (typically add-ins are required), all statistical software come with steep learning curves. Thus, part of my selling point to you is that learning to use R is at the cutting-edge in your field and, given that all of the software you could use can have have their challenges, it is best to work with something that will be around and is in wide use, without the burden of a financial investment.

Table 1. A selective list of statistical software.

see Wikipedia for list of additional software

Third, what about online sites like plot.ly where, for free, you can plot and, in some cases, calculate statistics? What about the web application at Brightstat, which claims to provide an SPSS-like experience online (Stricker 2008)? While it is true that there are many wonderful websites that can perform many of the statistical tests we will use this semester, these sites are not suitable for more than occasional use.

How to get started with R.

.

A few words about my approach to the software in my biostatistics course. I note to students during the semester that our course is not a programming course, but rather, programming skills are acquired along the way. Our focus is how to do statistical analysis, how to design experiments, etc. Coding skills without understanding of the why we’re doing it is not unlike learning to cook by following a great recipe as opposed to learning about the art and science of flavor (NPR Science Friday, March 15, 2024) — the homework is solved but the why we can support our conclusions may be obscure.

We set aside the first couple of class meetings simply to setup student’s Windows PC or Apple macOS computers to do data science; for students with access to Chromebooks or iPad tablets, we need to set them up to run R in the Cloud (eg., Google Colaboratory) — in other words, there’s a lot of pre-work that can present barriers to actual doing statistics in R (or other software). Once we have the environment set up for the students, we utilize a “tell-show-do” approach to learning how to code. The “tell” part includes sharing script; the “do” part is mostly done by students outside of class as part of completing homework.

First work with the language can be frustrating, so I also encourage a “five minute” rule — if it ain’t working, stop. So, by all means when Run through a troubleshooting checklist, but please don’t continue to struggle — the gap is between the instructions and where the student is with the material — it’s a me-as-instructor problem, not a student must work harder problem! 

All statistics work follows a similar workflow.

The basic workflow for a data project, using R as the primary statistical programming tool, looks like that illustrated in Figure 1.

Figure 1. A basic workflow with R.

The workflow begins with the assumption that a data set already exists. Getting the data into an R session depends on the source of the data, for example stored in spreadsheet file or from a table on a webpage. Data cleaning and processing of the data may be in order — for small and modest-sized data sets, my own preference is to do most of the data cleaning in my spreadsheet app, but R has full capabilities to handle or “wrangle” data. Next, the raw data should be stored in R session in a data frame object. Once the data frame object is available, exploratory data analysis, including summary statistics and data visualization can proceed. Higher order analytics, including statistical inference and modeling, can then begin.

The R statistical language, accompanied by additional packages to extend its capabilities beyond basic math and statistical functions, provided a complete statistical environment. R is best viewed as a programming language for statistics (data analysis), and data processing. Power users of R learn how to write scripts that do t-tests, ANOVA, regression, etc. The scripts are just lines of code that R understands and it provides the user tremendous control over analysis and inference of data sets. Because of this flexibility and power, however, R can be intimidating at first. So, we’ll start slowly with scripts, introducing just what we need to get started and build from there. We’ll be addressing R issues in more depth over the next several weeks, but for the first week, our goal(s) should be to make sure each of you knows how to start/exit R, how to create and utilize a working directory, and how to use R as a calculator. You obtain your copy of R from the R Project for Statistical Computing, available at https://www.r-project.org. Instructions to install R are provided in Install R.  A ten-part tutorial to get started using R is provided in Mike’s Workbook for Biostatistics.

 /Users/username/Desktop/BI311

Why R Commander, why not RStudio?

.

We utilize an R package that provides a menu-driven context to much of the typical statistics one needs to do biostatistics. The package is called R Commander (Rcmdr), which provides a graphic user interface or GUI. Rcmdr therefore significantly eases the learning curve for doing statistics with R. We use a package called R Commander, which provides drop down menus for most of the typical kinds of analyses. Rcmdr is in use in many courses across the world (more than 17K downloads in August 2025), and among the other GUI available for R, Rcmdr is among the best supported GUI available for R. R Commander function is extended by plug-ins; as of August 2025, there were 31 plugins that extend Rcmdr’s capabilities. Instructions to install R Commander are provided in Install R Commander. 

Wait! Why don’t we use Excel? My instructor in {insert course here} used Excel…

.

A very reasonable question for you to ask — why don’t we use Microsoft Excel or Google Sheets for statistics? Moreover, it is highly likely that you have gained at least some introduction to descriptive statistics and graphing with spreadsheets in former courses — shouldn’t we learn statistics within a software framework you are already familiar?

After all, “Can’t Microsoft Excel do statistics?” Mostly the answer is, no, not really (Fig 2).

Figure 2. “Spreadsheets,” xkcd.com no. 2180.

MS Excel, Google Sheets, Apple Numbers, and for that matter, Calc, the spreadsheet application in my favorite office app LibreOffice (LibreOffice is a free, open-source alternative to Microsoft Office), can be used to calculate many descriptive statistics. With some effort, these applications can be extended by use of either Analysis ToolPak or Solver Add-ins to do more complicated statistics like regression and analysis of variance, and curve fitting.

However, use of MS Excel for statistical analysis involves learning a number of commands, syntax, and developing workflows that are neither intuitive nor standard. Some publishers have provided add-ins that are reportedly designed to simplify this process (eg, XLStat, UNISTAT, Real Statistics using Excel). None of these options are free and none are in use in any major way by scientists (see Robert Muenchen’s The popularity of data analysis software). The free add-ins of Analysis ToolPak and Solver may work for you if you own a Windows PC, but only Solver is included for the Mac versions of Excel. MacOS users may download and install StatPlus:MacLE, which is a limited, but free alternative to the Analysis ToolPak add-in; for a complete package a Pro version is available (licenses started at $89, web site: www.analystsoft.com/en/products/statplusmacle/).

An additional caution: you should be aware that there have been reports over the years that algorithms selected by Microsoft for Excel have not always been to industry standards (eg, McCullogh and Wilson 2005). In short, the fit of Excel and other spreadsheet apps for use in statistics is not a simple one. To do the kinds of statistics we will use routinely in class, Excel would need to be modified with add-ins, and the add-ins would be the result of programming by someone. And you would still need to learn how to write the code.

What about graphics? You may like Microsoft Excel’s ability to do graphics. Indeed, Excel, Google Sheets, and LibreOffice Calc can be used to generate many typical kinds of statistical plots. But again, in comparison to R, spreadsheet app graphics are limited and require a deal of effort to generate acceptable plots. I think you’ll be surprised at how straight-forward R is. Here’s an example, first rendered in Microsoft Excel, then in base R. And importantly, the kinds of plots Excel does well at are not necessarily the plots suitable for research publication. For example, Excel allows you to make bar charts easily, but cannot do box plots. Box plots are preferred over bar (column) charts for ratio scale data.

Statistics comparisons between R and MS Excel.

.

About that learning curve. Let’s compare R and MS Excel for basic functions common in data analysis. Similar conclusions hold for comparisons to Google Sheets and LibreOffice Calcs.

Table 2 lists the observations we can use to conduct comparisons of the applications.

A simple data set

One of the first steps in data analysis is to produce what are called descriptive statistics. Common descriptive statistics are the mean and the sample standard deviation. Let’s compare Excel and R for retrieving these two statistics.

With Excel, to calculate the arithmetic mean of 24 numbers, enter the values into a single column of 24 rows, then enter “=average(A2:A25)“, without the quotes, into a new cell of the spreadsheet. “A2:A25” refers to where data would be contained in column A rows 2 through 25. Typically the first row in a worksheet would contain the name of the variable, eg, “A.” Depending on the significant figures set, the estimate returned by Excel for the mean of A is 59.58333333.

Similarly, to obtain the standard deviation, type =stdev(A2:A25), into a new cell of the spreadsheet. Again, depending on the significant figures set, Excel returns a value of 37.05215674 for the standard deviation of A.

In contrast, to obtain the mean and standard deviation for a variable in an R data set, all you would type at the R prompt (>), or in the script window

and then submit the code is:

myA <- c(12, 14, 20, 25, 28, 29, 32, 34, 35, 39, 47, 47, 50, 53, 54, 71, 79, 87, 89, 96, 105, 122, 130, 132)

where the “c” is a function to combine arguments into a vector and saved to the object myA, followed at the new line by

mean(myA)

Hit enter after entering the command) and R returns

[1] 59.58333

For the standard deviation, write the R base function sd()

sd(myA)

Hit enter after entering the command and R returns

[1] 37.05216

It’s not much of a difference, but note that to get the mean (arithmetic average) I typed seven characters in R, but 16 characters in Excel; similarly, for the standard deviation I typed in 5 characters in R, but 13 characters in Excel. That’s a savings of 56% and 62%, respectively. Excel tries to help by using AutoComplete to anticipate what you want to enter, but AutoComplete doesn’t always work properly (eg, see gene name errors generated by use of default Microsoft Excel settings, Ziemann et al 2016).

In conclusion, R is quicker for descriptive statistics.

Graphics comparison between R and MS Excel.

.

MS Excel is often cited for its graphics capabilities (Camões 2016). We can make the familiar scatter plots, bar charts, and pie charts in Excel. These plots and more are easily obtained in R. I won’t elaborate here about graphics, we talk at some length about graphics in Chapter 4. But here’s one example in R.

Let’s plot myB vs myA. We already provided the data for variable A, here’s the data for variable B.

17, 21, 21, 26, 27, 32, 28, 42, 40, 30, 71, 53, 56, 61, 55, 89, 82, 63, 116, 162, 116, 154, 137, 149

Don’t recall how to assign a set of numbers to an object, B, in R? See above and look again at how we assigned the numbers to object myA.

To get a simple scatter plot (Fig 3), I may write at the R prompt.

plot(myA,myB)

Figure 3. Basic scatter plot made in R, plot(myA,myB).

And here’s the comparable default plot (Fig 4) from Microsoft Excel, Office 365

Figure 4. Basic scatterplot made in Microsoft Excel.

Now, both graphs need some work, and to be fair, these are just the defaults. With some effort, you can make an Excel graph look pretty good. But note — the defaults in Excel don’t generate axis labels, while R default plot does. Excel adds a useless title and legend; both need to be removed. Excel also adds grid lines where typically one would not include these in a scientific plot; for example, in the case of scatterplots, grid lines can distract the viewer from the data points and reduce the data-to-ink ratio Tufte 1984).

Count the steps to generate an acceptable scatter plot (Table 3). I’ve also added R Commander (Rcmdr) steps for comparisons (Rcmdr lets you use drop-down menus like Excel or Google Sheets or LibreOffice Calcs).

Table 3. Steps needed to make a simple scatterplot in R, R Commander, or Microsoft Excel.

Conclusion? R is quicker for routine statistical plots like a scatter plot. And I didn’t even count the steps needed to change MS Excel’s dreadful diamond icon points.

That’s one step in R, four steps in Rcmdr, but eight steps for Microsoft Excel. LibreOffice Calc is a little better at only four steps, but like MS Excel, you’d need to change several components to the graph (Fig 5).

Figure 5. Basic scatterplot made in LibreOffice Calc.

plot(myA,myB, pch = 22)

So, you’re telling me I don’t need a spreadsheet application?

.

No, not at all. We use spreadsheets, and more generally, databases, to store data. Spreadsheets apps are designed to make data entry and data management approachable and efficient. They remain an important tool for researchers (Browman and Woo 2017).

R is not that great of a spreadsheet; packages are available to seamlessly tie your spreadsheet and database data to R via ODBC. We will routinely enter and manipulate data in MS Excel, then import the data into R for analysis.

Spreadsheet apps like MS Excel and Google Sheets (see also LibreOffice Calc) are great at being a spreadsheet program, R is great at being a statistical software program. You should take advantage of what the tools do best.

Still not convinced?

.

R is in use all over the world, by students and professionals alike and if one is going to spend the time to learn how to use a statistics software program, you should learn a standard program, like R.

And it’s not just me. Read about R in this 2009 New York Times piece, “Data analysts captivated by R’s power.” Look who purchased (April 2015) Revolution Analytics, a major player in the development of the R programming language.

Why install R on your computer?

.

Convenience. Control. Offline.

At the Biology department of Chaminade University, we have installed and maintain R, Rcmdr, and RStudio along with all required packages on our Macbook Pro® Lab computers for your use during class and during optional, proctored biostatistics work sessions. Since 2018, R is increasingly available “in the cloud” (eg, RStudio Cloud), which would mean you could run R in your browser and avoid installation on your computer. You can run significant analysis with R in the cloud via the free Google Colaboratory and CoCalc are now available: I encourage you to look into these platforms. Unfortunately, these services are not quite ready for the classroom. For example, RStudio in the Cloud is free to use on a limited basis, but quickly requires a significant subscription cost with increasing use. Google Colab and GoCalcs require use of Jupyter notebooks, which add yet another layer to the learning curve without focusing on learning statistics. Second, although access to their servers is easy, running simultaneous connections via Chaminade’s single public IP address is likely to lead to problems for us. Third, I want you to use R Commander (Rcmdr) to assist in the learning curve — Rcmdr cannot be run in the Cloud (ie, RStudio in the Cloud, Google Colaboratory, or CoCalc).

Therefore, you are encouraged to install R, Rcmdr, and even RStudio, onto your own computers, in part because of the convenience, but also because R is not generally available to students on campus, ie, only the Biology department’s computers have the up-to-date R software installed.

To get started, go to your Canvas website and view How to install R on your own computer.

An additional benefit to installing a version of R on your computer, you’ll understand more about the software if you take the time to install and if need be, troubleshoot your installation of the software. Moreover, there’s a considerable amount of help out there for R. For example, a simple Google search(keywords: tutorial “install R”), returns more than 700K hits, and more than 40K January 2023 alone (add “after:2023-01-01” to Google search box). In fact, there’s so much out there that you’ll want to sample from several sites and select the voice that works best for you.  

Questions

1. Write up three learning outcomes for this page. Hint: Point your favorite generative AI to this page and ask for help.
2. I listed several spreadsheet apps: Microsoft Excel, Google Sheets, Apple Numbers, LibreOffice Calc. Which of these are free to download and install to your computer? Which are freely available via the Cloud?
3. What level of confidence do you have to this statement: I am confident in my ability to use spreadsheet apps for my BI-311 work this semester. Response options: 
   1. Strongly disagree.
   2. Disagree.
   3. Agree.
   4. Strongly agree.
4. What level of confidence do you have to this statement: I am confident in my ability to use the R programming software for my BI-311 work this semester. Response options: 
   1. Strongly disagree.
   2. Disagree.
   3. Agree.
   4. Strongly agree.

.

Quiz Chapter 2.2

Why do we use R software?

Start Quiz

Question

✕

✓

Your answer:

Correct answer:

Next

Please wait..

↺

You got {{SCORE_CORRECT}} out of {{SCORE_TOTAL}}

Your Answers

Chapter 2 contents

.

• Introduction
• Why biostatistics?
• Why do we use R Software?
• A brief history of (bio)statistics
• Experimental Design and rise of statistics in medical research
• Scientific method and where statistics fits
• Statistical reasoning
• Chapter 2 – References and suggested reading