Visualising uncertainty in time-series using animations

May 13, 2017
By

(This article was first published on Bluecology blog, and kindly contributed to R-bloggers)

Visualising uncertainty in time-series using animations

I was reading an example of visualising election outcome uncertainty
using a wobbly
needle
and
wondered, could animating uncertainty help communicate for the kind of
work I do on population modelling?

One challenge when plotting uncertainty intervals for population trends
is that confidence intervals don’t capture the wandering nature of
population processes. The figure here shows the median of a poisson
random walk, with 90% quantiles.

In fact, each trend line will wander around, like the next image. So,
the median and intervals give a false impression that an individual
trend will experience relative stability.


A graph of the median and quantiles doesn’t represent the uncertainty in any single trend line very well.

I decided to try and animate the uncertainty, using R.

The first step was to try and install the R package animate.
Unfortunately I was in a world of Unix coding pain trying to get
animation to work. In particular to get R to talk to the
ImageMagick program. After cutting and
pasting numerous commands from StackOverflow to my terminal I gave up.
It seems I need to update my OS before I can get any further (from
Mavericks).

In the end I found I could use ImageMagick directly from the command
line. So now what we need is a folder full of .png files that we wish to
turn into the animation. That Unix code to do that is at the end of this
blog.

First we should set up a random walk process. I decided to go for a
Poisson random walk, because it will have integer numbers, like a real
population. I wrote a little function to achieve this

walkmod <- function(t, yinit, thetasd, rho){
  y <- numeric(t)
  y[1] <- yinit
  theta <- rnorm(t, sd = thetasd)
  for(i in 2:t){
    lambda <- (y[i-1]*rho + yinit*(1-rho))* exp(theta[i])
    y[i] <- rpois(1, lambda)
  }
  return(y)
}

The core equation is in the for loop. It just says the mean abundance
at a given time is dependent on abundance at the last time-step with
correlation rho plus some extra variation theta. We then sample the
observed abundance from a poisson distribution.

Now let’s use our function. We will load in the purrr package to help
plotting replicate runs of the random walk and the leaflet package
because it has some nice functions for making colour scales.

library(purrr)
library(leaflet)

t <- 50
yinit <- 10
thetasd <- 0.01
rho <- 0.9
n <- 30

y <- matrix(NA, nrow = t, ncol = n)
tmat <- matrix(rep(1:t), nrow = t, ncol = n)

for (i in 1:n){
    set.seed(i)
    y[,i] <- walkmod(t, yinit, thetasd, rho)
}
ymax <- max(y)
medy <- apply(y, 1 , median)

We just looped over the walkmod function n times so we could
generate n random walks.

Now let’s build a colour scale. I would like to colour the lines by the
abundance quantile they end on. So we apply colorQuantile to all
abundances in the last time-step to get a red to blue (RdBu) colour
palette for the lines:

quants <- c(0, 0.25, 0.5, 0.75, 1)
yquant <- quantile(y[t,], probs = quants)
mycols <- colorQuantile("RdBu", domain = y[t,], probs = quants)

The trick to making a .gif animation is just to make a bunch of figures
and save them as .png files in the order you want them to animate (or if
you can get animate package to work within R you can export the .gif
directly without having to save the .pngs).

So we just loop over walkmod plotting our lines each time. Note I use
walk from the purrr package to redraw the background lines in grey
each time.

thispath <- "/insert/your/path/here/"

for (i in 1:n){
    fname <- paste(thispath, i, ".png", sep = "")
    png(filename =fname, width = 480, height = 480)
    plot(0, 0, xlim = c(0, t), ylim = c(0, ymax),
         bty = 'n', type = "n", las = 1,
         ylab = "Abundance", xlab = "Time")
    lines(1:t, medy, col = "grey90", lwd = 2)
    walk(1:i, ~lines(tmat[,.x], y[,.x], lwd = 0.7, col = grey(0.5, 0.2)))
  lines(1:t, y[,i], lwd = 2, col = mycols(y[t, i]))
    legend("topleft", legend = quants,
        lwd = 2, col = mycols(yquant),
    title = "Quantile t=100", cex = 0.8)
  dev.off()
}

If you want to check this from R, without saving the files, then turn of
the png() function that saves the files and insert Sys.sleep(0.1) at
the end of the loop, so it will print each figure in turn but pause for
1/10 th of second each time.

We have one final figure to plot, that is the median line. We can do
this outside the for loop:

fname <- paste(thispath, i+1, ".png",  sep = "")
png(filename =fname, width = 480, height = 480)

plot(0, 0, xlim = c(0, t), ylim = c(0, ymax),
     bty = 'n', type = "n", las = 1,
     ylab = "Abundance", xlab = "Time", main = "Median")
walk(1:i, ~lines(tmat[,.x], y[,.x], lwd = 0.5, col = grey(0.5, 0.2)))
lines(1:t, medy, lwd = 3, col = "purple")
dev.off()

cd "your/path/here"
magick *.png images.gif
convert images.gif -set delay 30 \( +clone -set delay 100 \) +swap +delete images_pause.gif

The first line creates the gif and the second line adds a pause of 30
1/100ths of a second after each image and a pause of 1 second at the end
of the loop.

To leave a comment for the author, please follow the link and comment on their blog: Bluecology blog.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.

Sponsors

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)