http://pad.software-carpentry.org/2017-12-04-standrews

• Fundamentals of R and RStudio
• Fundamentals of programming
• Best practices for organising code
• Best practices for reproducibility
• Effective data analysis in R

• Understand what R is
• Understand what RStudio is
• Understand why R is different from something like Excel

• R is a programming language
• R is the software that interprets programs written in the R programming language

• R is free (commercial support can be bought)
• R is widely-used and interdisciplinary (sciences, humanities, engineering, etc.)
• R has many excellent packages for statistics and data analysis, visualisation and graphics
• R has an international, friendly user community

• Separates data from analysis
• Not point-and-click: every step is explicit and transparent
• Easy to share, adapt, reuse, publish analyses with new/modified data (GitHub)
• R can be run on supercomputers…

• RStudio is an integrated development environment (IDE) - all platforms

• Familiarity with the RStudio IDE
• Introduce R syntax
• Learn good project management practices
• Set up a working directory with version control (git) in RStudio

• Try out some of your own calculations in the interactive console

2 ** 16

## [1] 65536

15 %% 4

## [1] 3

15 %/% 4

## [1] 3

• An item of data goes in the box (called Name)
• When we refer to the box (variable) by its name, we really mean what’s in the box

name <- "Samia"
name

## [1] "Samia"

x <- 1 / 40
x

## [1] 0.025

x ^ 2

## [1] 0.000625

log(x)

## [1] -3.688879

• Function (log(), sin() etc.) ≈ “canned script”
  • automate complicated tasks
  • make code more readable and reusable
• Some functions are built-in (in base packages, e.g. sqrt(), lm(), plot())
• Some functions are imported in libraries

args(fname)            # arguments for fname
?fname                 # help page for fname
help(fname)            # help page for fname
??fname                # any mention of fname
help.search("text")    # any mention of "text"
vignette(fname)        # worked examples for fname
vignette()             # show all available vignettes

To remove variables from your workspace, use the rm() function

x <- 1
ls()

## [1] "name" "x"

rm(x, name)
ls()

## character(0)

To remove all variables use the broom in the Environment tab (Rstudio), or:

rm(list=ls())

What will be the value of each variable after each statement in the following program?

mass <- 47.5
age <- 122
mass <- mass * 2.3
age <- age - 20

• mass = 47.5, age = 102
• mass = 109.25, age = 102
• mass = 47.5, age = 122
• mass = 109.25, age = 122

• Descriptive (but not too long)
• Avoid existing names (e.g. mean, matrix, list, etc.)
• Consistent style
  • periods.between.words
  • underscores_between_words
  • camelCaseToSeparateWords

• letters, numbers, underscores, and periods ([a-zA-z0-9_.])
• cannot start with a number
• whitespace is not allowed

• Use a single working directory per project/analysis
  • easier to move, share, and find files
  • use relative paths to locate files
• Treat raw data as read-only
  • keep in separate subfolder (data)?
• Clean data ready for work
  • keep cleaned/modified data in separate folder?
• Consider output generated by analysis to be disposable
  • can be regenerated by running analysis/code
  • don’t place under version control

• RStudio tries to help you manage your projects
  • R Project concept - files and subdirectory structure
  • git integration
  • switching between multiple projects within RStudio
  • stores project history

• At the console (you’ve done this)
• In a script
• As an interactive notebook
• As a markdown file
• As a Shiny app

We’re going to create a new dataset and R script.

• Putting code in a script makes it easier to modify, share and run

• Load data into an R project
• Produce summary statistics of data
• Extract subsets of data
• Basic plotting in R (base graphics)

• Patients have been given a new treatment for arthritis
• We have measurements of inflammation over a period of days for each patient
• We want to produce a preliminary analysis and graphs for this data

• https://swcarpentry.github.io/r-novice-inflammation/files/r-novice-inflammation-data.zip

(the link is available on the course Etherpad page)

• You created data manually earlier, but this is rare
• Data are most commonly read in from plain text files

read.csv(file = "data/inflammation-01.csv", header = FALSE)

How would you open a similar data file that had:

• a comma (,) as the decimal point character
• semi-colon (;) as the field separator

using read.csv()

• We use indexing to refer to elements of a variable
  • square brackets: []
  • row, then column: [row, column]

data[1, 1]     # First value in dataset
data[30, 20]   # Middle value of dataset

• To get a range of values, use the : separator (meaning ‘to’)

data[1:4, 1:4]   # rows 1 to 4; columns 1 to 4

• To select a complete row or column, leave it blank

