Simple plotting with R

cornflowerblue

It is important to watch one’s figure. To obsess over it, inspecting every detail, carefully erasing unpleasant angles, making sure all elements are perfectly aligned. The hours I spend caring about my figures may one day be my downfall, but until that time I shall continue to worry about things like inner margins, font embedding and what shade of cornflower blue to use.*

A few weeks ago** I decided to use R for some data analysis. It transpires that the R I thought I knew was in fact a semi-incomprehensible amalgamation of gnuplot and Octave, so I was forced to go back to the beginning. A lot of language-learning boils down to reading documentation (still looking for the official documentation for Mandarin, though) but that’s inefficient when you have a specific task to do. I’ve waded around in the help files, and here’s what I dredged up for basic (but not ugly) plotting. Note: as far as I can tell, ggplot2 is the main R package for plotting. This post is not about that. That comes later.

Basics

data<-read.table("file_path_goes_here", sep=",", as.is=TRUE)

This’ll suck your datafile into R. The sep option should contain whatever separates fields in your file. Default is whitespace. The as.is=TRUE option stops R from turning fields containing text into ‘factors’, which may not be what you want. Note: Try to name your variables informatively, unlike me.

Let’s suppose we want to plot some function of the fourth column of the data frame  (data[,4]) against the first column (data[,1]).  Note: if you wanted the fifth row for example, you would write data[5,]. The i,jth element is just data[i,j].

len<-length(data[,1]) == nrow(data)
p_col<-2*(1-pnorm(abs(data[,4])))
data<-cbind(data,abs_col)
chro<-17

The first line is figuring out how many rows are in our data frame in two equivalent ways. Note that length(data) just gives you the number of columns. The second line creates a new column whose values are a function of the (standard normal) cumulative distribution function (pnorm) of the absolute value of whatever is in the fourth column of the data frame. The third line appends (columnwise, with cbind) it to the data frame. Note: I think this might mess up the column names of the data-frame, so there must be a better way to do this. The last line is just defining a variable we’ll use later.

Now let’s use some other functions to create (useless***) vectors.

x<-seq(min(data[,1]), max(data[,1]), length.out=1000)
y<-rep(runif(100), 10)

x is just a vector of real numbers starting at min(data[,1) and ending with max(data[,1]), with evenly spaced intervals so its length is 1000 (seq creates a sequence). y is also a vector of length 1000. runif(100) creates a vector of 100 numbers between 0 and 1 taken from a uniform distribution. rep then repeats this 10 times. If you wanted the repetition to be elementwise (eg [a,a,b,b] instead of [a,b,a,b]) just write each=10.

Basic plot

plot(x, y, xlab="X Axis Name", ylab="Y Axis Name", main="Boring Plot Title")

So boring and default. How about

hist(y, freq=FALSE, main="Hist takes similar options")

Surprisingly enough this makes a histogram. freq=FALSE makes the y-axis into densities (eg fractions) rather than frequencies.

Cool plots

Now the basics are covered, let’s go crazy and make something that looks like this:

A collection of meaningless plots.
A collection of meaningless plots.

layout(matrix(c(1,1,2,3), 2, 2, byrow=TRUE), widths=c(1,3),heights=c(1.75,1))

We’re creating the plot environment here. What is going on is as follows: we’ve got three plots to arrange. Each one of these corresponds to a number in the matrix, which determines where it will go. Note: you can clearly try to give nonsensical layouts here, like demanding the first plot exist in both upper left and bottom right corners. This produces… interesting results.

The matrix command takes a list (c(1,1,2,3)) and matrix…es it. Here we’re telling it to have 2 rows and 2 columns, and to fill the matrix hull with the elements of the list row by row (byrow=TRUE – the default is by column).

We then inform layout that we want the columns of the matrix to have relative widths c(1,3) (note the second plot on the second line is 3 times as wide as the first) and the rows to have relative heights c(1.75,1). Many joyous minutes can be spent tweaking these values.

par(family="Avenir",bty="l")

par lets you globally set graphical parameters. A lot of its options are the same as for plot or hist, but it saves you having to write them for every individual plot. Here we’re setting the font family and bty, which determines the box style for the plots (eg the presence of axis lines). The options for these are depictions of the results: L, 7, C, U, ], etc. Note: I love Avenir, but using ‘strange’ fonts like this can cause problems if/when you go to save these plots directly to pdf.

I thoroughly recommend perusing the options of par. They will open your eyes to the scope of basic R plotting. In particular mar and omi lend themselves well to tweaking, and mfrow is like a poor-man’s layout.

plot(data[,1]/(1e07), data[,4], type="h", col=ifelse(abs(data[,4])>2.5, "darkorange1","black"), ylim=c(1.03*min(data[,4], na.rm=TRUE), 1.03*max(data[,4], na.rm=TRUE)), main=paste("Some data on chromosome", chro), xlab="Physical Distance (e+07)", ylab="Test Statistic")

