Simulate and show an animation of correlated datasets

set.seed(2015-08-02)
rbivariate <- function(r, mean.x = 0, sd.x=1, mean.y=0, sd.y=1, n=1) {
   z1 <- rnorm(n)
   z2 <- rnorm(n)
   x <- sqrt(1-r^2)*sd.x*z1 + r*sd.x*z2 + mean.x
   y <- sd.y*z2 + mean.y
   return(list(x,y))
}
library(animation)
library(ggplot2)
library(gridExtra)
frames <- 200
saveGIF(
        replicate(frames, {
            r13 <- rbivariate(0.13, n=52)
            r23 <- rbivariate(0.23, n=52)
            r33 <- rbivariate(0.33, n=52)
            r41 <- rbivariate(0.41, n=52)

            p13 <- qplot(r13[[1]], r13[[2]], xlab="x", ylab="y", main="0.13") + coord_cartesian(xlim=c(-2,2), ylim=c(-2,2)) + stat_smooth(method="lm", se=FALSE)
            p23 <- qplot(r23[[1]], r23[[2]], xlab="x", ylab="y", main="0.23") + coord_cartesian(xlim=c(-2,2), ylim=c(-2,2)) + stat_smooth(method="lm", se=FALSE)
            p33 <- qplot(r33[[1]], r33[[2]], xlab="x", ylab="y", main="0.33") + coord_cartesian(xlim=c(-2,2), ylim=c(-2,2)) + stat_smooth(method="lm", se=FALSE)
            p41 <- qplot(r41[[1]], r41[[2]], xlab="x", ylab="y", main="0.41") + coord_cartesian(xlim=c(-2,2), ylim=c(-2,2)) + stat_smooth(method="lm", se=FALSE)

            grid.arrange(p13, p23, p33, p41)
           }),
    interval = 0.8, ani.width = 900, ani.height = 900,
    movie.name = "/home/gwern/yvain-correlates-visualized.gif")