data[5, ]     # row 5
data[, 16]    # column 16

• R provides useful functions to summarise data
• We can use indexing to get summary information on individual patients and days

max(data)           # largest value in dataset
max(data[2, ])      # largest value for patient 2
min(data[, 7])      # smallest value on day 7
mean(data[, 7])     # mean value on day 7
sd(data[, 7])       # standard deviation of values on day 7

Given the vector in R:

animal <- c("m", "o", "n", "k", "e", "y")

Can you generate slices that do the following:

• return the first three characters
• return the last three characters

• return the first three characters in reverse order

• What do animal[-1] and animal[-4] do?

• Can you explain what animal[-1:-4] does?

• We could calculate mean inflammation for every patient (or day) this way, but it’s tedious

• R has several ways to automate this process
• We’d like to apply a function (mean) to each row in the data:

apply(X = data, MARGIN = 1, FUN = mean)

rowMeans(data)
colMeans(data)

“The purpose of computing is insight, not numbers.” - Richard Hamming

• Base graphics are powerful tools for visualisation and understanding

plot(avg_inflammation_patient)

max_day_inflammation <- apply(dat, 2, max)
plot(max_day_inflammation)

plot(apply(dat,2,min))       # 3 functions in one!

Can you add a plot to your script showing:

• standard deviation, by day
• a histogram of inflammation by day

• Basic data types in R
• Common data structures in R
• How to find out the type/structure of R data
• Understand how R’s data types and structures relate to your own data

• R is mostly used for data analysis
• R has special types and structures to help you work with data
• Much of the focus is on tabular data (data frames)

• Data types in R are atomic
  • All data is one of these types
  • All data structures are built from these

Create examples of data with the following characteristics:

For each variable, test that it has the data type you intended

• vector
• factor
• list
• matrix
• data.frame

Vectors are atomic: they can contain only a single data type

What data type are the following vectors (xx)?

xx <- c(1.7, "a")
xx <- c(TRUE, 2)
xx <- c("a", TRUE)

• Coercion means changing data from one type to another
• R will perform implicit coercion on vectors to make them atomic

• Manual coercion with as.<type_name>()

Data often comes in one of two types:

• quantitative: e.g. integers or real numbers
  (weight <- 17.2; rooms <- 7)
• categorical: e.g. ordered or unordered classes
  (grade <- "8", coat <- "brindled")

• Factors are special vectors that represent categorical data
• Stored as vectors of labelled integers
• Cannot be treated as strings/text

Create a new factor, defining control and case experiments, and inspect the result:

f <- factor(c("case", "control", "case", "control", "case"))
str(f)

##  Factor w/ 2 levels "case","control": 1 2 1 2 1

• Using the help available to you in RStudio, can you create a factor with the same values, but where the control level is numbered 1?

• Matrices are 2D vectors of atomic values
  • An extremely important data type in numerical analyses

# Create matrix of zeroes
m1 <- matrix(0, ncol = 6, nrow = 3)

# Create matrix of numbers 1 and 2
m2 <- matrix(c(1, 2), ncol = 4, nrow = 3)

Can you create a matrix with:

• 5 columns
• 10 rows
• Containing the numbers 1:50

• lists are like vectors, but can hold any combination of datatype
  • elements in a list are denoted by [[]] and can be named

# create a list
l <- list(1, 'a', TRUE, matrix(0, nrow = 2, ncol = 2), f)
l_named <- list(a = "SWC", b = 1:4)

• We have used indexing, slicing and names to get data by ‘location’

> animal[c(2,4,6)]
[1] "o" "k" "y"
> m2[2:3, 3:4]
     [,1] [,2]
[1,]    2    1
[2,]    1    2
> l_named$b
[1] 1 2 3 4

• Logical indexes can select data that meets certain criteria

x <- c(5.4, 6.2, 7.1, 4.8, 7.5)
mask <- c(TRUE, FALSE, TRUE, FALSE, TRUE)
x[mask]
x[x > 7]

• Understand the concept of a data.frame
• Understand how a data.frame is built from R data structures
• Know how to access any element of a data.frame
• Read data into a data.frame
• Write data out from a data.frame

• The cats data is a data.frame

> class(cats)
[1] "data.frame"
> cats
    coat weight likes_string
1 calico    2.1            1
2  black    5.0            0
3  tabby    3.2            1

• A named list of vectors having identical lengths.
  • Each column is a vector
  • Each vector can be a different data type

