Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

## 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 <- 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.