Home > Uncategorized > Summary plots

Summary plots

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: ,
  1. Duleep Kumar Samuel
    August 3, 2010 at 4:24 am

    Thanks really useful and beautiful

  2. January 28, 2011 at 3:19 am

    Thanks for this code it looks really nice. I’m going to make a widget out of it!!

  3. Michael Bishop
    December 2, 2011 at 9:57 pm

    Hi, I should have mentioned this earlier, but I shared this code (slightly edited) on CrossValidated: http://stats.stackexchange.com/a/10447/3748

    • December 3, 2011 at 11:19 pm

      Thanks – glad you’ve found it useful. Thanks for adding the code to remove the NAs.

  1. No trackbacks yet.

Leave a comment