Key features

Can be used to easily create an interactive sankey diagram, as well as, other network layout such as dendrogram, radial and diagnonal networks.

Key R functions and options

Key R functions:

forceNetwork(). Creates a D3 JavaScript force directed network graph

forceNetwork(Links, Nodes, Source, Target,
             Value, NodeID, Nodesize, Group)

Key Arguments:

• Links: edges list. Edge IDs should start with 0
• Nodes: Nodes list. Node IDs should start with 0
• Source, Target: the names of the column, in the edges data, containing the network source and target variables, respectively.
• Value: the name of the column, in the edge data, containing the weight values for edges. Used to indicate how wide the links are.
• NodeID: the name of the column, in the nodes data, containing the node IDs. Used for labeling the nodes.
• Nodesize: the name of the column, in the nodes data, with some value to vary the node radius’s with.
• Group: the name of the column, in the nodes data, specifying the group of each node.

Prepare nodes and edes data

1. Prepare the nodes and the edges data:

nodes_d3 <- mutate(nodes, id = id - 1)
edges_d3 <- mutate(edges, from = from - 1, to = to - 1)

2. Create the interactive network:

library(networkD3)

forceNetwork(
  Links = edges_d3, Nodes = nodes_d3,  
  Source = "from", Target = "to",      # so the network is directed.
  NodeID = "label", Group = "id", Value = "weight", 
  opacity = 1, fontSize = 16, zoom = TRUE
  )

Create sankey diagram

You can create a d3-styled sankey diagram. A Sankey diagram is a good fit for the phone call data. There are not too many nodes in the data, making it easier to visualize the flow of phone calls.

Create a sankey diagram:

sankeyNetwork(
  Links = edges_d3, Nodes = nodes_d3, 
  Source = "from", Target = "to", 
  NodeID = "label", Value = "weight", 
  fontSize = 16, unit = "Letter(s)")

Other hierarchical layouts exist in the network3D package to visualize tree-like graphs. In the example below, we start by computing hierarchical clustering using a sample of the USArrests data set:

set.seed(123)
hc <- USArrests %>% sample_n(15) %>%
  scale() %>% dist() %>%
  hclust(method = "complete")

Other network layouts

• dendroNetwork:

dendroNetwork(hc, fontSize = 15)

Other alternatives are:

• radialNetwork:

radialNetwork(
  as.radialNetwork(hc), fontSize = 15
  )

• diagonalNetwork:

diagonalNetwork(
  as.radialNetwork(hc), fontSize = 15
  )

Key features

• Creates interactive network graphs.
• Possible to customize nodes and edge as you want.
• Can be used to directly visualize interactively a network generated with the igraph package.
• Can be used to visualize recursive partitioning and regression trees generated with the rpart package.
• Possible to use images and icons for node shapes.
• Supports igraph layouts

Key R function and options

Key R function:

visNetwork( 
  nodes = NULL, edges = NULL,
  width = NULL, height = NULL, 
  main = NULL, submain = NULL, footer = NULL
  )

Key Arguments:

• nodes: nodes list information. Should contain at least the column “id”. See visNodes() for more options to control nodes. Other colums can be included in the data, such as:
  • “id” : id of the node, needed in edges information
  • “label” : label of the node
  • “group” : group of the node. Groups can be configure with visGroups().
  • “value” : size of the node
  • “title” : tooltip of the node
  • …
• edges: edges list information. Required at least columns “from” and “to”. See visEdges() for more options to control edges.
  • “from” : node id of begin of the edge
  • “to” : node id of end of the edge
  • “label” : label of the edge
  • “value” : size of the node
  • “title” : tooltip of the node
  • …

Create a classic network graphs

To have always the same network, you can use the function visLayout(randomSeed = 12):

library("visNetwork")

visNetwork(nodes, edges) %>%
  visLayout(randomSeed = 12) 

If you want to control the width of edges according to a variable, you should include the column “width” in the edge list data. You should manually calculate and scale the edge width.

In the following R code, we’ll customize the visNetwork() output by using an igraph layout and changing the edges width.

First add the column width in the edges list data frame. Set the minimum width to 1:

edges <- mutate(edges, width = 1 + weight/5)

Create the network graph with the variable edge widths and the igraph layout = “layout_with_fr”.

visNetwork(nodes, edges) %>% 
  visIgraphLayout(layout = "layout_with_fr") %>% 
  visEdges(arrows = "middle") %>%
  visLayout(randomSeed = 1234)  

Visualize classification and regression trees

As mentioned above, you can visualize classification and regression trees generated using the rpart package.

Key function: visTree() [in visNetwork version >= 2.0.0].

For example, to visualize a classification tree, type the following R code:

# Compute
library(rpart)
res <- rpart(Species~., data=iris)
# Visualize
visTree(res, main = "Iris classification Tree",
        width = "80%",  height = "400px")

Use tbl_graph

• Create a tbl_graph network object using the phone call data:

library("navdata")
data("phone.call2")
phone.net <- tbl_graph(
  nodes = phone.call2$nodes, 
  edges = phone.call2$edges,
  directed = TRUE
  )

• Visualize:

ggraph(phone.net, layout = "graphopt") + 
  geom_edge_link(width = 1, colour = "lightgray") +
  geom_node_point(size = 4, colour = "#00AFBB") +
  geom_node_text(aes(label = label), repel = TRUE)+
  theme_graph()

Use as_tbl_graph R function

One can also use the as_tbl_graph() function to converts the following data structure and network objects:

• data.frame, list and matrix data [R base]
• igraph network objects [igraph package]
• network network objects [network pakage]
• dendrogram and hclust [stats package]
• Node [data.tree package]
• phylo and evonet [ape package]
• graphNEL, graphAM, graphBAM [graph package (in Bioconductor)]

In the following example, we’ll create a correlation matrix network graph. The mtcars data set will be used.

Compute the correlation matrix between cars using the corrr package:

• 1. Use the mtcars data set
• 2. Compute the correlation matrix: correlate()
• 3. Convert the upper triangle to NA: shave()
• 4. Stretch the correlation data frame into long format
• 5. Keep only high correlation

library(corrr)
res.cor <- mtcars [, c(1, 3:6)] %>%  # (1)
  t() %>% correlate() %>%            # (2)
  shave(upper = TRUE) %>%            # (3)
  stretch(na.rm = TRUE) %>%          # (4)
  filter(r >= 0.998)                 # (5)
res.cor

## # A tibble: 59 x 3
##               x             y     r
##                     
## 1     Mazda RX4 Mazda RX4 Wag 1.000
## 2     Mazda RX4      Merc 230 1.000
## 3     Mazda RX4      Merc 280 0.999
## 4     Mazda RX4     Merc 280C 0.999
## 5     Mazda RX4    Merc 450SL 0.998
## 6 Mazda RX4 Wag      Merc 230 1.000
## # ... with 53 more rows

Create the correlation network graph:

set.seed(1)
cor.graph <- as_tbl_graph(res.cor, directed = FALSE)
ggraph(cor.graph) + 
  geom_edge_link() + 
  geom_node_point() +
  geom_node_text(
    aes(label = name), size = 3, repel = TRUE
    ) +
  theme_graph()

Centrality

Centrality is an important concept when analyzing network graph. The centrality of a node / edge measures how central (or important) is a node or edge in the network.

We consider an entity important, if he has connections to many other entities. Centrality describes the number of edges that are connected to nodes.

There many types of scores that determine centrality. One of the famous ones is the pagerank algorithm that was powering Google Search in the beginning.

Examples of common approaches of measuring centrality include:

• betweenness centrality. The betweenness centrality for each nodes is the number of the shortest paths that pass through the nodes.

• closeness centrality. Closeness centrality measures how many steps is required to access every other nodes from a given nodes. It describes the distance of a node to all other nodes. The more central a node is, the closer it is to all other nodes.

• eigenvector centrality. A node is important if it is linked to by other important nodes. The centrality of each node is proportional to the sum of the centralities of those nodes to which it is connected. In general, nodes with high eigenvector centralities are those which are linked to many other nodes which are, in turn, connected to many others (and so on).

• Hub and authority centarlities are generalization of eigenvector centrality. A high hub node points to many good authorities and a high authority node receives from many good hubs.