# Create a data frame
df <- data.frame(a=c(1,2,3), b=c('eeny', 'meeny', 'miney'),
                 c=c(TRUE, FALSE, TRUE))
summary(df)

##        a           b         c          
##  Min.   :1.0   eeny :1   Mode :logical  
##  1st Qu.:1.5   meeny:1   FALSE:1        
##  Median :2.0   miney:1   TRUE :2        
##  Mean   :2.0                            
##  3rd Qu.:2.5                            
##  Max.   :3.0

I made some mistakes when defining this data.frame.

Can you spot and fix them?

author_book <- data.frame(author_first = c('Charles', 'Ernst', "Theodosius"),
                          author_last = c(Darwin, Mayr, Dobzhansky),
                          year = c(1942, 1970))

Can you predict the class for each column in the following example?

country_climate <- data.frame(country=c("Canada", "Panama",
                                        "South Africa", "Australia"),
                               climate=c("cold", "hot", "temperate",
                                         "hot/temperate"),
                               temperature=c(10, 30, 18, "15"),
                               northern_hemisphere=c(TRUE, TRUE, FALSE,
                                                     "FALSE"),
                               has_kangaroo=c(FALSE, FALSE, FALSE, 1))

Can you create the following data frame, but make b a vector of character elements, rather than a factor?

df_chr <- data.frame(a=c(1,2,3), b=c('eeny', 'meeny', 'miney'),
                     c=c(TRUE, FALSE, TRUE))

• We bind vectors to add columns, and lists to add rows
  • Adding a new factor level may be required

df <- cbind(df, vals = 3:1)
levels(df$b) <- c('eeny', 'meeny', 'miney', 'mo')
df <- rbind(df, list(4, 'mo', FALSE, 0))

write.table(df, "data/df_example.tab", sep="\t")

We need to provide

• the data.frame
• the path to the file being written
• a column separator

gapminder <- read.table("data/gapminder-FiveYearData.csv", sep=",", header=TRUE)

• R can also read data direct from the internet

url <- paste("https://raw.githubusercontent.com/resbaz/r-novice-gapminder-files/",
             "master/data/gapminder-FiveYearData.csv", sep = '')
gapminder <- read.table(url, sep=",", header=TRUE)

str(gapminder)              # structure of the data.frame
typeof(gapminder$year)      # data type of a column
length(gapminder)           # length of the data.frame
nrow(gapminder)             # number of rows in data.frame
ncol(gapminder)             # number of columns in data.frame
dim(gapminder)              # number of rows and columns in data.frame
colnames(gapminder)         # column names from data.frame
head(gapminder)             # first few rows of dataframe
summary(gapminder)          # summary of data in data.frame columns

• data.frames are lists and subset in the same way
• data.frames are 2D data (tabular) and subset like matrixes

gapminder[3]                  # single column (get dataframe)
gapminder[["lifeExp"]]        # single column (get vector/factor)
gapminder$year                # single column (get vector/factor)
gapminder[1:3,]               # row slice (get dataframe)
gapminder[3,]                 # row slice (get dataframe)
gapminder[, 3]                # column slice (get vector/factor)
gapminder[, 3, drop=FALSE]    # column slice (get dataframe)

I made some mistakes when subsetting gapminder, can you fix them?

# Extract observations collected for the year 1957
gapminder[gapminder$year = 1957,]

# Extract all columns except 1 through 4
gapminder[, -1:4]

# Extract all rows where life expectancy is greater than 80 years
gapminder[gapminder$lifeExp > 80]

# ADVANCED: Extract rows for years 2002 and 2007
gapminder[gapminder$year == 2002 | 2007]

• Understand what packages are
• How to install a desired package
• How to use a package in your code

In R:

• a package is the basic unit of reusable code
• many useful and specialist tools are distributed as packages
• over 10,000 packages are available at CRAN
• you can distribute your own code as a package

installed.packages()               # see installed packages
install.packages("packagename")    # install a new package
update.packages()                  # update installed packages
library(packagename)               # import a package for use in your code

Can you check if the following packages are installed on your system, and install them if necessary?

dplyr
ggplot2
knitr

• Be able to use ggplot2 to generate publication-quality graphics
• Understand the grammar of graphics

• ggplot2 is part of the Tidyverse, a collection of packages for data science
  • ggplot2 is the graphics package

• You can use ggplot2 like base graphics
  • qplot() ≈ plot()

library(ggplot2)
plot(gapminder$lifeExp, gapminder$gdpPercap, col=gapminder$continent)
qplot(lifeExp, gdpPercap, data=gapminder, colour=continent)

