• Science and humanities are increasingly data-driven
  • Early-career training has not prepared all researchers for this

• Enable systematic, replicable and reproducible work
  • Design principles
    • Best practices for data
  • Software development methods
    • Automation of repetitive calculations

• What a computer does
  • A series of instructions
  • Data is piped through programs, and a result emerges

• What a researcher does
  • Exploring data, developing hypotheses, writing code, interpreting results
• Outputs include:
  • datasets, methods, teaching materials, software, papers, etc.

• Reproduced from Stoudt et al. (2021)

• Fundamentals of R and RStudio
• Fundamentals of programming (in R)
• Data management with the tidyverse
• Publication-quality data visualisation with ggplot2
• Reporting with RMarkdown

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

• free (though commercial support can be bought)
• widely used
  • sciences, humanities, engineering, statistics, etc.
• has many excellent specialised packages for data analysis and visualisation
• international, friendly user community

• RStudio is an integrated development environment (IDE)

• 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, with extremely large datasets…

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

x <- 1 / 40
x

## [1] 0.025

x ^ 2

## [1] 0.000625

log(x)

## [1] -3.688879

name <- "Samia"
name

## [1] "Samia"

current_temperature = 28.6
subjectID = "GCF_00001236452.1"
GPS_Location = "54N, 36E"

• descriptive, but not too long
• letters, numbers, underscores, and periods ([a-zA-z0-9_.])
• cannot contain whitespace or start with a number (x2 is allowed, 2x is not)
• case sensitive (Weight is not the same as weight)
• do not reuse names of built-in functions
• Consistent style:
  • lower_snake, UPPER_SNAKE, lowerCamelCase, UpperCamelCase

• automate complicated tasks
• make code more readable and reusable

• Some functions are built-in (in base packages, e.g. sqrt(), lm(), plot())
• Groups of related functions can be imported as 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

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

• 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 a separate subfolder (data?)
• Clean data ready for work programmatically
  • keep cleaned/modified data in separate folder (clean_data?)
• Consider output generated by analysis to be disposable
  • can be regenerated by running analysis/code

• RStudio tries to help you manage your projects
  • R Project concept - files and subdirectory structure
  • integration with version control
  • 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

• 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://github.com/swcarpentry/r-novice-inflammation/raw/main/data/r-novice-inflammation-data.zip

(the link is also 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)

Someone gives you a data file that has:

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

How would you open it, using read.csv()

• We use indexing to refer to elements of a matrix
  • 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 row (patient) 2
min(data[, 7])      # smallest value on column (day) 7
mean(data[, 7])     # mean value on day 7
sd(data[, 7])       # standard deviation of values on day 7

• Calculating for every patient (or day) this way is tedious

• So apply a function (mean) to each row in the data:

• R has several ways to automate this process

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 plots to your script showing:

• scatterplot of standard deviation of inflammation across all patients, by day
• a histogram of average inflammation across all patients, 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 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
• data.frame

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

What data type are the following vectors (xx, yy, zz)?

xx <- c(1.7, "a")
yy <- c(TRUE, 2)
zz <- 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 comes as 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

• 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"
> l_named$b
[1] 1 2 3 4

• Logical indexes 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]

• 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   Length:3           Mode :logical  
##  1st Qu.:1.5   Class :character   FALSE:1        
##  Median :2.0   Mode  :character   TRUE :2        
##  Mean   :2.0                                     
##  3rd Qu.:2.5                                     
##  Max.   :3.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

In R:

• a package is a collection (or library) 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

• 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

• 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()

• 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:

• 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?

“Tidy datasets are all alike, but every messy dataset is messy in its own way”

• Data Cleaning is not just a first step
  • repeated when new data turns up, new ideas arrive, etc.
• About 80% of the effort of data analysis is cleaning and preparing data

df1 <- data.frame(treatment=c("treatmenta", "treatmentb"),
                  John.Smith=c(NA, 2),
                  Jane.Doe=c(16, 11),
                  Mary.Johnson=c(3, 1))
colnames(df1) = c("", "John Smith", "Jane Doe", "Mary Johnson")

df2 <- data.frame(name=c("John Smith", "Jane Doe", "Mary Johnson"),
                 treatmenta=c(NA, 16, 3),
                 treatmentb=c(2, 11, 1))
colnames(df2) = c("", "treatmenta", "treatmentb")

• A dataset is a collection of VALUES

• VARIABLES: can change or vary
  • values that are measured or decided by a researcher: height, temperature, duration, treatment, etc.
• OBSERVATIONS: values measured across all variables for the same individual/unit/group: person, reactor, religion, company, etc.

• Rows=OBERVATIONS, Columns=VARIABLES
• Rows=OBERVATIONS, Columns=NEITHER
• Rows=NEITHER, Columns=VARIABLES
• Rows=NEITHER, Columns=NEITHER

• VARIABLES
  • person (John, Mary, Jane)
  • treatment (a and b)
  • result (NA, 16, 3, 2, 11, 1)

• Tidy data is a standard way of structuring a dataset
  • not the only way
  • makes it easy to extract needed varaibles
  • well suited to R (supports vectorisation)

• 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 in Tidy Data
  • 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 in your script by including pipes) to produce a dataframe from gapminder containing:

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))

• 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))

• 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

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

• We often want to perform operations (or not) conditional on whether something is TRUE
  • 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
}

• 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
}

• 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

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\)

• 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 best

• 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 <- data %>% 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)
}

• ADVANCED: Make the facet wrapping optional

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