Archive

Archive for August, 2010

An HSV colour wheel in R

If you’ve read any of my previous posts, you’ll notice that they’re rather scanty on colour. There’s a reason for this. Mainly, that to get a good colour output takes some time. I recently read a commentary in Nature methods (sorry if you don’t have access to it, but this looks like it may be the first part of an interesting series of articles), which discusses colour in graphics. The author suggests a colour wheel, and I thought I’d have a go in R:

You have to click on it to read the text, sorry. There’s probably much easier ways to do it, and it takes a silly amount of time to render (several seconds! – all those nested loops), but this code below makes the colour wheel. If you set the variables t.hue, t.sat and t.val, the bottom right box is the resulting colour (the box just to the bottom right of the colour wheel is the hue with sat and val set to 1.0). Then on the right is the plot of val, and below is the plot of sat. As you go anti-clockwise from the x axis round your hue increases from 0.0 to 1.0.

So you can play around with colour, see what works and what doesn’t. This uses the HSV approach, which seemed okay for my purposes. rgb2hsv() converts rgb into hsv (obviously), if you are more familiar with the RGB approach. There are lots of other resources for colour in R, one of my favourites is here, and of course you can always search R-bloggers.

## colour plot

require(graphics)
t.hue <- 0.65     ## this is the user entered hue, sat and value
t.sat <- 0.5
t.val <- 0.9
def.par <- par(no.readonly = TRUE)
layout( matrix(c(1,1,2,1,1,2,3,3,4), 3, 3, byrow = TRUE))

## prepare the plot for the wheel 
x <- (-100:100)*0.01
y <- (-100:100)*0.01
## blank plot to prepare the axis
plot(x,y, pch = 20, col = 0, bty = "n", xaxt = "n", yaxt = "n", ann = F) 

## make the wheel
for (x in (-100:100)*0.01){
  for (y in (-100:100)*0.01){
    theta <- atan2(y,x)     # theta is the angle
    hue <-  Mod(theta/(pi)) # make the hue dependent upon the angle 
    sat <- (x^2 + y^2)      # make the saturation depend upon distance from origin
    if (x^2 + y^2 <= 1){
       if (y > 0) {points(x,y, pch = 19, col = hsv(h = hue/2, s = sat, v = 1))}
       if (y < 0) {points(-x,y, pch = 19, col = hsv(h = hue/2 + 0.5, s = sat, v = 1))}
      }
    }
  }
legend("center", "hue", bty = "n")
text(0.9,0, labels = "0.0")
text(0,0.9, labels = "0.25")
text(-0.9,0, labels = "0.5")
text(0,-0.9, labels = "0.75") 
## bottom right colour box inset into wheel
for (x in (80:100)*0.01){
  for (y in (-80:-100)*0.01){
    points (x,y, pch = 19, col = hsv(t.hue, s = 1, v = 1))
    }
  }

## right sided v scale 
x <- (0:100)*0.01
y <- (0:100)*0.01
plot(x,y, pch = 20, col = 0, xaxt = "n", yaxt = "n", bty = "n", ann = F)
for (x in (50:100)*0.01){
  for (y in (0:100)*0.01){
    hue <-  t.hue
    sat <- 1
    points(x,y, pch = 19, col = hsv(h = hue, s = sat, v = y))
    }
  }
legend("topleft", "value", bty = "n")
arrows(0.0, t.val, 0.5, t.val,length = 0.01, angle = 20)

  ## bottom saturation scale 
x <- (0:100)*0.01
y <- (0:100)*0.01
plot(x,y, pch = 20, col = 0, xaxt = "n", yaxt = "n", bty = "n", ann = F)
for (x in (0:100)*0.01){
  for (y in (0:50)*0.01){
    hue <-  t.hue
    points(x,y, pch = 19, col = hsv(h = hue, s = x, v = 1))
    }
  }
legend("topleft", "saturation", bty = "n")
arrows(t.sat,1.0, t.sat, 0.5, length = 0.01, angle = 20)

## bottom right plot
x <- (0:100)*0.01
y <- (0:100)*0.01
plot(x,y, pch = 20, col = 0, xaxt = "n", yaxt = "n", bty = "n", ann = F)
for (x in (0:25)*0.01){
  for (y in (0:100)*0.01){    
    points(x,y, pch = 19, col = hsv(h = t.hue, s = t.sat, v = t.val))
    }
  }