• Each observation in the data is a point
• A point’s aesthetics determine how it is rendered
  • co-ordinates on the image; size; shape; colour
• aesthetics can be constant or mapped to variables
• Many different plots can be generated from the same data by changing aesthetics

• aesthetics define datapoint representations as a new dataset

• If data are drawn as points: scatterplot
• If data are drawn as lines: line plot
• If data are drawn as bars: bar chart

# Generate plot of GDP per capita against life Expectancy
p <- ggplot(data=gapminder, aes(x=lifeExp, y=gdpPercap, color=continent))
p + geom_point()
p + geom_line()

Can you create another figure in your script showing how life expectancy changes as a function of time, as a scatterplot?

• We’ve just used another “Grammar of Graphics” concept: layers
  • ggplot2 plots are built as layers

• Data and aesthetics can be defined in a base ggplot object
  • values from the base are inherited by the other layers
  • the base can be overridden in other layers

• Data and aesthetics can be defined in a base ggplot object
  • values from the base are inherited by the other layers
  • the base can be overridden in other layers

p <- ggplot(data=gapminder, aes(x=lifeExp, y=gdpPercap, colour=continent))
p + geom_point()

• Data and aesthetics can be defined in a base ggplot object
  • values from the base are inherited by the other layers
  • the base can be overridden in other layers

p <- ggplot(data=gapminder, aes(x=lifeExp, y=gdpPercap, colour=continent))
p + geom_line(aes(group=country))

• We can use several layers of geoms to build a plot
  • alpha controls opacity for a layer

p <- ggplot(data=gapminder, aes(x=lifeExp, y=gdpPercap, color=continent))
p + geom_line(aes(group=country)) + geom_point(alpha=0.4)

Can you create another figure in your script showing how life expectancy changes as a function of time, coloured by continent, with two layers:

• a line plot, grouping points by country
• a scatterplot showing each data point, with 35% opacity

(commit your changes to the gapminder.R script)

• Data transformations are handled with scale layers

• Some geom layers transform the dataset
  • Usually this is a data summary (e.g. smoothing or binning)

• So far all our plots have all data in a single figure
• Comparisons can be clearer with multiple panels:
  • facets
  • “small multiples plots”

# Compare life expectancy over time by continent
p <- ggplot(data=gapminder, aes(x=year, y=lifeExp, colour=continent,
                                group=country))
p <- p + geom_line() + scale_y_log10()
p + facet_wrap(~continent)

Can you create a scatterplot and contour densities of GDP per capita against population size, with colour filled by continent?

ADVANCED: Transform the x axis to better visualise data spread, and use facets to panel density plots by year.

• How to manipulate data.frames with the six verbs of dplyr
  • a ‘grammar of data manipulation’

• dplyr is another package in the Tidyverse
• Facilitates analysis by groups
  • Helps avoid repetition

> mean(gapminder[gapminder$continent == "Africa", "gdpPercap"])
[1] 2193.755
> mean(gapminder[gapminder$continent == "Americas", "gdpPercap"])
[1] 7136.11
> mean(gapminder[gapminder$continent == "Asia", "gdpPercap"])
[1] 7902.15

library(dplyr)

select(gapminder, year, country, gdpPercap)
gapminder %>% select(year, country, gdpPercap)

• filter() selects rows on the basis of some condition

filter(gapminder, continent=="Europe")

# Select gdpPercap by country and year, only for Europe
eurodata <- gapminder %>%
              filter(continent == "Europe") %>%
              select(year, country, gdpPercap)

Can you write a single line (which may span multiple lines and include pipes) that produces a dataframe containing:

• life expectancy, country, and year data
• only for African nations

How many rows does the dataframe have?

group_by(gapminder, continent)
gapminder %>% group_by(continent)

# Produce table of mean GDP by continent
gapminder %>%
    group_by(continent) %>%
    summarize(meangdpPercap=mean(gdpPercap))

• Can you calculate the average life expectancy per country in the gapminder data?
• Which nation has longest life expectancy, and which the shortest?

• Two useful functions related to summarize()
  • count()/tally(): a function that reports a table of counts by group
  • n(): a function used within summarize(), filter() or mutate() to represent count by group

gapminder %>%
  filter(year == 2002) %>%
  count(continent, sort = TRUE)

gapminder %>%
  group_by(continent) %>%
  summarize(se_lifeExp = sd(lifeExp)/sqrt(n()))

• mutate() is a function allowing creation of new variables

# Calculate GDP in $billion
gdp_bill <- gapminder %>%
  mutate(gdp_billion = gdpPercap * pop / 10^9)

