R Crash Course
Mark Dunning
Last modified: 27 Feb 2023
Introduction to R - Part III
Lets make sure we have read the gapminder data into R and have the relevant packages loaded.
library(readr)
library(dplyr)
library(ggplot2)
gapminder <- read_csv("raw_data/gapminder.csv")Customising a plot
Now make a scatter plot of gdp versus life expectancy as we did in the previous session. One of the last topics we covered was how to add colour to a plot. This can make the plot more appealing, but also help with data interpretation. In this case, we can use different colours to indicate countries belonging to different continents.
ggplot(gapminder, aes(x = gdpPercap, y=lifeExp,col=continent)) + geom_point()
The shape and size of points can also be mapped from the data. However, it is easy to get carried away.
ggplot(gapminder, aes(x = gdpPercap, y=lifeExp,shape=continent,size=pop)) + geom_point()
Scales and their legends have so far been handled using
ggplot2 defaults. ggplot2 offers functionality
to have finer control over scales and legends using the scale
methods.
Scale methods are divided into functions by combinations of
the aesthetics they control.
the type of data mapped to scale.
scale_aesthetic_type
Try typing in scale_ then tab to autocomplete. This will
provide some examples of the scale functions available in
ggplot2.
Although different scale functions accept some variety in their arguments, common arguments to scale functions include -
name - The axis or legend title
limits - Minimum and maximum of the scale
breaks - Label/tick positions along an axis
labels - Label names at each break
values - the set of aesthetic values to map data values
We can choose specific colour palettes, such as those provided by the
RColorBrewer package. This package provides palettes for
different types of scale (sequential, diverging, qualitative).
library(RColorBrewer)
display.brewer.all(colorblindFriendly = TRUE)
When creating a plot, always check that the colour scheme is appropriate for people with various forms of colour-blindness
When experimenting with colour palettes and labels, it is useful to save the plot as an object
p <- ggplot(gapminder, aes(x = gdpPercap, y=lifeExp,col=continent)) + geom_point()# We can also change the text displayed above the legend with the name parameter.
p + scale_color_manual(values=brewer.pal(6,"Set2"))
Or we can even specify our own colours; such as The University of Sheffield branding colours
my_pal <- c(rgb(0,159,218,maxColorValue = 255),
rgb(31,20,93,maxColorValue = 255),
rgb(249,227,0,maxColorValue = 255),
rgb(0,155,72,maxColorValue = 255),
rgb(190,214,0,maxColorValue = 255))
p + scale_color_manual(values=my_pal)
NEW:- A set of palettes based on works in the Metropolitan Museum of Art (New York) has been made available.
https://github.com/BlakeRMills/MetBrewer
## this will check if MetBrewer is already installed, and will only install if it is not found
if(!require("MetBrewer")) install.packages("MetBrewer")
library(MetBrewer)
p + scale_color_manual(values=met.brewer(name = "Greek"))
Various labels can be modified using the labs
function.
p + labs(x="Wealth",y="Life Expectancy",title="Relationship between Wealth and Life Expectancy")
We can also modify the x- and y- limits of the plot so that any
outliers are not shown. ggplot2 will give a warning that
some points are excluded.
p + xlim(0,60000)Warning: Removed 5 rows containing missing values (`geom_point()`).
Saving is supported by the ggsave function. A variety of
file formats are supported (.png, .pdf,
.tiff, etc) and the format used is determined from the
extension given in the file argument. The height, width and
resolution can also be configured. See the help on ggsave
(?ggsave) for more information.
ggsave(p, file="my_ggplot.png")Saving 7 x 5 in imageMost aspects of the plot can be modified from the background colour
to the grid sizes and font. Several pre-defined “themes” exist and we
can modify the appearance of the whole plot using a
theme_.. function.
p + theme_bw()
More themes are supported by the ggthemes package. You
can make your plots look like the Economist, Wall Street Journal or
Excel (but please don’t do this!)
## this will check if ggthemes is already installed, and will only install if it is not found
if(!require("ggthemes")) install.packages("ggthemes")
library(ggthemes)
p + theme_excel()
Exercise
Exercise: Use a boxplot to compare the life
expectancy values of Australia and New Zealand. Use a Set2
palette from RColorBrewer to colour the boxplots and apply
a “minimal” theme to the plot.
Facets
One very useful feature of ggplot2 is faceting. This
allows you to produce plots subset by variables in your data. In the
scatter plot above, it was quite difficult to see if the relationship
between gdp and life expectancy was the same for each continent. To
overcome this, we would like a see a separate plot for each
continent.
To facet our data into multiple plots we can use the
facet_wrap (1 variable) or facet_grid (2
variables) functions and specify the variable(s) we split by.
p + facet_wrap(~continent)
The facet_grid function will create a grid-like plot
with one variable on the x-axis and another on the y-axis.
p + facet_grid(continent~year)
The previous plot was a bit messy as it contained all combinations of year and continent. Let’s suppose we want our analysis to be a bit more focused and disregard countries in Oceania (as there are only 2 in our dataset) and maybe years between 1997 and 2002.
We should know how to restrict the rows from the
gapminder dataset using the filter function.
Instead of filtering the data, creating a new data frame and
constructing the data frame from these new data we can use
the%>% operator to create the data frame on the fly and
pass directly to ggplot. Thus we don’t have to save a new
data frame or alter the original data.
filter(gapminder, continent!="Oceania", year %in% c(1997,2002,2007)) %>%
ggplot(aes(x = gdpPercap, y=lifeExp,col=continent)) + geom_point() + facet_grid(continent~year)
Summarising and grouping with dplyr
The summarise function can take any R function that
takes a vector of values (i.e. a column from a data frame) and returns a
single value. Some of the more useful functions include:
minminimum valuemaxmaximum valuesumsum of valuesmeanmean valuesdstandard deviationmedianmedian valueIQRthe interquartile rangen_distinctthe number of distinct valuesnthe number of observations (Note: this is a special function that doesn’t take a vector argument, i.e. column)
library(dplyr)
summarise(gapminder, min(lifeExp), max(gdpPercap), mean(pop))It is also possible to summarise using a function that takes more than one value, i.e. from multiple columns. For example, we could compute the correlation between year and life expectancy. Here we also assign names to the table that is produced.
gapminder %>%
summarise(MinLifeExpectancy = min(lifeExp),
MaximumGDP = max(gdpPercap),
AveragePop = mean(pop),
Correlation = cor(year, lifeExp))However, it is not particularly useful to calculate such values from
the entire table as we have different continents and years. The
group_by function allows us to split the table into
different categories, and compute summary statistics for each year (for
example).
gapminder %>%
group_by(year) %>%
summarise(MinLifeExpectancy = min(lifeExp),
MaximumGDP = max(gdpPercap),
AveragePop = mean(pop))gapminder %>%
group_by(year,continent) %>%
summarise(MinLifeExpectancy = min(lifeExp),
MaximumGDP = max(gdpPercap),
AveragePop = mean(pop))`summarise()` has grouped output by 'year'. You can override using the `.groups` argument.We can list as many summary functions as we like. Whilst this can make our code somewhat verbose there are many helper functions available. Consider an example where we want to average all the columns in our data:-
gapminder %>%
group_by(year) %>%
summarise(MeanLifeExpectancy = mean(lifeExp),
MeanGDP = mean(gdpPercap),
MeanPop = mean(pop))This wasn’t a huge effort to write this code. However, it would be
much more tedious for a dataset with many more columns. Recognising
this, we can use the convenient summarise_all function.
This will return NA values for columns that do not contain
numeric values.
gapminder %>%
group_by(continent) %>%
summarise_all(mean)The nice thing about summarise is that it can followed
up by any of the other dplyr verbs that we have met so far
(select, filter, arrange..etc).
As the country column of the previous output containing
missing values we can exclude it from further processing.
gapminder %>%
group_by(continent) %>%
summarise_all(mean) %>%
select(-country)Returning to the correlation between life expectancy and year, we can summarise as follows:-
gapminder %>%
group_by(country) %>%
summarise(Correlation = cor(year , lifeExp))We can then arrange the table by the correlation to see which countries have the lowest correlation
gapminder %>%
group_by(country) %>%
summarise(Correlation = cor(year , lifeExp)) %>%
arrange(Correlation)We can filter the results to find observations of interest
gapminder %>%
group_by(country) %>%
summarise(Correlation = cor(year , lifeExp)) %>%
filter(Correlation < 0)The countries we identify could then be used as the basis for a plot.
library(ggplot2)
filter(gapminder, country %in% c("Rwanda","Zambia","Zimbabwe")) %>%
ggplot(aes(x=year, y=lifeExp,col=country)) + geom_line()
Exercise
- Summarise the
gapminderdata in an appropriate manner to produce a plot to show the change in averagegdpPercapfor each continent over time. - see below for a suggestion
- HINT: you will need to use the
geom_colfunction to create the bar plot
- HINT: you will need to use the
Joining
In many real life situations, data are spread across multiple tables or spreadsheets. Usually this occurs because different types of information about a subject, e.g. a patient, are collected from different sources. It may be desirable for some analyses to combine data from two or more tables into a single data frame based on a common column, for example, an attribute that uniquely identifies the subject.
dplyr provides a set of join functions for combining two
data frames based on matches within specified columns. For those
familiar with such SQL, these operations are very similar to carrying
out join operations between tables in a relational database.
As a toy example, lets consider two data frames that some results of testing whether genes A, B and C are significant in our study (gene expression, mutations, etc.)
gene_results <- data.frame(Name=LETTERS[1:3], pvalue = c(0.001, 0.1,0.01))
gene_resultsWe might also have a data frame containing more data about the genes;
such as which chromosome they are located on. As part of our data
interpretation we might need to know where in the genome the genes are
located. Note that both data frames have a column called
Name. This column will be used to identify genes common to
both tables.
gene_anno <- data.frame(Name = c("A","B","D"), chromosome=c(1,1,3))
gene_annoThere are various ways in which we can join these two tables together. We will first consider the case of a “left join”.
Animated gif by Garrick Aden-Buie
left_join returns all rows from the first data frame
regardless of whether there is a match in the second data frame. Rows
with no match are included in the resulting data frame but have
NA values in the additional columns coming from the second
data frame.
Animations to illustrate other types of join are available at https://github.com/gadenbuie/tidy-animated-verbs
left_join(gene_results, gene_anno)Joining, by = "Name"right_join is similar but returns all rows from the
second data frame that have a match with rows in the first data frame
based on the specified column.
right_join(gene_results, gene_anno)Joining, by = "Name"inner_join only returns those rows where matches could
be made
inner_join(gene_results, gene_anno)Joining, by = "Name"Wrap-up
We have introduced a few of the essential packages from the R tidyverse that can help with data manipulation and visualisation.
Hopefully you will feel more confident about importing your data into R
and producing some useful visualisations. You will probably have
questions regarding the analysis of your own data. Some good starting
points to get help are listed below.
To finish the workshop we will look at the analysis of some relevant data that we can import into R and analyse with the tools from the workshop.
Introducing the COVID-19 data
Data for global COVID-19 cases are available online from CSSE at Johns Hopkins University on their github repository.
github is an excellent way of making your code and analysis available for others to reuse and share. Private repositories with restricted access are also available. Here is a useful beginners guide.
R is capable of downloading files to our own machine so we can analyse them. We need to know the URL (for the COVID data we can find this from github, or use the address below) and can specify what to call the file when it is downloaded.
download.file("https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv",destfile = "raw_data/time_series_covid19_confirmed_global.csv")We can use the read_csv function as before to import the
data and take a look. We can see the basic structure of the data is one
row for each country / region and columns for cases on each day.
covid <- read_csv("raw_data/time_series_covid19_confirmed_global.csv")Rows: 289 Columns: 1136
── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (2): Province/State, Country/Region
dbl (1134): Lat, Long, 1/22/20, 1/23/20, 1/24/20, 1/25/20, 1/26/20, 1/27/20, 1/28/20, 1/29/20, 1/30/20, 1/31/20, 2/1/20, 2/2/20, 2/3/20, 2/4/20, 2/5/20, 2/6/20, 2/7/20, 2/8/20, 2/9/20, 2/10/20, 2/11/20, 2/12/20, 2/13/...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.covidWe can potentially join these data to gapminder, and it
would be beneficial to have one column name in common between both
files. We can rename the Country/Region column
of our new data frame to match gapminder.
covid <- read_csv("raw_data/time_series_covid19_confirmed_global.csv") %>%
rename(country = `Country/Region`) Rows: 289 Columns: 1136
── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (2): Province/State, Country/Region
dbl (1134): Lat, Long, 1/22/20, 1/23/20, 1/24/20, 1/25/20, 1/26/20, 1/27/20, 1/28/20, 1/29/20, 1/30/20, 1/31/20, 2/1/20, 2/2/20, 2/3/20, 2/4/20, 2/5/20, 2/6/20, 2/7/20, 2/8/20, 2/9/20, 2/10/20, 2/11/20, 2/12/20, 2/13/...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.covidMuch of the analysis of this dataset has looked at trends over time
(e.g. increasing /decreasing case numbers, comparing trajectories). As
we know by now, the ggplot2 package allows us to map
columns (variables) in our dataset to aspects of the plot.
In other words, we would expect to create plots by writing code such as:-
ggplot(covid, aes(x = Date, y =...)) + ...Unfortunately such plots are not possible with the data in it’s
current format. Counts for each date are containing in a different
column. What we require is a column to indicate the date, and the
corresponding count in the next column. Such data arrangements are known
as long data; whereas we have wide data. Fortunately
we can convert between the two using the tidyr package
(also part of tidyverse).
## install tidyr if you don't already have it
install.packages("tidyr")For more information on tidy data, and how to convert between long and wide data, see
library(tidyr)
covid <- read_csv("raw_data/time_series_covid19_confirmed_global.csv") %>%
rename(country = `Country/Region`) %>%
pivot_longer(5:last_col(),names_to="Date", values_to="Cases")Rows: 289 Columns: 1136
── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (2): Province/State, Country/Region
dbl (1134): Lat, Long, 1/22/20, 1/23/20, 1/24/20, 1/25/20, 1/26/20, 1/27/20, 1/28/20, 1/29/20, 1/30/20, 1/31/20, 2/1/20, 2/2/20, 2/3/20, 2/4/20, 2/5/20, 2/6/20, 2/7/20, 2/8/20, 2/9/20, 2/10/20, 2/11/20, 2/12/20, 2/13/...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.covidAnother point to note is that the dates are not in an internationally recognised format, which could cause a problem for some visualisations that rely on date order. We can fix by explicitly converting to YYYY-MM-DD format.
For more ways of dealing with dates in R see the
lubridate package.
covid <- read_csv("raw_data/time_series_covid19_confirmed_global.csv") %>%
rename(country = `Country/Region`) %>%
pivot_longer(5:last_col(),names_to="Date", values_to="Cases") %>%
mutate(Date=as.Date(Date,"%m/%d/%y"))Rows: 289 Columns: 1136
── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (2): Province/State, Country/Region
dbl (1134): Lat, Long, 1/22/20, 1/23/20, 1/24/20, 1/25/20, 1/26/20, 1/27/20, 1/28/20, 1/29/20, 1/30/20, 1/31/20, 2/1/20, 2/2/20, 2/3/20, 2/4/20, 2/5/20, 2/6/20, 2/7/20, 2/8/20, 2/9/20, 2/10/20, 2/11/20, 2/12/20, 2/13/...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.covidAnother useful modification is to make sure only one row exists for each country. If we look at the data for some countries (e.g. China and UK) there are different entries for provinces and oversees territories.
## the count function tabulates the number of observations in a particular column
filter(covid, country == "China") %>%
count(`Province/State`)So we can change the Cases to be the sum of all cases
for that country on a particular day. We can do this using the
group_by and summarise functions from
above
covid <- read_csv("raw_data/time_series_covid19_confirmed_global.csv") %>%
rename(country = `Country/Region`) %>%
pivot_longer(5:last_col(),names_to="Date", values_to="Cases") %>%
mutate(Date=as.Date(Date,"%m/%d/%y")) %>%
group_by(country,Date) %>%
summarise(Cases = sum(Cases))Rows: 289 Columns: 1136
── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (2): Province/State, Country/Region
dbl (1134): Lat, Long, 1/22/20, 1/23/20, 1/24/20, 1/25/20, 1/26/20, 1/27/20, 1/28/20, 1/29/20, 1/30/20, 1/31/20, 2/1/20, 2/2/20, 2/3/20, 2/4/20, 2/5/20, 2/6/20, 2/7/20, 2/8/20, 2/9/20, 2/10/20, 2/11/20, 2/12/20, 2/13/...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
`summarise()` has grouped output by 'country'. You can override using the `.groups` argument.covidExercise
What plots and summaries can you make from these data?
- Plotting the number of cases over time for certain countries
- Which country in each continent currently has the highest number of cases?
- Normalise the number of cases for population size (using 2007
population figures as a population estimate)?
- e.g. cases per 100,000
- Which European countries have the highest number of cases per 100,000 population