lab_graphics.Rmd

---
title: "Graphics in R"
subtitle: "R Programming Foundation for Life Scientists"
output:
  bookdown::html_document2:
    highlight: textmate
    toc: true
    toc_float:
      collapsed: true
      smooth_scroll: true
      print: false
    toc_depth: 4
    number_sections: true
    df_print: default
    code_folding: none
    self_contained: false
    keep_md: false
    encoding: 'UTF-8'
    css: "assets/lab.css"
    include:
      after_body: assets/footer-lab.html
---

```{r,child="assets/header-lab.Rmd"}
```

```{r,include=FALSE}
library(lubridate)
```

# Introduction

In this lab, we will go step-by-step through manually building a scientific plot using base graphics in R.

# Generating data

First, we will produce some random data that we will later plot. Make a data frame with

- 20 random coordinates (x,y) and
- radius $r$ for each data point.
- The *x* coord takes random values from 1 to 25 and
- both the *y* and the radius *r* coord are samples from $N(0,1)$.
- Each point (row of the data frame) has a name *ind1 ... ind25*.

First, look at the defaults:

- plot the data in the simplest possible way.

```{r,fig.height=6,fig.width=10}
#20 random datapoints
x <- sample(c(1:25), size=20, replace=T)
y <- rnorm(n=20, mean=0, sd=1) # sample from normal
r <- rnorm(n=20, mean=0, sd=1) # radius from normal
names <- paste("ind", 1:20, sep="") # assign some names
data <- data.frame(cbind(X=x,Y=y, R=r), row.names=names)
plot(data[,1:2])
```

# Building a plot

As you see, the points are displayed in a simple way, axes are set automatically, the radius is not reflected on the plot in any way (3rd dimension).

## Hide markers

Build the plot from scratch, begin by displaying no points. You can do this by setting `type = 'n'`.   

* Other `type`s:   
  + "p" for points.   
  + "l" for lines.   
  + "b" for both points and lines.  
  + "c" for empty points joined by lines.  
  + "o" for overplotted points and lines.  
  + "s" and "S" for stair steps.   
  + "h" for histogram-like vertical lines.   
  + "n" does not produce any points or lines.   
  

```{r,fig.height=6,fig.width=10}
plot(data[,1:2], type='n')
```

## Hide border/axes

Remove the default box around the plot and axes.

```{r,fig.height=6,fig.width=10}
plot(data[,1:2], type='n',xaxt='n', yaxt='n',xlab="", ylab="", frame.plot=F)
```

## Axes labels

Create X and Y axis so that they cover the whole range of *x* and *y*. For the Y axis, set 10 equidistant tickmarks and set labels to their values rounded to two decimals. Turn the labels, so that they are parallel to the OX axis.

```{r,fig.height=6,fig.width=10,echo=FALSE}
## CODE EVALUATED BUT HIDDEN FROM OUTPUT ##
plot(data[,1:2], type='n',xaxt='n', yaxt='n',xlab="", ylab="", frame.plot=F)

#Create X axis
coords.x <- seq(min(data$X),max(data$X), by=1)
axis(side=1, # 1-left, 2-top, 3-right, 4-bottom
  at=coords.x # coordinates for tickmarks
)

#Create Y axis
#we want 10 tickmarks along the data range
coords.y <- seq(min(data$Y), max(data$Y), length.out=10)

#and our labels will be the rounded values of y
labels.y <- round(coords.y, digits=2)
axis(side=2,
  at=coords.y,
  labels=labels.y, # we want specific labels
  las=2 # turn the text so it is parallel to OX
)
```

```{r,eval=FALSE}
#Create X axis
coords.x <- seq(min(data$X),max(data$X), by=1)
axis(side=1, # 1-left, 2-top, 3-right, 4-bottom
  at=coords.x # coordinates for tickmarks
)

#Create Y axis
#we want 10 tickmarks along the data range
coords.y <- seq(min(data$Y), max(data$Y), length.out=10)

#and our labels will be the rounded values of y
labels.y <- round(coords.y, digits=2)
axis(side=2,
  at=coords.y,
  labels=labels.y, # we want specific labels
  las=2 # turn the text so it is parallel to OX
)
```

## Grid lines

Plot gridlines so that it is easier to read the plot. There should be a grey dashed line from each tickmark on both axes.

```{r,fig.height=6,fig.width=10,echo=FALSE}
## CODE EVALUATED BUT HIDDEN FROM OUTPUT ##

plot(data[,1:2], type='n',xaxt='n', yaxt='n',xlab="", ylab="", frame.plot=F)

#Create X axis
coords.x <- seq(min(data$X),max(data$X), by=1)
axis(side=1, # 1-left, 2-top, 3-right, 4-bottom
  at=coords.x # coordinates for tickmarks
)

#Create Y axis
#we want 10 tickmarks along the data range
coords.y <- seq(min(data$Y), max(data$Y), length.out=10)

#and our labels will be the rounded values of y
labels.y <- round(coords.y, digits=2)
axis(side=2,
  at=coords.y,
  labels=labels.y, # we want specific labels
  las=2 # turn the text so it is parallel to OX
)

abline(v=coords.x, col="darkgrey", lty=3)
abline(h=coords.y, col="darkgrey", lty=3)
```

```{r,fig.height=6,fig.width=10,eval=FALSE}
abline(v=coords.x, col="darkgrey", lty=3)
abline(h=coords.y, col="darkgrey", lty=3)
#you could also use grid()
```

## Add markers

Define a new *mycol* function that takes a color name and a transparency value as two arguments and returns the corresponding rgb color value. 
  
```{r,fig.height=6,fig.width=10}
#Function for adding transparency to a given color.
mycol <- function(colname="olivedrab", transparency=.5) {
#convert color name to its RGB value and add the desired
#transparency
color <- c(as.vector(col2rgb(colname))/255, transparency)
# and make a new color from the above
color <- rgb(color[1], color[2], color[3], color[4])
return(color)
}
```

Plot datapoints so that their size is proportional to $e^r$ where $r$ is the radius, points at even X should be round and blue and points at odd X square and grey.

```{r,fig.height=6,fig.width=10,echo=FALSE}
## CODE EVALUATED BUT HIDDEN FROM OUTPUT ##

plot(data[,1:2], type='n',xaxt='n', yaxt='n',xlab="", ylab="", frame.plot=F)

#Create X axis
coords.x <- seq(min(data$X),max(data$X), by=1)
axis(side=1, # 1-left, 2-top, 3-right, 4-bottom
  at=coords.x # coordinates for tickmarks
)

#Create Y axis
#we want 10 tickmarks along the data range
coords.y <- seq(min(data$Y), max(data$Y), length.out=10)

#and our labels will be the rounded values of y
labels.y <- round(coords.y, digits=2)
axis(side=2,
  at=coords.y,
  labels=labels.y, # we want specific labels
  las=2 # turn the text so it is parallel to OX
)

abline(v=coords.x, col="darkgrey", lty=3)
abline(h=coords.y, col="darkgrey", lty=3)

points(data[data$X%%2 == 0,], pch=19, cex=exp(r), col=mycol("slateblue", .5))
points(data[data$X%%2 != 0,], pch=15, cex=exp(r), col=mycol("grey", .5))
```

```{r,fig.height=6,fig.width=10,eval=FALSE}
#Plot radii
points(data[data$X%%2 == 0,], pch=19, cex=exp(r), col=mycol("slateblue", .5))
points(data[data$X%%2 != 0,], pch=15, cex=exp(r), col=mycol("grey", .5))
```

Plot centers of the points as a cross: grey for blue/even points and red for grey/odd points.

```{r,fig.height=6,fig.width=10,echo=FALSE}
## CODE EVALUATED BUT HIDDEN FROM OUTPUT ##

plot(data[,1:2], type='n',xaxt='n', yaxt='n',xlab="", ylab="", frame.plot=F)

#Create X axis
coords.x <- seq(min(data$X),max(data$X), by=1)
axis(side=1, # 1-left, 2-top, 3-right, 4-bottom
  at=coords.x # coordinates for tickmarks
)

#Create Y axis
#we want 10 tickmarks along the data range
coords.y <- seq(min(data$Y), max(data$Y), length.out=10)

#and our labels will be the rounded values of y
labels.y <- round(coords.y, digits=2)
axis(side=2,
  at=coords.y,
  labels=labels.y, # we want specific labels
  las=2 # turn the text so it is parallel to OX
)

abline(v=coords.x, col="darkgrey", lty=3)
abline(h=coords.y, col="darkgrey", lty=3)

points(data[data$X%%2 == 0,], pch=19, cex=exp(r), col=mycol("slateblue", .5))
points(data[data$X%%2 != 0,], pch=15, cex=exp(r), col=mycol("grey", .5))

points(data[data$X%%2 == 0,], pch=3, cex=1, col="darkgrey")
points(data[data$X%%2 != 0,], pch=3, cex=1, col="red")
```

```{r,fig.height=6,fig.width=10,eval=FALSE}
points(data[data$X%%2 == 0,], pch=3, cex=1, col="darkgrey")
points(data[data$X%%2 != 0,], pch=3, cex=1, col="red")
```

## Annotation

Add grey text 'Center' at the center of the plot.

```{r,fig.height=6,fig.width=10,echo=FALSE}
## CODE EVALUATED BUT HIDDEN FROM OUTPUT ##

plot(data[,1:2], type='n',xaxt='n', yaxt='n',xlab="", ylab="", frame.plot=F)

#Create X axis
coords.x <- seq(min(data$X),max(data$X), by=1)
axis(side=1, # 1-left, 2-top, 3-right, 4-bottom
  at=coords.x # coordinates for tickmarks
)

#Create Y axis
#we want 10 tickmarks along the data range
coords.y <- seq(min(data$Y), max(data$Y), length.out=10)

#and our labels will be the rounded values of y
labels.y <- round(coords.y, digits=2)
axis(side=2,
  at=coords.y,
  labels=labels.y, # we want specific labels
  las=2 # turn the text so it is parallel to OX
)

abline(v=coords.x, col="darkgrey", lty=3)
abline(h=coords.y, col="darkgrey", lty=3)

points(data[data$X%%2 == 0,], pch=19, cex=exp(r), col=mycol("slateblue", .5))
points(data[data$X%%2 != 0,], pch=15, cex=exp(r), col=mycol("grey", .5))

points(data[data$X%%2 == 0,], pch=3, cex=1, col="darkgrey")
points(data[data$X%%2 != 0,], pch=3, cex=1, col="red")

center.x <- mean(range(data[,1]))
center.y <- mean(range(data[,2]))
text(x=center.x, y=center.y, "Center", col="lightgrey")
```

```{r,fig.height=6,fig.width=10,eval=FALSE}
center.x <- mean(range(data[,1]))
center.y <- mean(range(data[,2]))
text(x=center.x, y=center.y, "Center", col="lightgrey")
```

## Titles

Add title 'Odds and Ends' and text 'X' and 'Y' on the margins of the appropriate axes.

