Overview of the data

We will primarily be using two data sets:msleep,WOLF.CSV. Some examples will also use the datasetMALIN_HEAD.TXT. These data are described athttp://www.ucd.ie/ecomodel/Resources/datasets.html

The wolf data set was used in the“Importing data into R”tutorial and the msleep data set was used in the“Getting started with R”worksheet.

Import data

Start by loading the two data sets. See the data import tutorial for more details on this.

# Import data ------------------

# Load the 'ggplot2' package (required to access msleep dataset)
library('ggplot2')

wolf = read.table('WOLF.CSV', header=T, sep=',') # Import wolf data set
data(msleep)                                     # Load msleep data set

ALWAYScheck that the data have been imported correctly. For this you can use functions such as:ls(),head(),tail()and `summary()

Once the data are imported we create a subset of the data that only contains data from populations 1 and 2. We also create a new qualitative variable (factor) that codes for hunting intensity in the different populations (light hunting intensity in population 1, heavy hunting intensity in population 2).

# Subset the wolf data frame and remove unwanted levels
wolf.sub = droplevels(subset(wolf, Population!=3))


# Create a new variable called Hunting
wolf.sub$Hunting = 'Heavy'
wolf.sub$Hunting[wolf.sub$Population==1] = 'Light'

# Make the variable Hunting a factor
wolf.sub$Hunting = as.factor(wolf.sub$Hunting)

Central tendency

(opens in a new window)Central tendencydescribes the middle of the data, where most observations lie.

Two measures of central tendency are the mean and the median.

# Measures of central tendency for all cortisol concentrations
# The argument na.rm=TRUE removes any missing values (NA's)
mean(wolf.sub$Cpgmg, na.rm=TRUE)               # Mean
median(wolf.sub$Cpgmg, na.rm=TRUE)             # Median
quantile(wolf.sub$Cpgmg, prob=0.5, na.rm=TRUE) # Median = 50% quantile

The mean is commonly used, but it is overly influenced by extreme outliers

Example: Does the mean of the numbers 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 1000 reflect the ‘middle’ of the data?

The median is amore robustestimator of central tendency that is less affected by extreme outliers (e.g. repeat the above example using the median instead of the mean).

Spread

(opens in a new window)Spreaddescribes how far the data are scattered around the central tendency.

Three measures of spread are the standard deviation, the inter-quartile range (IQR) and the median absolute deviation (MAD).

# Measures of spread for all cortisol concentrations
sd(wolf.sub$Cpgmg, na.rm=T)                  # (Sample) Standard deviation
IQR(wolf.sub$Cpgmg, na.rm=T)                 # Inter-quartile range
mad(wolf.sub$Cpgmg, na.rm=T)                 # Median absolute deviation

The standard deviation quantifies spread about the mean. Like the mean it is influenced by outliers. The square of the standard deviation is called the variance.

Both the inter-quartile range and the median absolute deviation quantify spread about the median. Both are more robust to outliers than the standard deviation.

The inter-quartile range is defined as the difference between the 25% and 75% quantiles.

summary()of a data frame

Applying thesummary()function to a data frame will automatically gives many of these descriptive statistics.

# Summarise the wolf.sub data frame
summary(wolf.sub)   

   Individual         Sex              Population       Colour         
 Min.   :  1.00   Length:148         Min.   :1.000   Length:148        
 1st Qu.: 37.75   Class :character   1st Qu.:1.000   Class :character  
 Median : 74.50   Mode  :character   Median :2.000   Mode  :character  
 Mean   : 74.50                      Mean   :1.696                     
 3rd Qu.:111.25                      3rd Qu.:2.000                     
 Max.   :148.00                      Max.   :2.000                     
                                                                       
     Cpgmg           Tpgmg            Ppgmg        Hunting   
 Min.   : 4.75   Min.   : 3.250   Min.   :12.76   Heavy:103  
 1st Qu.:12.16   1st Qu.: 4.378   1st Qu.:19.50   Light: 45  
 Median :15.38   Median : 5.030   Median :25.00              
 Mean   :16.61   Mean   : 5.617   Mean   :25.89              
 3rd Qu.:19.98   3rd Qu.: 6.067   3rd Qu.:30.01              
 Max.   :40.43   Max.   :15.130   Max.   :53.28              
                                  NA's   :79                 

Note: the 1st and 3rd quartiles (1st Qu. and 3rd Qu.) are the 25% and 75% quantiles respectively. The 2nd quartile is the median (50% quantile).

Other descriptive statistics

Some functions for other descriptive statistic:

# Other descriptive statistics
quantile(wolf.sub$Cpgmg, prob=0.25, na.rm=T) # Quantiles (here 25% quantile)
max(wolf.sub$Cpgmg, na.rm=T)                 # Maximum
min(wolf.sub$Cpgmg, na.rm=T)                 # Minimum
range(wolf.sub$Cpgmg, na.rm=T)               # Range = Maximum - Minimum

Themomentspackage must be loaded for descriptive statistics of skewness and kurtosis.

library(moments)                            # Load the moments package

skewness(wolf.sub$Cpgmg, na.rm=T)           # Skewness of cortisol data

Note: In the tutorial on distributions we’ll see that skew is the asymmetry of a distribution about its mean

Correlation

Acorrelationis a relationship betweentwo variables. The strength of a relationship can be quantified using a correlation coefficient (we will look at this in more detail in the tutorial on regression).

Two common correlation coefficients arePearson correlation coefficientandSpearman’s rank correlation coefficient.

Here is an example of calculating correlation coefficients between variablesCpgmgandPpgmgin thewolf.subdata frame.

# Calculate Pearson's correlation coefficient
cor(wolf.sub$Cpgmg, wolf.sub$Ppgmg, use='complete.obs', method='pearson')

# Calculate Spearman's rank correlation coefficient
cor(wolf.sub$Cpgmg, wolf.sub$Ppgmg, use='complete.obs', method='spearman')

Both these coefficients can take values between -1 and 1.

• Correlations close to -1 are very strong negative relationships
• Correlations close to 0 are weak relationships
• Correlations close to 1 are very strong positive relationships

Contingency Tables

Qualitative variables can be described by contingency tables. Given a qualitative variable (e.g. types of mammals: herbivore, omnivore, …) a contingency table records how many observations there are in each category. A contingency table can easily display up to two qualitative variables. Higher numbers of variables quickly become difficult to represent (e.g. three qualitative variables would require the equivalent of a three dimensional table).

The functiontable()can be used to create contingency tables. Below are examples with the msleep data frame

# A one dimensional contingency table
# The distribution of observations across mammalian foraging types
table(msleep$vore)

  carni   herbi insecti    omni 
     19      32       5      20 

A two dimensional contingency table showing the number of species from each feeding type in each category of conservation status (cd=critically endangered, en=endangered, lc=least concern, nt=near threatened, vu=vulnerable).

# A two dimensional contingency table
table(msleep$conservation, msleep$vore)

              
               carni herbi insecti omni
  cd               1     1       0    0
  domesticated     2     7       0    1
  en               1     2       1    0
  lc               5    10       2    8
  nt               1     3       0    0
  vu               4     3       0    0

Barplots

A barplot is a one way to displaying how data are distributed across aqualitative variable. Each bar represents one category of the variable and the height of each bar usually represents the number of data points in each category.

A barplot is a graphical representation of a contingency table.

Below we use the packageggplot2to produce a bar plot of the feeding types of species in themsleepdata.

(opens in a new window)Video Tutorial:The anatomy of a ggplot command (3 mins)

library(ggplot2)      # Load ggplot2 graphics library

# Define labels for the bars
labs = c('herbi'='Herbivore',
         'carni'='Carnivore',
         'omni'='Omnivore',
         'insecti'='Insectivore')

# Use ggplot2 to produce a ...
# Bar plot of the distribution of feeding types in the msleep data set
ggplot(data=msleep,                   # Define the data to plot
       aes(x=vore)) + 
  geom_bar() +                        # Draw a bar plot
  scale_x_discrete(labels=labs) +     # Set x-axis labels
  labs(x='Feeding Type',              # Set axis titles
       y='Number of Species') + 
  theme_bw()                          # Set background to white

Notice that we have relabeled the bars with more meaningful names. Graphical output should be styled so that it is easy to read.

We will see how the font size of the text can be increased later. For now we will tweak the plot to use different colours to fill the bars.

# Bar plot of the distribution of feeding types in the msleep data set
# with different colours for the different bars
ggplot(data=msleep,                             # Define the data to plot
       aes(x=vore, 
           fill=vore)) +                        # Define how bars are coloured
  geom_bar() +                                 
  scale_x_discrete(name = 'Feeding Type',       
                   labels=labs) + 
  scale_fill_discrete(name = 'Feeding Type',    # Define the legend's title and labels
                      labels=labs) + 
  labs(x='Feeding Type',                       
       y='Number of Species') + 
  theme_bw() 

ggplot2 help pagesare available at(opens in a new window)https://ggplot2.tidyverse.org

R’s basic plotting commands also has thebarplot()function to produce barplots. Try running the command below.

# R's base function for producing a barplot

# Bar plot of the distribution of feeding types in the msleep data set
barplot(table(msleep$vore), 
        names.arg=c('Carnivore','Herbivore','Insectivore','Omnivore'),
        col='darkgrey',
        xlab='Feeding Type',ylab='Number of Species')

Pie charts

Pie charts are sometimes used for displaying how a whole quantity is divided into proportions. They are another way of representing the data’s distribution across a qualitative variable. Each section of a pie chart is one category of the qualitative variable and the size of the section represents the proportion of data points out of the entire data set that have that value.

Displaying the above barplot as a pie chart using R’spie()function gives,

# Display a pie chart

# Plot the proportions of foraging types in the msleep dataset
pie(table(msleep$vore))

One disadvantage of this chart is that it does not communicate the number of data points used to calculate the proportions.

Make your visualisation have a big impact on people by making them accessible to as many people as possible.

When producing visualisations be aware that different people will see your visualisations differently. Some examples of reasons for perceptual differences are:

• Viewing a black & white print out of your visualisation
• Colour blindness:(opens in a new window)https://www.nature.com/articles/d41586-021-02696-z
• People with low vision often magnify images. Magnified images can appear pixelated and hard to understand if they are too low resolution.

Google have come up with some nice visualisation design principles:(opens in a new window)https://medium.com/google-design/redefining-data-visualization-at-google-9bdcf2e447c6

Using RStudio

If you are using RStudio you can use its GUI to save graphical output to a file. Here is a short video tutorial on this:

(opens in a new window)Video Tutorial:Saving R graphics to a file using RStudio (3 mins)

Using ggplot2

To save a plot to a file you can use theggsave()command. Here is the code to save the coloured bubble plot (above).

# Save the figure to a file
ggsave(filename = 'coloured_plot.png',      # Set filename    
       width=8,                             # Set width (inches)
       height=5)                            # Set height (inches)

Violin plots

Violin plots are similar to boxplots, but the boxes are replaced by violin shapes that represent the distribution of the data. Below is the violin plot version of the boxplot above for cortisol in populations 1 & 2 (this violin plot was produced using packageggplot2). The boxplot version is also shown in grey inside each violin.

# Create a violin plot and overlay a box and whiskers plot


# Define axis label and title font size to be 18 
format = theme(axis.title = element_text(size=18), 
               axis.text = element_text(size=18)) 

ggplot(data=wolf.sub,           # Define data to plot
       aes(x=Hunting, 
           y=Cpgmg)) + 
  geom_violin(fill='green4') +  # Draw violin plot in green
  geom_boxplot(colour='black',  # Draw box and whiskers plot 
               fill='grey', 
               width=0.1) + 
  theme_bw() +                  # Set background to white
  format +                      # Set the look of the plot
  labs(x='Population',          # Set axis titles
       y='Cortisol (pg per mg)')

## Heatmap code
library(tidyr)
library(ggplot2)
# Import rainfall data
rainfall = read.table('MALIN_HEAD.TXT', header=T, sep='\t')

# Create an array of the months in chronological order
monthList = c('Jan','Feb','Mar','Apr','May','Jun',
              'Jul','Aug','Sep','Oct','Nov','Dec')

# Put rainfall data in long format
rainfall_long = pivot_longer(data=rainfall, 
                             names_to = 'Month', 
                             values_to = 'rain', 
                             cols = any_of(monthList))

# Define the order of the months
rainfall_long$Month = ordered(rainfall_long$Month, 
                              levels=monthList)

# Create the heatmap plot
ggplot(data=rainfall_long, 
       aes(x=Year,
           y=Month, 
           fill=rain)) + 
  geom_tile(colour='white') +
  scale_fill_gradient(low = "white", 
                      high = "steelblue") +
  labs(x="Year", 
       y="Month", 
       fill = "Rainfall (mm)", 
       title="Monthly Rainfall at Malin Head") +
  theme_bw() +
  theme(legend.title = element_text(size = 10),
        legend.text = element_text(size = 12),
        plot.title = element_text(size=16),
        axis.title=element_text(size=14, face="bold"),
        axis.text = element_text(size=14),
        axis.text.x = element_text(angle = 90, hjust = 1)) 

# Display a correlation plot using the corrplot package

library(corrplot)

# First calculate correlation coefficients to be visualised
cor_matrix = cor(msleep[ ,c(6:9)], method='pearson',use='complete.obs')

# Now produce the plot
corrplot(cor_matrix, method='circle', type='upper', addCoef.col = "black")