# Calculate total/sd of GDP by continent and year
gdp_bycontinents_byyear <- gapminder %>%
    mutate(gdp_billion=gdpPercap*pop/10^9) %>%
    group_by(continent,year) %>%
    summarize(mean_gdpPercap=mean(gdpPercap),
              sd_gdpPercap=sd(gdpPercap),
              mean_gdp_billion=mean(gdp_billion),
              sd_gdp_billion=sd(gdp_billion))

• How to make data-dependent choices in R
• Use if() and else()
• Repeat operations in R
• Use for() loops
• vectorisation to avoid repeating operations

• We often want to perform operations (or not) conditional on a piece of data
  • The if() and else construct is useful for this

# if
if (condition is true) {
  PERFORM ACTION
}

# if ... else
if (condition is true) {
  PERFORM ACTION
} else {  # i.e. if the condition is false,
  PERFORM ALTERNATIVE ACTION
}

Can you use an if() statement to report whether there are any records from 2002 in the gapminder dataset?

Can you do the same for 2012?

• for() loops are a very common construct in programming
  • for each <item> in a group, <do something (with the item)>
• Not as useful in R as in some other languages

for(iterator in set of values){
  do a thing
}

• while() loops are useful when you need to do something while some condition is true

while(this condition is true){
  do a thing
}

Can you use a for() loop and an if() statement to print whether each letter in the alphabet is a vowel?

• for() and while() loops can be useful, but are not efficient
• Most functions in R are vectorised
  • When applied to a vector, apply to all elements in that vector
  • No need to loop

x < 1:4
x * 2
y <- 6:9
x + y
m <- matrix(1:12, nrow = 3, ncol = 4)
m * m
m %*% m

We want to sum the following series of fractions

\(\frac{1}{1^2} + \frac{1}{2^2} + \frac{1}{3^2} + \ldots + \frac{1}{n^2}\)

for large values of \(n\)

Can you do this using vectorisation for \(n = 10,000\)?

• Why functions are important
• How to write a new function
• Defining a function that takes arguments
• Returning a value from a function
• Set default values for function arguments

• Functions let us run a complex series of commands in one go
  • under a memorable/descriptive name
  • invoked with that name
  • with a defined set of inputs and outputs
  • to perform a logically coherent task

• Small functions with one obvious, clearly-defined task are good

• You will often need to write your own functions
• They take a standard form

<function_name> <- function(<arg1>, <arg2>) {
  <do something>
  return(<result>)
}

my_sum <- function(a, b) {
  the_sum <- a + b
  return(the_sum)
}

• So far, you’ve been able to use R’s built-in help to see function documentation
  • This isn’t available for your functions unless you write it

• State what the code does (and why)
• Define inputs and outputs
• Give an example

• We can define functions that take multiple arguments
• We can also define default values for arguments

# Calculate total GDP in gapminder data
calcGDP <- function(data, year_in=NULL, country_in=NULL) {
  gdp <- gapminder %>% mutate(gdp=(pop * gdpPercap))
  if (!is.null(year_in)) {
    gdp <- gdp %>% filter(year %in% year_in)
  }
  if (!is.null(country_in)) {
    gdp <- gdp %>% filter(country %in% country_in)
  }
  return(gdp)
}

Can you write a function that takes an optional argument called letter, which:

• Plots the life expectancy per year for each country
• Only for countries whose name starts with a letter in letter

• Uses facet_wrap() to produce a grid of output graphs

• ADVANCED: Make the facet wrapping optional

starts.with <- substr(gapminder$country, start = 1, stop = 1)
az.countries <- gapminder[starts.with %in% c("A", "Z"), ]

• Create dynamic, reproducible reports
• Markdown syntax
• Inline R code in documents
• Producing documents in .pdf, .html, etc.

• A programming paradigm introduced by Donald Knuth
• The program (or analysis) is explained in natural language
  • The source code is interspersed
• The whole document is executable

• R Markdown files embody Literate Programming in R
• File \(\rightarrow\) New File \(\rightarrow\) R Markdown
• Enter a title
• Save the file (gets the extension .Rmd)

• Header information is fenced by ---

---
title: "Literate Programming"
author: "Leighton Pritchard"
date: "04/12/2017"
output: html_document
---

• Natural language is written as plain text

This is an R Markdown document. Markdown is a simple formatting syntax

• R code (which is executable) is fenced by backticks

• We’re going to create an R Markdown report on the gapminder data

You’ve learned: