06-scatterplots.Rmd

```{r, echo=FALSE}
knitr::opts_chunk$set(
  tidy = FALSE
  )
```
# Scatter-plots 

The basic plotting functions are very easy, but we first need some data to work with. For this example we are going to use one of the inbuilt datasets.

There are a number of datasets which come readily loaded with the `datasets` package. In this example we are using the `trees` dataset, which contains 'measurements of the girth, height and volume of timber in 31 felled black cherry trees' (See `?trees` for more information)

First, we will look at the data.

```{r}
trees
```

So far so good. To produce a basic scatter plot in R is very simple. We use the `plot` command. In this example, we need to specify which parts of the trees dataset we are going to use. Just as revision, we will try three different ways of doing this - all will produce the same plot, but each one has a subtly different use.

In the first example, we select a subscript from the dataframe using square brackets `[,1]`. Refer back to <a href="Tutorial/getting_data_into_R.html">Getting data into R</a> if you are unsure. This can be a useful way to select data from a larger dataset, if for instance we are dealing with matrices.

```{r Scatterplots-1}
plot(
     x = trees[,1],
     y = trees[,2]
     )
```

The second way, we will use the `$` operator to select a column from within the dataframe:

```{r Scatterplots-2}
plot(
     x = trees$Girth,
     y = trees$Height
     )
```

And finally, we can try the `with()` function. Using `with()` dispenses with the need for the `$` operator, and allows us to call objects within a dataframe without specifying the name of the dataframe every time.

```{r Scatterplots-3}

with(
     trees,
     plot(
          x = Girth,
          y = Height
          )
     )

```

It doesn't really matter which method you use in this example, but it is useful to revise different ways of getting at data within a dataframe

This plot is not particularly inspiring, so let's look at some of the ways we can make it nicer.

Within the `plot` command, there are a number of arguments that we can change in order to alter the appearance of the plot, these are [here](plot_arguments.html), lifted from the help file for `?plot`. You can get a similar list of arguments for any of the other plot functions (or indeed any function) by calling up the help file for that function (just place a "?" in front of the function: e.g. `?boxplot`).


### Point colour and shape 

So lets change some general things first. We'll label the axes (`xlab` or `ylab`)  and change the shape (`pch`) and colour (`col`) of the points.

```{r Scatterplots-4}

with(
     trees,
     plot(
          x = Girth,
          y = Height,
          xlab = "Girth (in)",
          ylab = "Height (ft)",
          pch = 16,
          col = "red"
          )
     )


```

A little bit better...

Note that there a wide range of symbol shapes available...