```{r,fig.height=6,fig.width=10,echo=FALSE}
## CODE EVALUATED BUT HIDDEN FROM OUTPUT ##

plot(data[,1:2], type='n',xaxt='n', yaxt='n',xlab="", ylab="", frame.plot=F)

#Create X axis
coords.x <- seq(min(data$X),max(data$X), by=1)
axis(side=1, # 1-left, 2-top, 3-right, 4-bottom
  at=coords.x # coordinates for tickmarks
)

#Create Y axis
#we want 10 tickmarks along the data range
coords.y <- seq(min(data$Y), max(data$Y), length.out=10)

#and our labels will be the rounded values of y
labels.y <- round(coords.y, digits=2)
axis(side=2,
  at=coords.y,
  labels=labels.y, # we want specific labels
  las=2 # turn the text so it is parallel to OX
)

abline(v=coords.x, col="darkgrey", lty=3)
abline(h=coords.y, col="darkgrey", lty=3)

points(data[data$X%%2 == 0,], pch=19, cex=exp(r), col=mycol("slateblue", .5))
points(data[data$X%%2 != 0,], pch=15, cex=exp(r), col=mycol("grey", .5))

points(data[data$X%%2 == 0,], pch=3, cex=1, col="darkgrey")
points(data[data$X%%2 != 0,], pch=3, cex=1, col="red")

center.x <- mean(range(data[,1]))
center.y <- mean(range(data[,2]))
text(x=center.x, y=center.y, "Center", col="lightgrey")

title("Odds and Ends")
mtext("Y", side=2, line=3, cex.lab=1,las=2, col="blue")
mtext("X", side=1, line=3, cex.lab=1,las=1, col="blue")
```

```{r,fig.height=6,fig.width=10,eval=FALSE}
title("Odds and Ends")
mtext("Y", side=2, line=3, cex.lab=1,las=2, col="blue")
mtext("X", side=1, line=3, cex.lab=1,las=1, col="blue")
```

## Legend

Add a legend for 'odd' and 'even' points. Place it in the top-right corner.

```{r,fig.height=6,fig.width=10,echo=FALSE}
## CODE EVALUATED BUT HIDDEN FROM OUTPUT ##

plot(data[,1:2], type='n',xaxt='n', yaxt='n',xlab="", ylab="", frame.plot=F)

#Create X axis
coords.x <- seq(min(data$X),max(data$X), by=1)
axis(side=1, # 1-left, 2-top, 3-right, 4-bottom
  at=coords.x # coordinates for tickmarks
)

#Create Y axis
#we want 10 tickmarks along the data range
coords.y <- seq(min(data$Y), max(data$Y), length.out=10)

#and our labels will be the rounded values of y
labels.y <- round(coords.y, digits=2)
axis(side=2,
  at=coords.y,
  labels=labels.y, # we want specific labels
  las=2 # turn the text so it is parallel to OX
)

abline(v=coords.x, col="darkgrey", lty=3)
abline(h=coords.y, col="darkgrey", lty=3)

points(data[data$X%%2 == 0,], pch=19, cex=exp(r), col=mycol("slateblue", .5))
points(data[data$X%%2 != 0,], pch=15, cex=exp(r), col=mycol("grey", .5))

points(data[data$X%%2 == 0,], pch=3, cex=1, col="darkgrey")
points(data[data$X%%2 != 0,], pch=3, cex=1, col="red")

center.x <- mean(range(data[,1]))
center.y <- mean(range(data[,2]))
text(x=center.x, y=center.y, "Center", col="lightgrey")

title("Odds and Ends")
mtext("Y", side=2, line=3, cex.lab=1,las=2, col="blue")
mtext("X", side=1, line=3, cex.lab=1,las=1, col="blue")

legend('topright',
	legend=c("odd", "even"),
	col=c(mycol("slateblue", .5), mycol("grey", .5)),
	pch=c(19,15),
	cex=1,
	pt.cex=1.2,
	title="Legend",
	bty='n'
)
```

```{r,fig.height=6,fig.width=10,eval=FALSE}
legend('topright',
	legend=c("odd", "even"),
	col=c(mycol("slateblue", .5), mycol("grey", .5)),
	pch=c(19,15),
	cex=1,
	pt.cex=1.2,
	title="Legend",
	bty='n'
)
```

# Visualizing growth data

A female child was measured at the following dates:

- '30-09-2015', '12-10-2015', '19-10-2015', '26-10-2015', '07-11-2015', '16-11-2015', '30-11-2015', '11-01-2016', '08-02-2016', '14-03-2016', '05-04-2016', '14-04-2016', '31-05-2016', '14-07-2016',

- the measured weights in grams were: 3300, 3540, 3895, 4070, 4230, 4385, 4855, 5865, not taken, 6736, 7065, 7080, 7530, 7640 and
- the measured lengths: 43, no measurement taken, 53, 54, 55, 56, 58, 62.5, 65, 67, 67.5, 67.5, 70.5, 71.5.
- The headcircumference for the same datapoints was (in cm): 34, 35.5, 36.1, 36.8, 36.8, 37.3, 38, 40.2, 41.4, 42.1, not taken, 43, 44, 45.

Your task is to plot these data on the WHO centile grids. Choose weight/length/circumference depending on the month you was born:

  - weight: Jan, Apr, Jul, Oct
  - length: Feb, May, Aug, Nov
  - circumference: Mar, Jun, Sep, Dec

## Prepare input data

We save the timepoints and measurements in vectors.  
```{r}
library(lubridate)
timepoints <- dmy(c('30-09-2015', '12-10-2015','19-10-2015', '26-10-2015', '07-11-2015', '16-11-2015','30-11-2015', '11-01-2016', '08-02-2016', '14-03-2016', '05-04-2016', '14-04-2016', '31-05-2016', '14-07-2016'))

weight <- c(3300, 3540, 3895, 4070, 4230, 4385, 4855, 5865, NA, 6736, 7065, 7080, 7530, 7640)
length <- c(43,NA,53,54,55,56,58,62.5,65,67,67.5,67.5,70.5,71.5)
head <- c(34,35.5,36.1,36.8,36.8,37.3,38,40.2,41.4,42.1,NA,43,44,45)
```

We can calculate the interval between the dates of measurements by two approaches:  
**Approach1:**  
Simply calculate the interval by subtracting each date from the the first date of data collection (e.g. '2016-05-31' - '2015-09-30' )  and then convert it to months.  

```{r}
xpoints <- (as.Date(timepoints) - as.Date(timepoints[1]) )/ 30
```

**Approach2:**  
Calculate the intervals by seconds and then use WHO standard months day length which is 30.4375 to calculate by month.   
Use function `dmy()` from the **lubridate** package to create a vector of timepoints.

HINTS:   
- We can define an Interval using the `%--%` operator.  
- check `as.duration()` and `ddays()` functions.  



```{r, echo = T, eval = F}
who.month <- 30.4375 #days
xpoints <- as.duration(timepoints[1] %--% timepoints) / ddays(1) / who.month
```

## Prepare reference data