The tidygraph package contains more than 10 centrality measures, prefixed with the term centrality_. These measures include:

centrality_authority()

centrality_betweenness()

centrality_closeness()

centrality_hub()

centrality_pagerank()

centrality_eigen()

centrality_edge_betweenness()

All of these centrality functions returns a numeric vector matching the nodes (or edges in the case of `centrality_edge_betweenness()).

In the following examples, we’ll use the phone call network graph. We’ll change the color and the size of nodes according to their values of centrality.

set.seed(123)
phone.net %>%
  activate(nodes) %>%
  mutate(centrality = centrality_authority()) %>% 
  ggraph(layout = "graphopt") + 
  geom_edge_link(width = 1, colour = "lightgray") +
  geom_node_point(aes(size = centrality, colour = centrality)) +
  geom_node_text(aes(label = label), repel = TRUE)+
  scale_color_gradient(low = "yellow", high = "red")+
  theme_graph()

Clustering

Clustering is a common operation in network analysis and it consists of grouping nodes based on the graph topology.

It’s sometimes referred to as community detection based on its commonality in social network analysis.

Many clustering algorithms from are available in the tidygraph package and prefixed with the term group_. These include:

• Infomap community finding. It groups nodes by minimizing the expected description length of a random walker trajectory. R function: group_infomap()
• Community structure detection based on edge betweenness. It groups densely connected nodes. R function: group_edge_betweenness().

In the following example, we’ll use the correlation network graphs to detect clusters or communities:

set.seed(123)
cor.graph %>%
  activate(nodes) %>%
   mutate(community = as.factor(group_infomap())) %>% 
  ggraph(layout = "graphopt") + 
  geom_edge_link(width = 1, colour = "lightgray") +
  geom_node_point(aes(colour = community), size = 4) +
  geom_node_text(aes(label = label), repel = TRUE)+
  theme_graph()

Demo data set

We’ll use a fake demo data set containing the number of phone calls between the president of some EU countries.

library("navdata")
data("phone.call")
head(phone.call, 3)

## # A tibble: 3 x 3
##    source destination n.call
##              
## 1  France     Germany      9
## 2 Belgium      France      4
## 3  France       Spain      3

To visualize the network graph, we need to create two data frames from the demo data sets:

• nodes list: containing nodes labels and other nodes attributes
• edges list: containing the relationship between the nodes. It consists of the edge list and any additional edge attributes.

In the following sections, we start by creating nodes and edges lists. Next, we’ll use the different packages to create network graphs.

Create nodes list

First, load the tidyverse R package for data manipulation:

library(tidyverse)

Then, compute the following key steps to create nodes list:

1. Take the distinct countries from both the “source” and “destination” columns
2. Change the column name to label
3. Join the information from the two columns together.

#  Get distinct source names
sources <- phone.call %>%
  distinct(source) %>%
  rename(label = source)

# Get distinct destination names
destinations <- phone.call %>%
  distinct(destination) %>%
  rename(label = destination)

# Join the two data to create node
# Add unique ID for each country
nodes <- full_join(sources, destinations, by = "label") 
nodes <- nodes %>%
  mutate(id = 1:nrow(nodes)) %>%
  select(id, everything())

head(nodes, 3)

## # A tibble: 3 x 2
##      id   label
##      
## 1     1  France
## 2     2 Belgium
## 3     3 Germany

Create edges list

Key steps:

1. Take the phone.call data, which are already in edges list format, showing the connection between nodes. Rename the column “n.call” to “weight”.
2. Join the node IDs to the edges list data
   1. Do this for the “source” column and rename the id column that are brought over from nodes. New name: “from”.
   2. Do this for the “destination” column and rename the id column. New name: “to”
   3. Select only the columns “from” and “to” in the edge data. We don’t need to keep the column “source” and “destination” containing the names of countries. These information are already present in the node data.

# Rename the n.call column to weight
phone.call <- phone.call %>%
  rename(weight = n.call)

# (a) Join nodes id for source column
edges <- phone.call %>% 
  left_join(nodes, by = c("source" = "label")) %>% 
  rename(from = id)

# (b) Join nodes id for destination column
edges <- edges %>% 
  left_join(nodes, by = c("destination" = "label")) %>% 
  rename(to = id)

# (c) Select/keep only the columns from and to
edges <- select(edges, from, to, weight)
head(edges, 3)

## # A tibble: 3 x 3
##    from    to weight
##      
## 1     1     3      9
## 2     2     1      4
## 3     1     8      3

igraph

1. Create an igraph network object:

• Key R function: graph_from_data_frame().

• Key arguments:
  • d: edge list
  • vertices: node list
  • directed: can be either TRUE or FALSE depending on whether the data is directed or undirected.

library(igraph)
net.igraph <- graph_from_data_frame(
  d = edges, vertices = nodes, 
  directed = TRUE
  )

2. Create a network graph with igraph

set.seed(123)
plot(net.igraph, edge.arrow.size = 0.2,
     layout = layout_with_graphopt)

tidygraph and ggraph

tidygraph and ggraph are modern R packages for network data manipulation (tidygraph) and visualization (ggraph). They leverage the power of igraph.

1. Create a network object using tidygraph:

• Key function: tbl_graph().
• key arguments: nodes, edges and directed.

library(tidygraph)
net.tidy <- tbl_graph(
  nodes = nodes, edges = edges, directed = TRUE
  )

2. Visualize network using ggraph

Key functions:

• geom_node_point(): Draws node points.

• geom_edge_link(): Draws edge links. To control the width of edge line according to the weight variable, specify the option aes(width = weight), where, the weight specify the number of phone.call sent along each route. In this case, you can control the maximum and minimum width of the edges, by using the function scale_edge_width() to set the range (minimum and maximum width value). For example: scale_edge_width(range = c(0.2, 2)).

• geom_node_text(): Adds text labels for nodes, by specifying the argument aes(label = label). To avoid text overlapping, indicate the option repel = TRUE.

• labs(): Change main titles, axis labels and legend titles.

Create a classic node-edge diagrams. Possible values for the argument layout include: 'star', 'circle', 'gem', 'dh', 'graphopt', 'grid', 'mds', 'randomly', 'fr', 'kk', 'drl', 'lgl'.

library(ggraph)
ggraph(net.tidy, layout = "graphopt") + 
  geom_node_point() +
  geom_edge_link(aes(width = weight), alpha = 0.8) + 
  scale_edge_width(range = c(0.2, 2)) +
  geom_node_text(aes(label = label), repel = TRUE) +
  labs(edge_width = "phone.call") +
  theme_graph()

Arc diagram layout

In the following example, we’ll:

• Layout the nodes linearly (horizontal line) using layout = "linear".
• Create an arc diagram by drawing the edges as arcs
• Add only the label names, instead of including node points.

# Arc diagram
ggraph(net.tidy, layout = "linear") + 
  geom_edge_arc(aes(width = weight), alpha = 0.8) + 
  scale_edge_width(range = c(0.2, 2)) +
  geom_node_text(aes(label = label), repel = TRUE) +
  labs(edge_width = "Number of calls") +
  theme_graph()+
  theme(legend.position = "top") 

# Coord diagram, circular
ggraph(net.tidy, layout = "linear", circular = TRUE) + 
  geom_edge_arc(aes(width = weight), alpha = 0.8) + 
  scale_edge_width(range = c(0.2, 2)) +
  geom_node_text(aes(label = label), repel = TRUE) +
  labs(edge_width = "Number of calls") +
  theme_graph()+
  theme(legend.position = "top")

Treemap layout

A treemap is a visual method for displaying hierarchical data that uses nested rectangles to represent the branches of a tree diagram. Each rectangles has an area proportional to the amount of data it represents.

To illustrate this layout, we’ll use the france.trade demo data set [ in navdata package]. It contains the trading percentage between France and different countries.

data("france.trade")
france.trade

## # A tibble: 251 x 3
##   source destination trade.percentage
##                       
## 1 France       Aruba            0.035
## 2 France Afghanistan            0.035
## 3 France      Angola            0.035
## 4 France    Anguilla            0.035
## 5 France     Albania            0.035
## 6 France     Finland            0.035
## # ... with 245 more rows

# Distinct countries
countries <- c(
  france.trade$source, france.trade$destination
) %>%
  unique()
# Create nodes list
nodes <- data_frame(
  id = 1:length(countries),
  label = countries
)

nodes <- nodes %>%
  left_join(
    france.trade[, c("destination", "trade.percentage")],
    by = c("label" = "destination" )
    ) %>%
  mutate(
    trade.percentage = ifelse(
      is.na(trade.percentage), 0, trade.percentage
      )
  )
head(nodes, 3)

## # A tibble: 3 x 3
##      id       label trade.percentage
##                      
## 1     1      France            0.000
## 2     2       Aruba            0.035
## 3     3 Afghanistan            0.035

per_route <- france.trade %>%  
  select(source, destination)

# (a) Join nodes id for source column
edges <- per_route %>% 
  left_join(nodes, by = c("source" = "label")) %>% 
  rename(from = id)

# (b) Join nodes id for destination column
edges <- edges %>% 
  left_join(nodes, by = c("destination" = "label")) %>% 
  rename(to = id)

# (c) Select/keep only the columns from and to
edges <- select(edges, from, to)
head(edges, 3)

## # A tibble: 3 x 2
##    from    to
##    
## 1     1     2
## 2     1     3
## 3     1     4

trade.graph <- tbl_graph(
  nodes = nodes, edges = edges, directed = TRUE
  )

# Generate colors for each country
require(ggpubr)
cols <- get_palette("Dark2", nrow(france.trade)+1)

# Visualize
set.seed(123)
ggraph(trade.graph, &#39;treemap&#39;, weight = "trade.percentage") + 
    geom_node_tile(aes(fill = label), size = 0.25, color = "white")+
  geom_node_text(
    aes(label = paste(label, trade.percentage, sep = "\n"),
        size = trade.percentage), color = "white"
    )+
  scale_fill_manual(values = cols)+
  scale_size(range = c(0, 6) )+
  theme_void()+
  theme(legend.position = "none")

require("map")
# (1)
world_map <- map_data("world")
france.trade <- france.trade %>%
  left_join(world_map, by = c("destination" = "region")) 
# (2)
ggplot(france.trade, aes(long, lat, group = group))+
  geom_polygon(aes(fill = trade.percentage ), color = "white")+
  geom_polygon(data = subset(world_map, region == "France"),
               fill = "red")+ # draw france in red
  scale_fill_gradientn(colours =c("lightgray", "yellow", "green"))+
  theme_void()

Dendrogram layout

Dendrogram layout are suited for hierarchical graph visualization, that is graph structures including trees and hierarchies.

In this section, we’ll compute hierarchical clustering using the USArrests data set. The output is visualized as a dendrogram tree.

1. compute hierarchical clustering using the USArrests data set;
2. then convert the result into a tbl_graph.

res.hclust <- scale(USArrests) %>%
  dist() %>% hclust() 
res.tree <- as.dendrogram(res.hclust)

3. Visualize the dendrogram tree. Key function: geom_edge_diagonal() and geom_edge_elbow()

# Diagonal layout
ggraph(res.tree, layout = "dendrogram") + 
  geom_edge_diagonal() +
  geom_node_text(
    aes(label = label), angle = 90, hjust = 1,
    size = 3
    )+
  ylim(-1.5, NA)+
  theme_minimal()

# Elbow layout
ggraph(res.tree, layout = "dendrogram") + 
  geom_edge_elbow() +
  geom_node_text(
    aes(label = label), angle = 90, hjust = 1,
    size = 3
    )+
  ylim(-1.5, NA)+theme_minimal()

R base graphs

R comes with simple functions to create many types of graphs. For example:

In the most cases, you can use the following arguments to customize the plot:

• pch: change point shapes. Allowed values comprise number from 1 to 25.
• cex: change point size. Example: cex = 0.8.
• col: change point color. Example: col = “blue”.
• frame: logical value. frame = FALSE removes the plot panel border frame.
• main, xlab, ylab. Specify the main title and the x/y axis labels -, respectively
• las: For a vertical x axis text, use las = 2.

In the following R code, we’ll use the iris data set to create a:

• 1. Scatter plot of Sepal.Length (on x-axis) and Sepal.Width (on y-axis).
• 2. Box plot of Sepal.length (y-axis) by Species (x-axis)

# (1) Create a scatter lot
plot(
  x = iris$Sepal.Length, y = iris$Sepal.Width,
  pch = 19, cex = 0.8, frame = FALSE,
  xlab = "Sepal Length",ylab = "Sepal Width"
  )
# (2) Create a box plot
boxplot(Sepal.Length ~ Species, data = iris,
        ylab = "Sepal.Length", 
        frame = FALSE, col = "lightgray")

Lattice graphics

The lattice R package provides a plotting system that aims to improve on R base graphs. After installing the package, whith the R command install.packages("lattice"), you can test the following functions.

• Main functions in the lattice package:

• Create a basic scatter plot of y by x. Syntax: y ~ x. Change the color by groups and use auto.key = TRUE to show legends:

library("lattice")
xyplot(
  Sepal.Length ~ Petal.Length, group = Species, 
  data = iris, auto.key = TRUE, pch = 19, cex = 0.5
  )

• Multiple panel plots by groups. Syntax: y ~ x | group.

xyplot(
  Sepal.Length ~ Petal.Length | Species, 
  layout = c(3, 1),               # panel with ncol = 3 and nrow = 1
  group = Species, data = iris,
  type = c("p", "smooth"),        # Show points and smoothed line
  scales = "free"                 # Make panels axis scales independent
  )

ggplot2 graphics

GGPlot2 is a powerful and a flexible R package, implemented by Hadley Wickham, for producing elegant graphics piece by piece. The gg in ggplot2 means Grammar of Graphics, a graphic concept which describes plots by using a “grammar”. According to the ggplot2 concept, a plot can be divided into different fundamental parts: Plot = data + Aesthetics + Geometry

• data: a data frame
• aesthetics: used to indicate the x and y variables. It can be also used to control the color, the size and the shape of points, etc…..
• geometry: corresponds to the type of graphics (histogram, box plot, line plot, ….)

After installing and loading the ggplot2 package, you can use the following key functions:

The main function in the ggplot2 package is ggplot(), which can be used to initialize the plotting system with data and x/y variables.

For example, the following R code takes the iris data set to initialize the ggplot and then a layer (geom_point()) is added onto the ggplot to create a scatter plot of x = Sepal.Length by y = Sepal.Width:

library(ggplot2)
ggplot(iris, aes(x = Sepal.Length, y = Sepal.Width))+
  geom_point()
# Change point size, color and shape
ggplot(iris, aes(x = Sepal.Length, y = Sepal.Width))+
  geom_point(size = 1.2, color = "steelblue", shape = 21)

Note that, in the code above, the shape of points is specified as number. To display the different point shape available in R, type this:

ggpubr::show_point_shapes()

It’s also possible to control points shape and color by a grouping variable (here, Species). For example, in the code below, we map points color and shape to the Species grouping variable.

# Control points color by groups
ggplot(iris, aes(x = Sepal.Length, y = Sepal.Width))+
  geom_point(aes(color = Species, shape = Species))
# Change the default color manually.
# Use the scale_color_manual() function
ggplot(iris, aes(x = Sepal.Length, y = Sepal.Width))+
  geom_point(aes(color = Species, shape = Species))+
  scale_color_manual(values = c("#00AFBB", "#E7B800", "#FC4E07"))

You can also split the plot into multiple panels according to a grouping variable. R function: facet_wrap(). Another interesting feature of ggplot2, is the possibility to combine multiple layers on the same plot. For example, with the following R code, we’ll:

• Add points with geom_point(), colored by groups.
• Add the fitted smoothed regression line using geom_smooth(). By default the function geom_smooth() add the regression line and the confidence area. You can control the line color and confidence area fill color by groups.
• Facet the plot into multiple panels by groups
• Change color and fill manually using the function scale_color_manual() and scale_fill_manual()

ggplot(iris, aes(x = Sepal.Length, y = Sepal.Width))+
  geom_point(aes(color = Species))+               
  geom_smooth(aes(color = Species, fill = Species))+
  facet_wrap(~Species, ncol = 3, nrow = 1)+
  scale_color_manual(values = c("#00AFBB", "#E7B800", "#FC4E07"))+
  scale_fill_manual(values = c("#00AFBB", "#E7B800", "#FC4E07"))

Note that, the default theme of ggplots is theme_gray() (or theme_grey()), which is theme with grey background and white grid lines. More themes are available for professional presentations or publications. These include: theme_bw(), theme_classic() and theme_minimal().

To change the theme of a given ggplot (p), use this: p + theme_classic(). To change the default theme to theme_classic() for all the future ggplots during your entire R session, type the following R code:

theme_set(
  theme_classic()
)

Now you can create ggplots with theme_classic() as default theme:

ggplot(iris, aes(x = Sepal.Length, y = Sepal.Width))+
  geom_point()

ggpubr for publication ready plots

The ggpubr R package facilitates the creation of beautiful ggplot2-based graphs for researcher with non-advanced programming backgrounds (Kassambara 2017).

For example, to create the density distribution of “Sepal.Length”, colored by groups (“Species”), type this:

library(ggpubr)
# Density plot with mean lines and marginal rug
ggdensity(iris, x = "Sepal.Length",
   add = "mean", rug = TRUE,             # Add mean line and marginal rugs
   color = "Species", fill = "Species",  # Color by groups
   palette = "jco")                      # use jco journal color palette

• Create a box plot with p-values comparing groups:

# Groups that we want to compare
my_comparisons <- list(
  c("setosa", "versicolor"), c("versicolor", "virginica"),
  c("setosa", "virginica")
)
# Create the box plot. Change colors by groups: Species
# Add jitter points and change the shape by groups
ggboxplot(
  iris, x = "Species", y = "Sepal.Length",
  color = "Species", palette = c("#00AFBB", "#E7B800", "#FC4E07"),
  add = "jitter"
  )+
  stat_compare_means(comparisons = my_comparisons, method = "t.test")

Bar plot of counts

• Plot types: Bar plot of the count of group levels
• Key function: geom_bar()
• Key arguments: alpha, color, fill, linetype and size

Demo data set: diamonds [in ggplot2]. Contains the prices and other attributes of almost 54000 diamonds. The column cut contains the quality of the diamonds cut (Fair, Good, Very Good, Premium, Ideal).

The R code below creates a bar plot visualizing the number of elements in each category of diamonds cut.

ggplot(diamonds, aes(cut)) +
  geom_bar(fill = "#0073C2FF") +
  theme_pubclean()

Compute the frequency of each category and add the labels on the bar plot:

• dplyr package used to summarise the data
• geom_bar() with option stat = "identity" is used to create the bar plot of the summary output as it is.
• geom_text() used to add text labels. Adjust the position of the labels by using hjust (horizontal justification) and vjust (vertical justification). Values should be in [0, 1].

# Compute the frequency
library(dplyr)
df <- diamonds %>%
  group_by(cut) %>%
  summarise(counts = n())
df

## # A tibble: 5 x 2
##         cut counts
##         
## 1      Fair   1610
## 2      Good   4906
## 3 Very Good  12082
## 4   Premium  13791
## 5     Ideal  21551

# Create the bar plot. Use theme_pubclean() [in ggpubr]
ggplot(df, aes(x = cut, y = counts)) +
  geom_bar(fill = "#0073C2FF", stat = "identity") +
  geom_text(aes(label = counts), vjust = -0.3) + 
  theme_pubclean()

Pie charts

Pie chart is just a stacked bar chart in polar coordinates.

First,

• Arrange the grouping variable (cut) in descending order. This important to compute the y coordinates of labels.
• compute the proportion (counts/total) of each category
• compute the position of the text labels as the cumulative sum of the proportion. To put the labels in the center of pies, we’ll use cumsum(prop) - 0.5*prop as label position.

df <- df %>%
  arrange(desc(cut)) %>%
  mutate(prop = round(counts*100/sum(counts), 1),
         lab.ypos = cumsum(prop) - 0.5*prop)
head(df, 4)

## # A tibble: 4 x 4
##         cut counts  prop lab.ypos
##              
## 1     Ideal  21551  40.0     20.0
## 2   Premium  13791  25.6     52.8
## 3 Very Good  12082  22.4     76.8
## 4      Good   4906   9.1     92.5

• Create the pie charts using ggplot2 verbs. Key function: coord_polar().

ggplot(df, aes(x = "", y = prop, fill = cut)) +
  geom_bar(width = 1, stat = "identity", color = "white") +
  geom_text(aes(y = lab.ypos, label = prop), color = "white")+
  coord_polar("y", start = 0)+
  ggpubr::fill_palette("jco")+
  theme_void()

• Alternative solution to easily create a pie chart: use the function ggpie()[in ggpubr]:

ggpie(
  df, x = "prop", label = "prop",
  lab.pos = "in", lab.font = list(color = "white"), 
  fill = "cut", color = "white",
  palette = "jco"
)

Dot charts

Dot chart is an alternative to bar plots. Key functions:

• geom_linerange():Creates line segments from x to ymax
• geom_point(): adds dots
• ggpubr::color_palette(): changes color palette.

ggplot(df, aes(cut, prop)) +
  geom_linerange(
    aes(x = cut, ymin = 0, ymax = prop), 
    color = "lightgray", size = 1.5
    )+
  geom_point(aes(color = cut), size = 2)+
  ggpubr::color_palette("jco")+
  theme_pubclean()

Easy alternative to create a dot chart. Use ggdotchart() [ggpubr]:

ggdotchart(
  df, x = "cut", y = "prop",
  color = "cut", size = 3,      # Points color and size
  add = "segment",              # Add line segments
  add.params = list(size = 2), 
  palette = "jco",
  ggtheme = theme_pubclean()
)

Data format

Create some data (wdata) containing the weights by sex (M for male; F for female):

set.seed(1234)
wdata = data.frame(
        sex = factor(rep(c("F", "M"), each=200)),
        weight = c(rnorm(200, 55), rnorm(200, 58))
        )

head(wdata, 4)

##   sex weight
## 1   F   53.8
## 2   F   55.3
## 3   F   56.1
## 4   F   52.7

Compute the mean weight by sex using the dplyr package. First, the data is grouped by sex and then summarized by computing the mean weight by groups. The operator %>% is used to combine multiple operations:

library("dplyr")
mu <- wdata %>% 
  group_by(sex) %>%
  summarise(grp.mean = mean(weight))
mu

## # A tibble: 2 x 2
##      sex grp.mean
##       
## 1      F     54.9
## 2      M     58.1

Basic plots

We start by creating a plot, named a, that we’ll finish in the next section by adding a layer.

a <- ggplot(wdata, aes(x = weight))

Possible layers include: geom_density() (for density plots) and geom_histogram() (for histogram plots).

Key arguments to customize the plots:

• color, size, linetype: change the line color, size and type, respectively
• fill: change the areas fill color (for bar plots, histograms and density plots)
• alpha: create a semi-transparent color.

Density plots

Key function: geom_density().

1. Create basic density plots. Add a vertical line corresponding to the mean value of the weight variable (geom_vline()):

# y axis scale = ..density.. (default behaviour)
a + geom_density() +
  geom_vline(aes(xintercept = mean(weight)), 
             linetype = "dashed", size = 0.6)
  
# Change y axis to count instead of density
a + geom_density(aes(y = ..count..), fill = "lightgray") +
  geom_vline(aes(xintercept = mean(weight)), 
             linetype = "dashed", size = 0.6,
             color = "#FC4E07")

2. Change areas fill and add line color by groups (sex):

• Add vertical mean lines using geom_vline(). Data: mu, which contains the mean values of weights by sex (computed in the previous section).
• Change color manually:
  • use scale_color_manual() or scale_colour_manual() for changing line color
  • use scale_fill_manual() for changing area fill colors.

# Change line color by sex
a + geom_density(aes(color = sex)) +
  scale_color_manual(values = c("#868686FF", "#EFC000FF"))

# Change fill color by sex and add mean line
# Use semi-transparent fill: alpha = 0.4
a + geom_density(aes(fill = sex), alpha = 0.4) +
      geom_vline(aes(xintercept = grp.mean, color = sex),
             data = mu, linetype = "dashed") +
  scale_color_manual(values = c("#868686FF", "#EFC000FF"))+
  scale_fill_manual(values = c("#868686FF", "#EFC000FF"))

3. Simple solution to create a ggplot2-based density plots: use ggboxplot() [in ggpubr].

library(ggpubr)

# Basic density plot with mean line and marginal rug
ggdensity(wdata, x = "weight", 
          fill = "#0073C2FF", color = "#0073C2FF",
          add = "mean", rug = TRUE)
     
# Change outline and fill colors by groups ("sex")
# Use a custom palette
ggdensity(wdata, x = "weight",
   add = "mean", rug = TRUE,
   color = "sex", fill = "sex",
   palette = c("#0073C2FF", "#FC4E07"))

Histogram plots

An alternative to density plots is histograms, which represents the distribution of a continuous variable by dividing into bins and counting the number of observations in each bin.

Key function: geom_histogram(). The basic usage is quite similar to geom_density().

1. Create a basic plots. Add a vertical line corresponding to the mean value of the weight variable:

a + geom_histogram(bins = 30, color = "black", fill = "gray") +
  geom_vline(aes(xintercept = mean(weight)), 
             linetype = "dashed", size = 0.6)

2. Change areas fill and add line color by groups (sex):

• Add vertical mean lines using geom_vline(). Data: mu, which contains the mean values of weights by sex.
• Change color manually:
  • use scale_color_manual() or scale_colour_manual() for changing line color
  • use scale_fill_manual() for changing area fill colors.
• Adjust the position of histogram bars by using the argument position. Allowed values: “identity”, “stack”, “dodge”. Default value is “stack”.

# Change line color by sex
a + geom_histogram(aes(color = sex), fill = "white", 
                   position = "identity") +
  scale_color_manual(values = c("#00AFBB", "#E7B800")) 

# change fill and outline color manually 
a + geom_histogram(aes(color = sex, fill = sex),
                         alpha = 0.4, position = "identity") +
  scale_fill_manual(values = c("#00AFBB", "#E7B800")) +
  scale_color_manual(values = c("#00AFBB", "#E7B800"))

3. Combine histogram and density plots:

• Plot histogram with density values on y-axis (instead of count values).
• Add density plot with transparent density plot

# Histogram with density plot
a + geom_histogram(aes(y = ..density..), 
                   colour="black", fill="white") +
  geom_density(alpha = 0.2, fill = "#FF6666") 
     

# Color by groups
a + geom_histogram(aes(y = ..density.., color = sex), 
                   fill = "white",
                   position = "identity")+
  geom_density(aes(color = sex), size = 1) +
  scale_color_manual(values = c("#868686FF", "#EFC000FF"))

4. Simple solution to create a ggplot2-based histogram plots: use gghistogram() [in ggpubr].

library(ggpubr)

# Basic histogram plot with mean line and marginal rug
gghistogram(wdata, x = "weight", bins = 30, 
            fill = "#0073C2FF", color = "#0073C2FF",
            add = "mean", rug = TRUE)
     
# Change outline and fill colors by groups ("sex")
# Use a custom palette
gghistogram(wdata, x = "weight", bins = 30,
   add = "mean", rug = TRUE,
   color = "sex", fill = "sex",
   palette = c("#0073C2FF", "#FC4E07"))

Alternative to density and histogram plots

1. Frequency polygon. Very close to histogram plots, but it uses lines instead of bars.
   • Key function: geom_freqpoly().
   • Key arguments: color, size, linetype: change, respectively, line color, size and type.
2. Area plots. This is a continuous analog of a stacked bar plot.
   • Key function: geom_area().
   • Key arguments:
     • color, size, linetype: change, respectively, line color, size and type.
     • fill: change area fill color.

In this section, we’ll use the theme theme_pubclean() [in ggpubr]. This is a theme without axis lines, to direct more attention to the data. Type this to use the theme:

theme_set(theme_pubclean())

• Create a basic frequency polygon and basic area plots:

# Basic frequency polygon
a + geom_freqpoly(bins = 30) 

# Basic area plots, which can be filled by color
a + geom_area( stat = "bin", bins = 30,
               color = "black", fill = "#00AFBB")

• Change colors by groups (sex):

# Frequency polygon: 
# Change line colors and types by groups
a + geom_freqpoly( aes(color = sex, linetype = sex),
                   bins = 30, size = 1.5) +
  scale_color_manual(values = c("#00AFBB", "#E7B800"))

# Area plots: change fill colors by sex
# Create a stacked area plots
a + geom_area(aes(fill = sex), color = "white", 
              stat ="bin", bins = 30) +
  scale_fill_manual(values = c("#00AFBB", "#E7B800"))

3. Dot plots. Represents another alternative to histograms and density plots, that can be used to visualize a continuous variable. Dots are stacked with each dot representing one observation. The width of a dot corresponds to the bin width.

• Key function: geom_dotplot().
• Key arguments: alpha, color, fill and dotsize.

Create a dot plot colored by groups (sex):

a + geom_dotplot(aes(fill = sex), binwidth = 1/4) +
  scale_fill_manual(values = c("#00AFBB", "#E7B800"))

4. Box plot:
   • Create a box plot of one continuous variable: geom_boxplot()
   • Add jittered points, where each point corresponds to an individual observation: geom_jitter(). Change the color and the shape of points by groups (sex)

ggplot(wdata, aes(x = factor(1), y = weight)) +
  geom_boxplot(width = 0.4, fill = "white") +
  geom_jitter(aes(color = sex, shape = sex), 
              width = 0.1, size = 1) +
  scale_color_manual(values = c("#00AFBB", "#E7B800")) + 
  labs(x = NULL)   # Remove x axis label

5. Empirical cumulative distribution function (ECDF). Provides another alternative visualization of distribution. It reports for any given number the percent of individuals that are below that threshold.

For example, in the following plots, you can see that:

• about 25% of our females are shorter than 50 inches
• about 50% of males are shorter than 58 inches

# Another option for geom = "point"
a + stat_ecdf(aes(color = sex,linetype = sex), 
              geom = "step", size = 1.5) +
  scale_color_manual(values = c("#00AFBB", "#E7B800"))+
  labs(y = "f(weight)")

6. Quantile-quantile plot (QQ plots). Used to check whether a given data follows normal distribution.

• Key function: stat_qq().
• Key arguments: color, shape and size to change point color, shape and size.

Create a qq-plot of weight. Change color by groups (sex)

# Change point shapes by groups
ggplot(wdata, aes(sample = weight)) +
  stat_qq(aes(color = sex)) +
  scale_color_manual(values = c("#00AFBB", "#E7B800"))+
  labs(y = "Weight")

Alternative plot using the function ggqqplot() [in ggpubr]. The 95% confidence band is shown by default.

library(ggpubr)
ggqqplot(wdata, x = "weight",
   color = "sex", 
   palette = c("#0073C2FF", "#FC4E07"),
   ggtheme = theme_pubclean())

Density ridgeline plots

The density ridgeline plot is an alternative to the standard geom_density() function that can be useful for visualizing changes in distributions, of a continuous variable, over time or space. Ridgeline plots are partially overlapping line plots that create the impression of a mountain range.

This functionality is provided in the R package ggridges (Wilke 2017).

1. Installation:

install.packages("ggridges")

2. Load and set the default theme to theme_ridges() [in ggridges]:

library(ggplot2)
library(ggridges)
theme_set(theme_ridges())

3. Example 1: Simple distribution plots by groups. Distribution of Sepal.Length by Species using the iris data set. The grouping variable Species will be mapped to the y-axis:

ggplot(iris, aes(x = Sepal.Length, y = Species)) +
  geom_density_ridges(aes(fill = Species)) +
  scale_fill_manual(values = c("#00AFBB", "#E7B800", "#FC4E07"))

ggplot(iris, aes(x = Sepal.Length, y = Species)) +
  geom_density_ridges(scale = 0.9) 

4. Example 4: Visualize temperature data.

• Data set: lincoln_weather [in ggridges]. Weather in Lincoln, Nebraska in 2016.

• Create the density ridge plots of the Mean Temperature by Month and change the fill color according to the temperature value (on x axis). A gradient color is created using the function geom_density_ridges_gradient()

ggplot(
  lincoln_weather, 
  aes(x = `Mean Temperature [F]`, y = `Month`)
  ) +
  geom_density_ridges_gradient(
    aes(fill = ..x..), scale = 3, size = 0.3
    ) +
  scale_fill_gradientn(
    colours = c("#0D0887FF", "#CC4678FF", "#F0F921FF"),
    name = "Temp. [F]"
    )+
  labs(title = 'Temperatures in Lincoln NE') 

For more examples, type the following R code:

browseVignettes("ggridges")

Bar plot and modern alternatives

In this section, we’ll describe how to create easily basic and ordered bar plots using ggplot2 based helper functions available in the ggpubr R package. We’ll also present some modern alternatives to bar plots, including lollipop charts and cleveland’s dot plots.

• Load required packages:

library(ggpubr)

• Load and prepare data:

# Load data
dfm <- mtcars
# Convert the cyl variable to a factor
dfm$cyl <- as.factor(dfm$cyl)
# Add the name colums
dfm$name <- rownames(dfm)
# Inspect the data
head(dfm[, c("name", "wt", "mpg", "cyl")])

##                                name   wt  mpg cyl
## Mazda RX4                 Mazda RX4 2.62 21.0   6
## Mazda RX4 Wag         Mazda RX4 Wag 2.88 21.0   6
## Datsun 710               Datsun 710 2.32 22.8   4
## Hornet 4 Drive       Hornet 4 Drive 3.21 21.4   6
## Hornet Sportabout Hornet Sportabout 3.44 18.7   8
## Valiant                     Valiant 3.46 18.1   6

• Create an ordered bar plot of the mpg variable. Change the fill color by the grouping variable “cyl”. Sorting will be done globally, but not by groups.

ggbarplot(dfm, x = "name", y = "mpg",
          fill = "cyl",               # change fill color by cyl
          color = "white",            # Set bar border colors to white
          palette = "jco",            # jco journal color palett. see ?ggpar
          sort.val = "asc",          # Sort the value in dscending order
          sort.by.groups = TRUE,     # Don't sort inside each group
          x.text.angle = 90,           # Rotate vertically x axis texts
          ggtheme = theme_pubclean()
          )+
  font("x.text", size = 8, vjust = 0.5)

• Create a Lollipop chart:
  • Color by groups and set a custom color palette.
  • Sort values in ascending order.
  • Add segments from y = 0 to dots. Change segment color and size.

ggdotchart(dfm, x = "name", y = "mpg",
           color = "cyl",                                
           palette = c("#00AFBB", "#E7B800", "#FC4E07"), 
           sorting = "asc", sort.by.groups = TRUE,                      
           add = "segments",                            
           add.params = list(color = "lightgray", size = 2), 
           group = "cyl",                                
           dot.size = 4,                                 
           ggtheme = theme_pubclean()
           )+
  font("x.text", size = 8, vjust = 0.5)

Read more: Bar Plots and Modern Alternatives

Data format

• Demo dataset: ToothGrowth
  • Continuous variable: len (tooth length). Used on y-axis
  • Grouping variable: dose (dose levels of vitamin C: 0.5, 1, and 2 mg/day). Used on x-axis.

First, convert the variable dose from a numeric to a discrete factor variable:

data("ToothGrowth")
ToothGrowth$dose <- as.factor(ToothGrowth$dose)
head(ToothGrowth)

##    len supp dose
## 1  4.2   VC  0.5
## 2 11.5   VC  0.5
## 3  7.3   VC  0.5
## 4  5.8   VC  0.5
## 5  6.4   VC  0.5
## 6 10.0   VC  0.5

Box plots

• Key function: geom_boxplot()
• Key arguments to customize the plot:
  • width: the width of the box plot
  • notch: logical. If TRUE, creates a notched box plot. The notch displays a confidence interval around the median which is normally based on the median +/- 1.58*IQR/sqrt(n). Notches are used to compare groups; if the notches of two boxes do not overlap, this is a strong evidence that the medians differ.
  • color, size, linetype: Border line color, size and type
  • fill: box plot areas fill color
  • outlier.colour, outlier.shape, outlier.size: The color, the shape and the size for outlying points.

1. Create basic box plots:

• Standard and notched box plots:

# Default plot
e <- ggplot(ToothGrowth, aes(x = dose, y = len))
e + geom_boxplot()

# Notched box plot with mean points
e + geom_boxplot(notch = TRUE, fill = "lightgray")+
  stat_summary(fun.y = mean, geom = "point",
               shape = 18, size = 2.5, color = "#FC4E07")

• Change box plot colors by groups:

# Color by group (dose)
e + geom_boxplot(aes(color = dose))+
  scale_color_manual(values = c("#00AFBB", "#E7B800", "#FC4E07"))

# Change fill color by group (dose)
e + geom_boxplot(aes(fill = dose)) +
  scale_fill_manual(values = c("#00AFBB", "#E7B800", "#FC4E07"))

Note that, it’s possible to use the function scale_x_discrete() for:

• choosing which items to display: for example c(“0.5”, “2”),
• changing the order of items: for example from c(“0.5”, “1”, “2”) to c(“2”, “0.5”, “1”)

For example, type this:

# Choose which items to display: group "0.5" and "2"
e + geom_boxplot() + 
  scale_x_discrete(limits=c("0.5", "2"))

# Change the default order of items
e + geom_boxplot() +
  scale_x_discrete(limits=c("2", "0.5", "1"))

2. Create a box plot with multiple groups:

Two different grouping variables are used: dose on x-axis and supp as fill color (legend variable).

The space between the grouped box plots is adjusted using the function position_dodge().

e2 <- e + geom_boxplot(
  aes(fill = supp),
  position = position_dodge(0.9) 
  ) +
  scale_fill_manual(values = c("#999999", "#E69F00"))
e2

Split the plot into multiple panel. Use the function facet_wrap():

e2 + facet_wrap(~supp)

Violin plots

Violin plots are similar to box plots, except that they also show the kernel probability density of the data at different values. Typically, violin plots will include a marker for the median of the data and a box indicating the interquartile range, as in standard box plots.

Key function:

• geom_violin(): Creates violin plots. Key arguments:
  • color, size, linetype: Border line color, size and type
  • fill: Areas fill color
  • trim: logical value. If TRUE (default), trim the tails of the violins to the range of the data. If FALSE, don’t trim the tails.
• stat_summary(): Adds summary statistics (mean, median, …) on the violin plots.

1. Create basic violin plots with summary statistics:

# Add mean points +/- SD
# Use geom = "pointrange" or geom = "crossbar"
e + geom_violin(trim = FALSE) + 
  stat_summary(
    fun.data = "mean_sdl",  fun.args = list(mult = 1), 
    geom = "pointrange", color = "black"
    )
    
# Combine with box plot to add median and quartiles
# Change color by groups
e + geom_violin(aes(fill = dose), trim = FALSE) + 
  geom_boxplot(width = 0.2)+
  scale_fill_manual(values = c("#00AFBB", "#E7B800", "#FC4E07"))+
  theme(legend.position = "none")

2. Create violin plots with multiple groups:

e + geom_violin(
  aes(color = supp), trim = FALSE,
  position = position_dodge(0.9) 
  ) +
  geom_boxplot(
    aes(color = supp), width = 0.15,
    position = position_dodge(0.9)
    ) +
  scale_color_manual(values = c("#00AFBB", "#E7B800"))

Dot plots

• Key function: geom_dotplot(). Creates stacked dots, with each dot representing one observation.
• Key arguments:
  • stackdir: which direction to stack the dots. “up” (default), “down”, “center”, “centerwhole” (centered, but with dots aligned).
  • stackratio: how close to stack the dots. Default is 1, where dots just just touch. Use smaller values for closer, overlapping dots.
  • color, fill: Dot border color and area fill
  • dotsize: The diameter of the dots relative to binwidth, default 1.

As for violin plots, summary statistics are usually added to dot plots.

1. Create basic dot plots:

# Violin plots with mean points +/- SD
e + geom_dotplot(
  binaxis = "y", stackdir = "center",
  fill = "lightgray"
  ) + 
  stat_summary(
    fun.data = "mean_sdl", fun.args = list(mult=1), 
    geom = "pointrange", color = "red"
    )

# Combine with box plots
e + geom_boxplot(width = 0.5) + 
  geom_dotplot(
    binaxis = "y", stackdir = "center",
    fill = "white"
    ) 

# Dot plot + violin plot + stat summary
e + geom_violin(trim = FALSE) +
  geom_dotplot(
    binaxis='y', stackdir='center',
    color = "black", fill = "#999999"
    ) +
  stat_summary(
    fun.data="mean_sdl",  fun.args = list(mult=1), 
    geom = "pointrange", color = "#FC4E07", size = 0.4
    )

2. Create dot plots with multiple groups:

# Color dots by groups
e + geom_boxplot(width = 0.5, size = 0.4) +
  geom_dotplot(
    aes(fill = supp), trim = FALSE,
    binaxis='y', stackdir='center'
  )+
  scale_fill_manual(values = c("#00AFBB", "#E7B800"))

# Change the position : interval between dot plot of the same group
e + geom_boxplot(
  aes(color = supp), width = 0.5, size = 0.4,
  position = position_dodge(0.8)
  ) +
  geom_dotplot(
    aes(fill = supp, color = supp), trim = FALSE,
    binaxis='y', stackdir='center', dotsize = 0.8,
    position = position_dodge(0.8)
  )+
  scale_fill_manual(values = c("#00AFBB", "#E7B800"))+
  scale_color_manual(values = c("#00AFBB", "#E7B800"))

Stripcharts

Stripcharts are also known as one dimensional scatter plots. These plots are suitable compared to box plots when sample sizes are small.

• Key function: geom_jitter()
• key arguments: color, fill, size, shape. Changes points color, fill, size and shape

1. Create a basic stripchart:

• Change points shape and color by groups
• Adjust the degree of jittering: position_jitter(0.2)
• Add summary statistics:

e + geom_jitter(
  aes(shape = dose, color = dose), 
  position = position_jitter(0.2),
  size = 1.2
  ) +
  stat_summary(
    aes(color = dose),
    fun.data="mean_sdl",  fun.args = list(mult=1), 
    geom = "pointrange",  size = 0.4
    )+
  scale_color_manual(values =  c("#00AFBB", "#E7B800", "#FC4E07"))

2. Create stripcharts for multiple groups. The R code is similar to what we have seen in dot plots section. However, to create dodged jitter points, you should use the function position_jitterdodge() instead of position_dodge().

e + geom_jitter(
  aes(shape = supp, color = supp), 
  position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.8),
  size = 1.2
  ) +
  stat_summary(
    aes(color = supp),
    fun.data="mean_sdl",  fun.args = list(mult=1), 
    geom = "pointrange",  size = 0.4,
    position = position_dodge(0.8)
    )+
  scale_color_manual(values =  c("#00AFBB", "#E7B800"))

Sinaplot

sinaplot is inspired by the strip chart and the violin plot. By letting the normalized density of points restrict the jitter along the x-axis, the plot displays the same contour as a violin plot, but resemble a simple strip chart for small number of data points (Sidiropoulos et al. 2015).

In this way the plot conveys information of both the number of data points, the density distribution, outliers and spread in a very simple, comprehensible and condensed format.

Key function: geom_sina() [ggforce]:

library(ggforce)
# Create some data
d1 <- data.frame(
  y = c(rnorm(200, 4, 1), rnorm(200, 5, 2), rnorm(400, 6, 1.5)),
  group = rep(c("Grp1", "Grp2", "Grp3"), c(200, 200, 400))
  )

# Sinaplot
ggplot(d1, aes(group, y)) +
  geom_sina(aes(color = group), size = 0.7)+
  scale_color_manual(values =  c("#00AFBB", "#E7B800", "#FC4E07"))

Mean and median plots with error bars

In this section, we’ll show how to plot summary statistics of a continuous variable organized into groups by one or multiple grouping variables.

Note that, an easy way, with less typing, to create mean/median plots, is provided in the ggpubr package. See the associated article at: ggpubr-Plot Means/Medians and Error Bars

Set the default theme to theme_pubr() [in ggpubr]:

theme_set(ggpubr::theme_pubr())

1. Basic mean/median plots. Case of one continuous variable and one grouping variable:

• Prepare the data: ToothGrowth data set.

df <- ToothGrowth
df$dose <- as.factor(df$dose)
head(df, 3)

##    len supp dose
## 1  4.2   VC  0.5
## 2 11.5   VC  0.5
## 3  7.3   VC  0.5

• Compute summary statistics for the variable len organized into groups by the variable dose:

library(dplyr)
df.summary <- df %>%
  group_by(dose) %>%
  summarise(
    sd = sd(len, na.rm = TRUE),
    len = mean(len)
  )
df.summary

## # A tibble: 3 x 3
##     dose    sd   len
##     
## 1    0.5  4.50  10.6
## 2      1  4.42  19.7
## 3      2  3.77  26.1

• Create error plots using the summary statistics data. Key functions:
  • geom_crossbar() for hollow bar with middle indicated by horizontal line
  • geom_errorbar() for error bars
  • geom_errorbarh() for horizontal error bars
  • geom_linerange() for drawing an interval represented by a vertical line
  • geom_pointrange() for creating an interval represented by a vertical line, with a point in the middle.

Start by initializing ggplot with the summary statistics data: - Specify x and y as usually - Specify ymin = len-sd and ymax = len+sd to add lower and upper error bars. If you want only to add upper error bars but not the lower ones, use ymin = len (instead of len-sd) and ymax = len+sd.

# Initialize ggplot with data
f <- ggplot(
  df.summary, 
  aes(x = dose, y = len, ymin = len-sd, ymax = len+sd)
  )

Possible error plots:

Create simple error plots:

# Vertical line with point in the middle
f + geom_pointrange()

# Standard error bars
f + geom_errorbar(width = 0.2) +
  geom_point(size = 1.5)

Create horizontal error bars. Put dose on y axis and len on x-axis. Specify xmin and xmax.

# Horizontal error bars with mean points
# Change the color by groups
ggplot(
  df.summary, 
  aes(x = len, y = dose, xmin = len-sd, xmax = len+sd)
  ) +
  geom_point(aes(color = dose)) +
  geom_errorbarh(aes(color = dose), height=.2)+
  theme_light()

• Add jitter points (representing individual points), dot plots and violin plots. For this, you should initialize ggplot with original data (df) and specify the df.summary data in the error plot function, here geom_pointrange().

# Combine with jitter points
ggplot(df, aes(dose, len)) +
  geom_jitter(
    position = position_jitter(0.2), color = "darkgray"
    ) + 
  geom_pointrange(
    aes(ymin = len-sd, ymax = len+sd),
    data = df.summary
    )

# Combine with violin plots
ggplot(df, aes(dose, len)) +
  geom_violin(color = "darkgray", trim = FALSE) + 
  geom_pointrange(
    aes(ymin = len-sd, ymax = len+sd),
    data = df.summary
    )

• Create basic bar/line plots of mean +/- error. So we need only the df.summary data.
  • 1. Add lower and upper error bars for the line plot: ymin = len-sd and ymax = len+sd.
  • 2. Add only upper error bars for the bar plot: ymin = len (instead of len-sd) and ymax = len+sd.

# (1) Line plot
ggplot(df.summary, aes(dose, len)) +
  geom_line(aes(group = 1)) +
  geom_errorbar( aes(ymin = len-sd, ymax = len+sd),width = 0.2) +
  geom_point(size = 2)

# (2) Bar plot
ggplot(df.summary, aes(dose, len)) +
  geom_bar(stat = "identity", fill = "lightgray", 
           color = "black") +
  geom_errorbar(aes(ymin = len, ymax = len+sd), width = 0.2) 

For line plot, you might want to treat x-axis as numeric:

df.sum2 <- df.summary
df.sum2$dose <- as.numeric(df.sum2$dose)
ggplot(df.sum2, aes(dose, len)) +
  geom_line() +
  geom_errorbar( aes(ymin = len-sd, ymax = len+sd),width = 0.2) +
  geom_point(size = 2)

• Bar and line plots + jitter points. We need the original df data for the jitter points and the df.summary data for the other geom layers.
  • 1. For the line plot: First, add jitter points, then add lines + error bars + mean points on top of the jitter points.
  • 2. For the bar plot: First, add the bar plot, then add jitter points + error bars on top of the bars.

# (1) Create a line plot of means + 
# individual jitter points + error bars 
ggplot(df, aes(dose, len)) +
  geom_jitter( position = position_jitter(0.2),
               color = "darkgray") + 
  geom_line(aes(group = 1), data = df.summary) +
  geom_errorbar(
    aes(ymin = len-sd, ymax = len+sd),
    data = df.summary, width = 0.2) +
  geom_point(data = df.summary, size = 2)

# (2) Bar plots of means + individual jitter points + errors
ggplot(df, aes(dose, len)) +
  geom_bar(stat = "identity", data = df.summary,
           fill = NA, color = "black") +
  geom_jitter( position = position_jitter(0.2),
               color = "black") + 
  geom_errorbar(
    aes(ymin = len-sd, ymax = len+sd),
    data = df.summary, width = 0.2) 

2. Mean/median plots for multiple groups. Case of one continuous variable (len) and two grouping variables (dose, supp).

• Compute the summary statistics of len grouped by dose and supp:

library(dplyr)
df.summary2 <- df %>%
  group_by(dose, supp) %>%
  summarise(
    sd = sd(len),
    len = mean(len)
  )
df.summary2

## # A tibble: 6 x 4
## # Groups:   dose [?]
##     dose   supp    sd   len
##      
## 1    0.5     OJ  4.46 13.23
## 2    0.5     VC  2.75  7.98
## 3      1     OJ  3.91 22.70
## 4      1     VC  2.52 16.77
## 5      2     OJ  2.66 26.06
## 6      2     VC  4.80 26.14

• Create error plots for multiple groups:
  • 1. pointrange colored by groups (supp)
  • 2. standard error bars + mean points colored by groups (supp)

# (1) Pointrange: Vertical line with point in the middle
ggplot(df.summary2, aes(dose, len)) +
  geom_pointrange(
    aes(ymin = len-sd, ymax = len+sd, color = supp),
    position = position_dodge(0.3)
    )+
  scale_color_manual(values = c("#00AFBB", "#E7B800"))


# (2) Standard error bars
ggplot(df.summary2, aes(dose, len)) +
  geom_errorbar(
    aes(ymin = len-sd, ymax = len+sd, color = supp),
    position = position_dodge(0.3), width = 0.2
    )+
  geom_point(aes(color = supp), position = position_dodge(0.3)) +
  scale_color_manual(values = c("#00AFBB", "#E7B800")) 

• Create simple line/bar plots for multiple groups.
  • 1. Line plots: change linetype by groups (supp)
  • 2. Bar plots: change fill color by groups (supp)

# (1) Line plot + error bars
ggplot(df.summary2, aes(dose, len)) +
  geom_line(aes(linetype = supp, group = supp))+
  geom_point()+
  geom_errorbar(
    aes(ymin = len-sd, ymax = len+sd, group = supp),
     width = 0.2
    )

# (2) Bar plots + upper error bars.
ggplot(df.summary2, aes(dose, len)) +
  geom_bar(aes(fill = supp), stat = "identity",
           position = position_dodge(0.8), width = 0.7)+
  geom_errorbar(
    aes(ymin = len, ymax = len+sd, group = supp),
    width = 0.2, position = position_dodge(0.8)
    )+
  scale_fill_manual(values = c("grey80", "grey30"))

• Create easily plots of mean +/- sd for multiple groups. Use the ggpubr package, which will automatically calculate the summary statistics and create the graphs.

library(ggpubr)
# Create line plots of means
ggline(ToothGrowth, x = "dose", y = "len", 
       add = c("mean_sd", "jitter"),
       color = "supp", palette = c("#00AFBB", "#E7B800"))

# Create bar plots of means
ggbarplot(ToothGrowth, x = "dose", y = "len", 
          add = c("mean_se", "jitter"),
          color = "supp", palette = c("#00AFBB", "#E7B800"),
          position = position_dodge(0.8))

• Use the standard ggplot2 verbs, to reproduce the line plots above:

# Create line plots
ggplot(df, aes(dose, len)) +
  geom_jitter(
    aes(color = supp),
    position = position_jitter(0.2)
    ) + 
  geom_line(
    aes(group = supp, color = supp),
    data = df.summary2
    ) +
  geom_errorbar(
    aes(ymin = len-sd, ymax = len+sd, color = supp),
    data = df.summary2, width = 0.2
    )+
  scale_color_manual(values = c("#00AFBB", "#E7B800"))

Add p-values and significance levels

In this section, we’ll describe how to easily i) compare means of two or multiple groups; ii) and to automatically add p-values and significance levels to a ggplot (such as box plots, dot plots, bar plots and line plots, …).

Key functions:

• compare_means() [ggpubr package]: easy to use solution to performs one and multiple mean comparisons.
• stat_compare_means() [ggpubr package]: easy to use solution to automatically add p-values and significance levels to a ggplot.

The most common methods for comparing means include:

1. Compare two independent groups:

• Compute t-test:

library(ggpubr)
compare_means(len ~ supp, data = ToothGrowth,
              method = "t.test")

## # A tibble: 1 x 8
##     .y. group1 group2      p  p.adj p.format p.signif method
##                     
## 1   len     OJ     VC 0.0606 0.0606    0.061       ns T-test

• Create a box plot with p-values. Use the option method = "t.test" or method = "wilcox.test". Default is wilcoxon test.

# Create a simple box plot and add p-values
p <- ggplot(ToothGrowth, aes(supp, len)) +
  geom_boxplot(aes(color = supp)) +
  scale_color_manual(values = c("#00AFBB", "#E7B800"))
p + stat_compare_means(method = "t.test")

# Display the significance level instead of the p-value
# Adjust label position
p + stat_compare_means(
  aes(label = ..p.signif..), label.x = 1.5, label.y = 40
  )

2. Compare two paired samples. Use ggpaired() [ggpubr] to create the paired box plot.

ggpaired(ToothGrowth, x = "supp", y = "len",
         color = "supp", line.color = "gray", line.size = 0.4,
         palette = "jco")+
  stat_compare_means(paired = TRUE)

3. Compare more than two groups. If the grouping variable contains more than two levels, then pairwise tests will be performed automatically. The default method is “wilcox.test”. You can change this to “t.test”.

# Perorm pairwise comparisons
compare_means(len ~ dose,  data = ToothGrowth)

## # A tibble: 3 x 8
##     .y. group1 group2        p    p.adj p.format p.signif   method
##                           
## 1   len    0.5      1 7.02e-06 1.40e-05  7.0e-06     **** Wilcoxon
## 2   len    0.5      2 8.41e-08 2.52e-07  8.4e-08     **** Wilcoxon
## 3   len      1      2 1.77e-04 1.77e-04  0.00018      *** Wilcoxon

# Visualize: Specify the comparisons you want
my_comparisons <- list( c("0.5", "1"), c("1", "2"), c("0.5", "2") )
ggboxplot(ToothGrowth, x = "dose", y = "len",
          color = "dose", palette = "jco")+ 
  stat_compare_means(comparisons = my_comparisons)+ 
  stat_compare_means(label.y = 50)                   

4. Multiple grouping variables:

• (1/2). Create a multi-panel box plots facetted by group (here, “dose”):

# Use only p.format as label. Remove method name.
ggplot(ToothGrowth, aes(supp, len)) +
  geom_boxplot(aes(color = supp))+
  facet_wrap(~dose) +
  scale_color_manual(values = c("#00AFBB", "#E7B800")) +
  stat_compare_means(label = "p.format")

• (2/2). Create one single panel with all box plots. Plot y = “len” by x = “dose” and color by “supp”. Specify the option group in stat_compare_means():

ggplot(ToothGrowth, aes(dose, len)) +
  geom_boxplot(aes(color = supp))+
  scale_color_manual(values = c("#00AFBB", "#E7B800")) +
  stat_compare_means(aes(group = supp), label = "p.signif")

• Paired comparisons for multiple groups:

# Box plot facetted by "dose"
p <- ggpaired(ToothGrowth, x = "supp", y = "len",
          color = "supp", palette = "jco", 
          line.color = "gray", line.size = 0.4,
          facet.by = "dose", short.panel.labs = FALSE)
# Use only p.format as label. Remove method name.
p + stat_compare_means(label = "p.format", paired = TRUE)

Read more at: Add P-values and Significance Levels to ggplots