legtr <- paste( "hue=", t.hue, sep = "")
legr  <- paste( "sat=", t.sat, sep = "")
legbr <- paste("val=", t.val, sep = "")
legend("topright", legtr, bty = "n")
legend("right", legr, bty = "n")
legend("bottomright", legbr, bty = "n")

## reset the graphics display to default
par(def.par)
Categories: Uncategorized Tags: ,

Summary plots

August 2, 2010 5 comments

So, when you first look at some data, it’s helpful to get a feel of it. One way to do this is to do a plot or two. I’ve found myself continuously doing the same series of plots for different datasets, so in the end I wrote this short code to put all the plots together as a time saving device. Not pretty, but gets the job done.

The output looks like this:

So on the top a histogram with a normal distribution plot. On the right a QQ normal plot, with an Anderson Darling p value. Then in the middle on the left is the same data put into different numbers of bins, to see how this affects the look of the data. And on the right, we pretend that each value is the next one in a time series with equal time intervals between readings, and plot these. Below this is the ACF and PACF plots.

Hope someone else finds this useful. If there’s easier ways to do this, let me know. To use the code – put your data into a text file as a series of numbers called data.txt in the working directory, and run this code:

## univariate data summary
require(nortest)
data <- as.numeric(scan ("data.txt"))
# first job is to save the graphics parameters currently used
def.par <- par(no.readonly = TRUE)
par("plt" = c(.2,.95,.2,.8))
layout( matrix(c(1,1,2,2,1,1,2,2,4,5,8,8,6,7,9,10,3,3,9,10), 5, 4, byrow = TRUE))

#histogram on the top left
h <- hist(data, breaks = "Sturges", plot = FALSE)
xfit<-seq(min(data),max(data),length=100)
yfit<-yfit<-dnorm(xfit,mean=mean(data),sd=sd(data))
yfit <- yfit*diff(h$mids[1:2])*length(data)
plot (h, axes = TRUE, main = "Sturges")
lines(xfit, yfit, col="blue", lwd=2)
leg1 <- paste("mean = ", round(mean(data), digits = 4))
leg2 <- paste("sd = ", round(sd(data),digits = 4)) 
legend(x = "topright", c(leg1,leg2), bty = "n")

## normal qq plot
qqnorm(data, bty = "n", pch = 20)
qqline(data)
p <- ad.test(data)
leg <- paste("Anderson-Darling p = ", round(as.numeric(p[2]), digits = 4))
legend(x = "topleft", leg, bty = "n")

## boxplot (bottom left)
boxplot(data, horizontal = TRUE)
leg1 <- paste("median = ", round(median(data), digits = 4))
lq <- quantile(data, 0.25)
leg2 <- paste("25th quantile =  ", round(lq,digits = 4)) 
uq <- quantile(data, 0.75)
leg3 <- paste("75th quantile = ", round(uq,digits = 4)) 
legend(x = "top", leg1, bty = "n")
legend(x = "bottom", paste(leg2, leg3, sep = "; "), bty = "n")


## the various histograms with different bins
h2 <- hist(data,  breaks = (0:12 * (max(data) - min (data))/12)+min(data), plot = FALSE)
plot (h2, axes = TRUE, main = "12 bins")

h3 <- hist(data,  breaks = (0:10 * (max(data) - min (data))/10)+min(data), plot = FALSE)
plot (h3, axes = TRUE, main = "10 bins")
 
h4 <- hist(data,  breaks = (0:8 * (max(data) - min (data))/8)+min(data), plot = FALSE)
plot (h4, axes = TRUE, main = "8 bins")

h5 <- hist(data,  breaks = (0:6 * (max(data) - min (data))/6)+min(data), plot = FALSE)
plot (h5, axes = TRUE,main = "6 bins")

## the time series, ACF and PACF
plot (data, main = "Time series", pch = 20)
acf(data, lag.max = 20)
pacf(data, lag.max = 20)

## reset the graphics display to default
par(def.par)
Categories: Uncategorized Tags: ,