Download manually and then load:   
Go to WHO website (http://www.who.int/childgrowth/standards/en/) and find out the link to the dataset of your concern, e.g. Weight for age, percentiles for girls have the following address: http://www.who.int/entity/childgrowth/standards/wfa_girls_p_0_5.xlsx

Or using direct link to the excel file:  

```{r,eval=T}
library(readxl)
uri <- "https://cdn.who.int/media/docs/default-source/child-growth/child-growth-standards/indicators/weight-for-age/tab_wfa_girls_p_0_5.xlsx?sfvrsn=666fe445_7"
local_file_path <- "wfa_girls_p_0_5.xlsx" ## give the local path of the downloaded file
download.file(url = uri, destfile = local_file_path, mode = "wb")
myData <-read_excel(local_file_path)
```

## Build empty plot

Create an empty plot to show your and WHO data,

```{r,accordion=TRUE,fig.height=6,fig.width=10,eval=T}
plot(1, xlim=c(0, max(myData$Month)), type='n', bty='n', ylim=c(0, max(myData[,c(5:19)])), las=1, xlab='Month', ylab='kg')
grid()
```

## Plot reference data

Plot WHO mean and percentiles: P25, P75, P0.1 and P99.9, use different colors and line types to make the plot pretty.

```{r,fig.height=6,fig.width=10,echo=FALSE,eval=FALSE}
## CODE EVALUATED BUT NOT DISPLAYED ##

plot(1, xlim=c(0, max(myData$Month)), type='n', bty='n', ylim=c(0, max(myData[,c(5:19)])), las=1, xlab='Month', ylab='kg')
grid()

lines(myData$M, col='grey', lty=1)
lines(myData$P25, col='blue', lty=2)
lines(myData$P75, col='blue', lty=2)
lines(myData$P01, col='tomato', lty=2)
lines(myData$P999, col='tomato', lty=2)
```

```{r,accordion=TRUE,eval=FALSE}
lines(myData$M, col='grey', lty=1)
lines(myData$P25, col='blue', lty=2)
lines(myData$P75, col='blue', lty=2)
lines(myData$P01, col='tomato', lty=2)
lines(myData$P999, col='tomato', lty=2)
```

## Plot input data

Plot your data on top of the percentiles, mind the units so that they match with the WHO ones

```{r,fig.height=6,fig.width=10,echo=FALSE,eval=T}
## CODE EVALUATED BUT NOT DISPLAYED ##

plot(1, xlim=c(0, max(myData$Month)), type='n', bty='n', ylim=c(0, max(myData[,c(5:19)])), las=1, xlab='Month', ylab='kg')
grid()

lines(myData$M, col='grey', lty=1)
lines(myData$P25, col='blue', lty=2)
lines(myData$P75, col='blue', lty=2)
lines(myData$P01, col='tomato', lty=2)
lines(myData$P999, col='tomato', lty=2)

points(xpoints, weight/1000, pch=3, type='l', cex=.5)
points(xpoints, weight/1000, pch=3, type='p', cex=.5)
```

```{r,accordion=TRUE,eval=FALSE}
points(xpoints, weight/1000, pch=3, type='l', cex=.5)
points(xpoints, weight/1000, pch=3, type='p', cex=.5)
```

## Add annotation

Add descriptions of the confidence lines in the margins

```{r,fig.height=6,fig.width=10,echo=FALSE,eval=T}
## CODE EVALUATED BUT NOT DISPLAYED ##

plot(1, xlim=c(0, max(myData$Month)), type='n', bty='n', ylim=c(0, max(myData[,c(5:19)])), las=1, xlab='Month', ylab='kg')
grid()

lines(myData$M, col='grey', lty=1)
lines(myData$P25, col='blue', lty=2)
lines(myData$P75, col='blue', lty=2)
lines(myData$P01, col='tomato', lty=2)
lines(myData$P999, col='tomato', lty=2)

points(xpoints, weight/1000, pch=3, type='l', cex=.5)
points(xpoints, weight/1000, pch=3, type='p', cex=.5)

mtext(text = c('P0.1','P25','P75','P99.9'), side = 4, at=myData[dim(myData)[1], c('P01','P25','P75','P999')], las=1, cex=.8)
```

```{r,accordion=TRUE,eval=FALSE}
mtext(text = c('P0.1','P25','P75','P99.9'), side = 4, at=myData[dim(myData)[1], c('P01','P25','P75','P999')], las=1, cex=.8)
```

# Visualizing Gapminder data

You task here is to use the already acquired R knowledge to plot an interesting relationship between two freely selected variables available at Hans Rosling's Gapminder Foundation page.

- Go to http://www.gapminder.org/data/
- Select a dataset of interest,
- Load data to R, take care of missing values etc.,
- Find a nice way of visualizing the relationship between some selected variables,
- Think of scales (linear, logarythmic), axes labels etc.,
- Be creative,
- Visualize a selected variables using boxplot and histogram on one plot (HINT: parameter mfrow),
- Discuss the result with your colleagues and TAs.

```{r}
unlink(local_file_path)
```