[![](http://rgraphics.limnology.wisc.edu/images/miscellaneous/pch.png)](http://rgraphics.limnology.wisc.edu/images/miscellaneous/pch.png)

...and pretty much any colour you can imagine. Colours can either be specified as a name (e.g. "red"), as a hexadecimal code (e.g. `col = "#FF0000"`) or in RGB format (e.g. `col = rgb(1, 0, 0)`). Try [http://www.colorpicker.com/](http://www.colorpicker.com) if you get stuck. Note that unless you are going to pay for colour figures in journal plates, or you can be sure that no colour blind people will read your work, it is good practice to use a combination of colours and symbols

We can use the package `RColorBrewer` to make picking complementary colour palettes a little easier. See `?RColorBrewer` or [http://colorbrewer2.org](http://colorbrewer2.org) for more information. There are a number of palettes for different applications. We will just use the first of them - this palette is good for colouring lines and points, but should be avoided for blocks of colour.


```{r, include=FALSE, echo=FALSE}
library(RColorBrewer)
```

We can see the various palettes with:

```{r}
display.brewer.pal(
     n = 9,
     name = "Set1"
     ) 
```

And create our palette with:

```{r}
pal <- brewer.pal(
     n = 9,
     name = "Set1"
     )
```

Note that `RColorBrewer` specifies colours with hexadecimals:

```{r}
pal
```

### Axes and outlines 

We are also able to customise the axes in any way we want, but it is usually necessary to remove the axes altogether with `xaxt="n"` and/or `yaxt="n"` and replace them with the `axis` command.

This example also demonstrates the use of the `bty` argument.

```{r Scatterplots-5}
with(
     trees,
     plot(
          x = Girth,
          y = Height,
          xlab = "Girth (in)",
          ylab = "Height (ft)",
          pch = 16,
          col = pal[2],         # second colour from the palette we produced.
          yaxt = "n",           # remove the y axis
          xaxt = "n",           # remove the y axis
          cex.lab = 1.2,      # scaling factor relating to the size of axis titles.
          cex = 1.4,          # scaling factor for the size of points
          bty = "l"             # this can be either "o", "l" or "u"
          )                   # and affects the box bordering the plot
     )

axis(
     side = 1,           # 1=below, 2=left, 3=above and 4=right.
     at = seq(0,20,4),   # a sequence denoting the location of axis tick marks
     cex.axis = 1.2      # a scaling factor of the default text size
     )
axis(
     side = 2, 
     at = seq(50,100,5), 
     las = 2,            # Rotate axis labels
     cex.axis = 1.2
     )


```

### Additional observations 

If we want to add additional observations to a plot, it is quite simple. We can use the `points` command, which is specified in a very similar way to the original `plot` command.

In this example, we will use the same `trees` dataset and plot as the previous example, but we will add a second set of observations which we have simulated using `jitter`. The `jitter` command adds a small amount of 'noise' to an observation by adding or removing a small amount from it.

```{r Scatterplots-6}

# Repeat the same code as the last example:

with(
     trees,
     plot(
          x = Girth,
          y = Height,
          xlab = "Girth (in)",
          ylab = "Height (ft)",
          pch = 16,
          col = pal[2],         # second colour from the palette we produced.
          yaxt = "n",           # remove the y axis
          xaxt = "n",           # remove the y axis
          cex.lab = 1.2,      # scaling factor relating to the size of axis titles.
          cex = 1.4,          # scaling factor for the size of points
          bty = "l"             # this can be either "o", "l" or "u"
          )                   # and affects the box bordering the plot
     )

axis(
     side = 1,           # 1=below, 2=left, 3=above and 4=right.
     at = seq(0,20,4),   # a sequence denoting the location of axis tick marks
     cex.axis = 1.2      # a scaling factor of the default text size
     )
axis(
     side = 2, 
     at = seq(50,100,5), 
     las = 2,            # Rotate axis labels
     cex.axis = 1.2
     )

# Add an additonal set of points:

with(
     trees,
     points(
          x = jitter(Girth, amount=2),       # See ?jitter for details
          y = jitter(Height, amount = 2),
          pch = 15,
          col = pal[1],
          cex = 1.4
          )
     )

# We also now need to add a legend:

legend(
     "bottomright",      # this can also take the form of x and y coordinates
     legend = c("Jittered","Original"), # labels for the legend
     col = pal[1:2],     # colour in order of appearance
     pch = c(15:16),     # pch in order of appearance
     pt.cex = 1.4,       # size of legend points
     bty="n"             # remove the black outline from the legend
     )

```

Note that in the above example we would need to adjust the limits of the x and y axes in order to accommodate all of the extra points which we added. This can be done with the arguments `xlim=c(1,2)` or `ylim=c(1,2)`, substituting the lowest and highest numbers for the 1 and 2.

### Log transformations 

It is sometimes necessary to transform your data when plotting. This is quite easy to achieve in R.

In this example, we will import directly from an online text file (also from Crawley's R book).

```{r}
sapdecay <- read.table(
     "http://www.bio.ic.ac.uk/research/mjcraw/therbook/data/sapdecay.txt", 
     header = T
     )
```

 We'll call the data to examine it. Note this data is dealt with on p70 of 'The R Book'. remember that typing the name of the object `sapdecay` is the 'short-hand' equivalent of typing `print(sapdecay)` - there are certain occasions when it will be necessary to use the full command - more of that later.

```{r}
sapdecay
```

<div class="ex">

### Exercise

Now plot the data producing the graph below. Remember to check the arguments listed above if you don't know how to do something.
</div>

```{r Scatterplots-7}
plot(
  x = sapdecay$x,
  y = sapdecay$y,
  bty = "l",
  ylab = "Radioactive decay",
  xlab = "Days",
  main = "Radioactive Decay Plot",
  col = pal[2],
  pch = 1,
  cex = 1.4,
  cex.axis = 1.2,
  cex.lab = 1.2
)

```

It is possible for us to fit a non-linear model to this data. We're not going to cover non-linear modelling in this material, but it is covered elsewhere in the R book.

The code below (.70 in the R book) fits an exponential model to this data

```{r, Scatterplots-8, echo=FALSE, echo=4:5}
plot(
  x = sapdecay$x,
  y = sapdecay$y,
  bty = "l",
  ylab = "Radioactive decay",
  xlab = "Days",
  main = "Radioactive Decay Plot",
  col = pal[2],
  pch = 1,
  cex = 1.4,
  cex.axis = 1.2,
  cex.lab = 1.2
)

xv<-seq(
  from = 0,
  to = 50,
  by = 0.1
)

lines(
  x = xv,
  y = exp(-0.055*xv)
)

```

In general however, it is easier if we can linearise our data somehow. For plots, this can be achieved in two ways in R. Either we add the command `log="x"` to convert the x-axis to a base 10 logarithmic scale (or you can specify `"y"` or `"xy"` for both), or we can convert the data itself before plotting using the `log10()` command.

<div class="ex">

### Exercise 

Reproduce the above graph again, but this time try specifying `log="y"`.Then try by converting `y` to a logarithmic scale with `log10(y)`. You should produce the two plots below:


</div>


```{r Scatterplots-9, echo=FALSE, fig.width=12}
par(mfrow=c(2,1))

plot(
  x = sapdecay$x,
  y = sapdecay$y,
  bty = "l",
  ylab = "Radioactive decay",
  xlab = "Days",
  main = "Radioactive Decay Plot",
  col = pal[2],
  pch = 1,
  cex = 1.4,
  cex.axis = 1.2,
  cex.lab = 1.2,
  log = "y"
)

plot(
  x = sapdecay$x,
  y = log10(sapdecay$y),
  bty = "l",
  ylab = "Radioactive decay",
  xlab = "Days",
  main = "Radioactive Decay Plot",
  col = pal[2],
  pch = 1,
  cex = 1.4,
  cex.axis = 1.2,
  cex.lab = 1.2
)

par(mfrow = c(1,1))
```

Notice that when we log the axes rather than the data, the axis labels remain a lot more comprehensible.

We are now able to fit a linear model to this data, which is mathematically simpler than specifying a non-linear model.

```{r Scatterplots-10, echo=4:5}

plot(
  x = sapdecay$x,
  y = sapdecay$y,
  bty = "l",
  ylab = "Radioactive decay",
  xlab = "Days",
  main = "Radioactive Decay Plot",
  col = pal[2],
  pch = 1,
  cex = 1.4,
  cex.axis = 1.2,
  cex.lab = 1.2,
  log = "y"
)


abline(
  lm(
    formula = log10(y)~x, 
    data = sapdecay
  )
)
```


### Customising white space 

So by now, you may be realising that pretty much everything in R plots can be customised. This includes the amount of whitespace you leave around your plots. This can very useful to change when you are creating two plots next to each other.

I'm not going to dwell on how to do this, because there are already some excellent guides written on the internet (see: <http://research.stowers-institute.org/efg/R/Graphics/Basics/mar-oma/>).

One of the plots from this page is included below, showing which commands control whitespace around a plot. Maybe take a deep breath first before looking at the code on the website!

```{r Scatterplots-11, echo=FALSE}

GenericFigure <- function(ID, size1, size2)
{
  plot(0:10, 0:10, type="n", xlab="X", ylab="Y")
  text(5,5, ID, col="red", cex=size1)
  box("plot", col="red")
  mtext("Figure",
        SOUTH<-1, line=3, adj=1.0, cex=size2,   col="blue")
}


# Figure 1.  Default R Plot Area / Figure Area
Figure1 <- function()
{
  oldpar <- par(oma=c(0,0,0,0))  # default values
  GenericFigure("Plot Area", 3,2)
  box("figure", col="blue")
  par(oldpar)
}

Figure2A <- function()
{
  par(oma=c(2,2,2,2))
  GenericFigure("Plot Area", 3,2)

  Margins <- capture.output( par()$mar )
  Margins <- substr(Margins, 5, nchar(Margins))
  Margins <-
       paste("mar = c(", gsub(" ",",",Margins), ")", sep="")
  mtext(Margins,
        NORTH<-3, line=2, adj=0.0, cex=1.5, col="blue")

  # "figure" box and "inner" margin box same for single figure plot
  box("figure",lty="dashed", col="blue")
  box("inner", lty="dotted", col="green")

  mtext("Outer Margin Area",
        SOUTH<-1, line=0.4, adj=1.0, cex=1.5, col="green", outer=TRUE)
  box("outer", lty="solid", col="green")

  OuterMargins <- capture.output( par()$oma )
  OuterMargins <- substr(OuterMargins, 5, nchar(OuterMargins))
  OuterMargins <-
       paste("oma = c(", gsub(" ",",",OuterMargins), ")", sep="")
  mtext(OuterMargins,
        NORTH<-3, line=0.4, adj=0.0, cex=1.5, col="green", outer=TRUE)
}

# Figure 2B.  Additional Annotations.
Figure2B <- function()
{
  Figure2A()

  # Text:  one string per mtext call
  mtext("mar[SOUTH<-1]", SOUTH<-1, at=5, line=2,    cex=1.2, col="blue")
  mtext("mar[WEST<-2]",  WEST <-2, at=5, line=2,    cex=1.2, col="blue")
  mtext("mar[NORTH<-3]", NORTH<-3, at=5, line=0.25, cex=1.2, col="blue")
  mtext("mar[EAST<-4]",  EAST <-4, at=5, line=0.25, cex=1.2, col="blue")

  # Text:  vector of strings in mtext call
  mtext(c("oma[SOUTH<-1]", "oma[WEST<-2]", "oma[NORTH<-3]", "oma[EAST<-4]"),
        c(SOUTH<-1, WEST<-2, NORTH<-3, EAST<-4),
        line=0.4, cex=1.2, col="green", outer=TRUE)
}

Figure2B()
```