abline(h=c(-2.5,2.5),col="darkorange1")

Right, so. What’s going on here:

  • I’m plotting data[,4] against data[,1]/(1e07). In my case, the elements of data[,1] are all of the order 10 million, so the division just rescales the x-axis to make it easier on the eye. R understands scientific notation like this, which is nice.
  • type: "p" plots points, "h" plots vertical bars (as in a histogram).
  • pch: (see later) is point style (if we had them). Tiny dots are ".", bigger dots are "20", see a diagram for more options.
  • col: (colour) you could write something simple like col="cornflowerblue" to make all of your points (or vertical lines as the case is here) that colour. But that’s boring. What we have here is "darkorange1" if abs(data[,4])>2.5, and "black" otherwise. Note: you may have noticed that my data[,4] column consists of z-scores, and I’m highlighting the ‘unusual’ ones.
  • ylim: we’re manually setting the limits on the y axis here, so it’s ylim = c(ymin, ymax). I’m going a little bit beyond the actual range of the data, and making sure to remove “nan“s in the min/max functions. Note: if your data is nice unlike mine, you shouldn’t need to do this.
  • paste: lets me combine text with other variables in the environment. In this case, it’s chro, which we defined earlier.
  • abline: this just adds a line to the plot. In this case, because a y-position was specified, it’s horizontal (actually, I added two lines, because I specified two positions). Alternately, abline(v=a) for a vertical line.

Let’s give it a legend.

legend(0,1.05*max(data[,4],na.rm=TRUE),"Orange : |Score| >2.5",bg="cornflowerblue",cex=0.9)

The first two options given to legend are the location of the upper-left corner of its box. cex is an overall scaling factor. See the options to par again.

I told layout to expect 3 plots, so this goes in the first place. I need to supply two more plots. The first one will be a histogram of our allegedly-normally-distributed dataset, column 4.

hist(data[,4],breaks=50,freq=FALSE,main="Scores",xlab="Test Statistic")
z<-seq(-4,4,length=200)
lines(z,dnorm(z),col="darkorange1")

The second two lines produce a normal distribution and overlay them onto the hist plot. (Both points and lines will do this, if you need to add more data to an existing plot). breaks is a histogram option saying how many individual bars we want, so you can adjust the ‘resolution’ of the histogram (so tweakable!). dnorm gives the probability density function for a standard normal, which we compare on the histogram.

And now for the final plot. Assuming the values in data[,4] come from a random variable with (standard) normal distribution, then p_col = 2*(1-pnorm(abs(data[,4])) is the probability of getting a more extreme value. (pnorm(abs(data[,4])) is just the cumulative distribution function of the absolute value of the variable, so 1-pnorm is the probability of getting a higher value. The multiplication by 2 compensates for our having taken an absolute value). We can then take the -log (base 10) of this, as is the style. Plotting is simple enough:

plot(data[,1]/(1e07), -log10(p_col), type="p", pch=20, cex=0.65 col=ifelse(p_col<0.05, "darkorange1","black"), ylim=c(0,max(-log10(p_col),na.rm=TRUE)), main=paste("Homebrew Manhattan plot for chromosome", chro), xlab="Physical Distance (e+07)", ylab="-log10(p)")

Here I’ve used cex to make the plot points a bit smaller.

Statistical note: depending on the nature of the analysis you’re doing, a ‘p-value’ of 0.05 may not be significant. You may be doing multiple hypothesis testing and will need to account for this. The example shown here is purely for demonstrative purposes and has no scientific legitimacy.

And we are done! For this particular example, I would be inclined to spend additional time tweaking the various inner/outer/inter-plot margins to reduce the amount of whitespace in the figure, but that’s just me.

*technical point: R only has one shade of cornflower blue.

**It seems I started this post over a month ago. Oops!

***I was originally planning on including the full complexity of my plotting script here, which involved vectors like this. Now they serve only as examples of what is possible.

Advertisements

Sketch – Barbarian Warrior

I have recently started drawing again. Here is my first proper attempt in approximately four years.

This is pretty basic, drawn from reference: http://mjranum-stock.deviantart.com/art/Barbarian-Warrior-27-77829696. I didn’t focus much on realism when I was growing up, and I almost never drew male figures so this is an attempt to remedy that.

I can’t say how long it took because I was doing it on and off for a few hours.

Things I learned:

  • Feet are hard. Must not become Rob Liefeld and deal with this by avoidance.
  • I don’t know how to draw realistic-looking hair.
  • Everything in this drawing is approximately the same tone despite this not being the case in the reference. This makes the drawing somewhat flat. When I started the sketch the plan was to do it very quickly, however, so I didn’t want to spend ages carefully shading. It might have been better to roughly shade in the dark areas first, and then add detail before I got bored.