Data Fun - Inspired by Darasaurus

After my recent post on Anscombe's Quartet in which I demonstrated how to efficiently adjust any data set to match mean, variance, correlation (x,y), as well as regression coefficients. Philip Waggoner tuned me onto Justin Matejka and George Fitzmaurice's Datasaurus R package/paper in which the authors demonstrate an alternative method of modifying existing data to fit a specified mean and variance for X and Y. Their method randomly applies small disturbances to individual observations to gradually move the data to match a target preference set.

Inspired by their use of point images which match a specific parameter set, I have done generated some of my own. For all of them X has a mean and variance of 1 and 11. While y has a mean and variance of 1 and 1. The correlation between X and Y is set to 0 which causes some distortion in the images. More on that in the post.

Figure 1: Shows a graph of 12 data sets each with 15 transitional data sets. The mean, variance, and correlations of X and Y are held constant throughout the sets and transitions.

Data Source

I generated the data myself using Mobilefish's upload photo and record clicks webapp. The source images are from images I found online.

The only slight trick to using the data generated by Mobilefish was that the y cooridates are typically tracked from the top of the page with software, yet most statistical graphing software plots with y starting from the bottom of the graph.

Raw Images

The raw data when plotted look like their source material..

New Images: Force Cor(X,Y)=0

When we force the correlation of X and Y to be zero certain point distributions become distorted.

For Bart and the Cat forcing cor(X,Y) has noticable distortions while for the flower minimal distortions seem to have been introduced.

New Images: Force Cor(X,Y)<>0

It gets even worse when we impose a constant correlation between X and Y. The following shows the distortions to the flower when we change b1, keeping Cor(X,Y) constant and fixing the Y plot limits.

Figure 8: Shows the effect on Var(Y) that changing b1, when all other factors are held constant.

Slight changes to the Anscombe-Generator Code

In order to generate graphs that had cor(X,Y)=0 I had to modify my previous code to allow variation in Y that was completely independent of X. The problem with my code was that if b1=1, my calculation used SSE = (b1^2 * var(X))*n in order to infer how large the variation in u needed to be (varianceu = (SSE/corXY^2 - SSE)/n). This backwards inference does not work if b1=0.

So, just for the special case of corXY=0 I have included an additional error term E which is helpful in the even that b1=0.

Summary

The thought of use points to make recognizable images had not occurred to me until I viewed Justin Matejka and George Fitzmaurice's Datasaurus work. I hope that in making a slightly more efficient distribution manipulator I will allow new and better datasets to be generated which will help students understand the importance of graphical exploration of their data.

Code

My code as well as the slight change to Anscombe's code can be found here. The standarized data sets can be found here. They are the ones that end with std.csv

Fun with Bordered Cubes

I am interested in generating 3D reasoning items in R. To this end I have adapted some of the awesome functions built in the rgl library to my ends. My new function is 'cube' and it takes position and automatically sizes itself as a 1x1x1 cube though this can be adjusted through the scale argument.

Find the code on github 

# Easy Bordered Cubes in R
 
library('rgl'); library('magrittr')
 
cube <- function(x=0,y=0,z=0, bordered=TRUE, 
                 filled = TRUE, lwd=2, scale=1,
                 fillcol = gray(.95),
                 bordercol ='black', ...) {
 
  mycube <- cube3d()
 
  # Reduce size to unit
  mycube$vb[4,] <- mycube$vb[4,]/scale*2
 
  for (i in 1:length(x)) {
    # Add cube border
    if (bordered) {
      bcube <- mycube
      bcube$material$lwd <- lwd
      bcube$material$front <- 'line'
      bcube$material$back <- 'line'
      bcube %>% translate3d(x[i], y[i], z[i]) %>% shade3d
    }
    # Add cube fill
    if (filled) {
      fcube <- mycube
      fcube$vb[4,] <- fcube$vb[4,]*1.01
      fcube$material$col <- fillcol
      fcube %>% translate3d(x[i], y[i], z[i]) %>% shade3d
    }
  }
}
 
clear3d()
cube(0,0,0)
cube(1,0,0, filled=F)
cube(-1,0,0, bordered=F)
movie3d(spin3d(axis=c(0,0,1), rpm=20), duration=2.95) 

  
# I mapped R using an excel spreadsheet which 
# translated Xs into 2D locations points
clear3d()
y <- c(1,1,1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,4,5,5,
       5,5,6,6,6,6,7,7,7,7,7,8,8,8,8,9,9,9,9,10,
       10,10,10,11,11,11,11,11,12,12,12,12,12)
 
x <- c(8,7,6,3,2,1,7,6,3,2,6,5,3,2,6,5,4,3,2,5,4,
       3,2,5,4,3,2,6,5,4,3,2,7,6,3,2,7,6,3,2,7,6,
       3,2,6,5,4,3,2,5,4,3,2,1)
 
z <- cummax(y)*.5
 
length(x)==length(y)
cube(x,y,z)
movie3d(spin3d(axis=c(0,0,1), rpm=20), duration=2.95) 

 # Let's see how sin and cos can work together
z <- seq(0,6,.1)
x <- sin(pi*z)*z
y <- cos(pi*z)*z
 
clear3d()
cube(x,y,z*2, scale=.75)
movie3d(spin3d(axis=c(0,0,1), rpm=20), duration=2.95) 

# Now let's look at some of the prototyping for
# my 3D reasoning items.
clear3d()
vlist <- c(0,0,0)
for (i in 1:15) {
  cube(vlist[1],vlist[2],vlist[3])
  step <- sample(1:3, 1)
  vlist[step] <- vlist[step]+(-1)^rbinom(1,1,.25)
}
rgl.material(shininess=1)
bg3d("white")
clear3d(type = "lights")
 
rgl.viewpoint(theta = 90, phi = 0, fov = 0)
rgl.snapshot("2014-09-09angle1.png") 

rgl.viewpoint(theta = 0, phi = 0, fov = 0)
rgl.snapshot("2014-09-09angle2.png") 

rgl.viewpoint(theta = 0, phi = 90, fov = 0)
rgl.snapshot("2014-09-09angle3.png")
 


rgl.light()
rgl.viewpoint(theta = 180, phi = 20, fov = 60)
rgl.snapshot("2014-09-09cubes3d1.png") 

rgl.viewpoint(theta = -20, phi = -20, fov = 60)
rgl.snapshot("2014-09-09cubes3d2.png")

Created by Pretty R at inside-R.org

Alpha testing shinyapps.io - first impressions

ShinyApps.io is a new server which is currently in alpha testing to host Shiny applications.  It is being designed by the RStudio team and provides some distinct features different from that of the ShinyApps.io is intended for larger applications and I am guessing commercial applications in the long run.  Right now it is only allowing users by invitation into their alpha program, but they are accepting applicants for beta testing.
spark.rstudio server which is intended primarily for pilot testing Shiny apps. 

The way the user interacts with the system is extremely powerful.  In the standard package for shiny it is extremely easy to experiment with applications under development by simply navigating to the directory of interest with the setwd() command then trying your app with the runApp() function in the shiny library.  The new server has its own functions including the new and extremely powerful deployApp() function which acts the same as the runApp() function except that it contacts the ShinyApps server and immediately sets up a connection and deploys your app.  This makes the entire process of developing shiny applications and deploying them much easier.

I hope that this new service will further increase the user base of shiny! I think an R based interface for generating graphs will provide an invaluable tool for teaching and demonstrating empirical analysis methods.

I have been playing around with Shiny's seeming animation functionality.  I have created an animation like interface simulating the bouncing of a ball.

http://econometricsbysimulation.shinyapps.io/bounce/

GitHub Source

Spatial Critter Swarming Simulation

# I am interested in how small bits of individualized instructions can create collective action.

# In this simulation I will give a single instruction to each individual in the swarm.

# Choose another individual who is not too close, then accelerate towards that individual.

# I also control momentum causing the previous movement and direction to only decay at a small rate.

# TO SEE Original Script


# Critters are initially distributed randomly on a 1 x 1 grid.

ncritters = 40

xypos = matrix(runif(ncritters*2),ncol=2)
plot(xypos, main="Critters are Initially Distributed Randomly"
          , xlab="X", ylab="Y")

# Now let's imagine that each critter has an ideal safe distance from each other critter.

safe.dist = .3

critter.speed = .001

# If another critter is not at that safe distance than the critter will move towards the closest nearby critter.

# Let's see how this works.

# First let's check how close each critter is to each other critter.
# We will accomplish this by going through each critter and checking how far away each other critter is.
distances = NULL
for (i in 1:ncritters) distances = rbind(distances, apply((xypos[i,]-t(xypos))^2,2,sum))

# In order to prevent critters from always chasing whatever is closest to them (and themselves) we drop anything which is closer than the safe.distance.
distances[abs(distances)closest =  matrix(1:ncritters, ncol=ncritters, nrow=ncritters)[apply(abs(distances), 1, order)[1,]]
  # The apply command will apply the order command to each row whiel the [1,] selects only the critter that is closes.

# Plot the
plot(xypos, xlab = "X", ylab = "Y")
for (i in 1:ncritters) arrows(x0=xypos[i,1], y0=xypos[i,2],
                              x1=xypos[closest,][i,1],
                              y1=xypos[closest,][i,2],
                              length=.1)

# This calculates the difference between the current position of each critter and that of the closest critter.
ab = xypos-xypos[closest,]

# To see how this is calculated, see my previous post simulating a werewolf attack.

# Now calculate the difference in the horizontal and vertical axes that the critters will move as a projection into the direction of the closest critter outside of the safe zone.
a.prime = critter.speed/(1 + (ab[,2]^2)/(ab[,1]^2))^.5
b.prime = (critter.speed^2-a.prime^2)^.5

# This corrects the movement to ensure that the critters are flying at each other rather than away from each other.
movement = cbind(a.prime * sign(ab[,2]), b.prime * sign(ab[,1]))
between = function(xy1,xy2,point) (point>xy1&pointxy2&pointmovement = movement*(-1)^between(xypos,xypos[closest,], xypos-movement)

# Set the new xypos
xypos1 = xypos+movement

points(xypos1, col="red")

# ------------------------------------------------------
# Let's turn this into an animation.

library(animation)

# loopnum = 100; ncritters=40; inertia = .5; show.grid=T; ani.pause=F; plot.fixed=F; plot.centered=F; brownian = F; arrow = T
flocking <- ani.pause="F," arrow="T)" brownian="F," function="" inertia=".5," loopnum="100," ncritters="40," p="" plot.centered="F," plot.fixed="F," show.grid="T,">
  # Generate xy initial positions.
  # xypos will hold the current critter position while
  # xypos0 will hold the position of the critters the previous time.
  xypos = xypos0 = matrix(runif(ncritters*2),ncol=2)-.5

  movement0 = 0

  # Loop though all of the loops.
  for (i in 1:loopnum) {

  # This specifies the range to be graphed.
  if (plot.fixed) rangex=rangey = -.5:.5
  if (!plot.fixed) {
    rangex = c(min(xypos[,1]), max(xypos[,1]))
    rangey = c(min(xypos[,2]), max(xypos[,2]))
  }

  #  This handles the grid size when
  if (plot.centered&!plot.fixed) {
    rangex=c(-max(abs(xypos[,1])), max(abs(xypos[,1])))
    rangey=c(-max(abs(xypos[,2])), max(abs(xypos[,2])))
  }

  # This centers the plot at the middle (0,0) if the plot width is also set to be fixed.
  if (plot.centered&plot.fixed) rangex=rangex-mean(rangex)
  if (plot.centered&plot.fixed) rangey=rangey-mean(rangey)

  # Draw critters
  plot(xypos, main="Swarming Animation", xlab="X", ylab="Y", axes=F, ylim=rangey, xlim=rangex, type="p")
  # Draw arrows.  The start of the arrows is the previous periods location.
  if (arrow&i>1) arrows(x0=xypos0[,1],y0=xypos0[,2],x1=xypos[,1],y1=xypos[,2], length = .1)

  # Show the grid in the background.
  if (show.grid) {
    abline(v=seq(-10,10,.1))
    abline(h=seq(-10,10,.1))
  }

  # Show the origin
  text(0,0, "(0,0)")

  # Calculate each critters distance from each other
  distances = NULL
  for (i in 1:ncritters) distances = rbind(distances, apply((xypos[i,]-t(xypos))^2,2,sum))

  # Drop those within the safe zone.
  distances[abs(distances)
  # This selects the critter closest to the selected critter.
  closest =  matrix(1:ncritters, ncol=ncritters, nrow=ncritters)[apply(abs(distances), 1, order)[1,]]

#   distances[as.apply(!apply(distances, 1, is.na),1,sum)==0,]=0

  # As done above
  ab = xypos-xypos[closest,]

  a.prime = critter.speed/(1 + (ab[,2]^2)/(ab[,1]^2))^.5
  b.prime = (critter.speed^2-a.prime^2)^.5

  movement = cbind(a.prime * sign(ab[,2]), b.prime * sign(ab[,1]))

  between = function(xy1,xy2,point) (point>xy1&pointxy2&point
  movement = movement*(-1)^between(xypos,xypos[closest,], xypos-movement)

  movement[is.na(movement)]=0

  movement0 = movement0*inertia + movement

  # This fancy dodad allows half of the change in movement to be due to random variation.
  if (brownian) movement0=movement0+matrix(rnorm(ncritters*2),ncol=2)*critter.speed/2

  # Set the previous round's xy position to be equal to the current round's.
  xypos0 = xypos

  # Update the current round's.
  xypos = xypos+movement0

  # This is only used in the event that the animate package is in use.
  if (ani.pause) ani.pause()
  }

}

# This generates a GIF animation demonstrating smoothly how these GIFs can be incorper
  ani.options(ani.width=400, ani.height=400, interval=.1)

# You must have imagemagick installed for this to work.
  saveGIF(flocking(300,100,.999, ani.pause=T), movie.name = "Swarming.gif", replace=T)

# Here are two different graphs generated by the previous command(though the one on the bottom uses 200 frames while the one on the top uses 300)
->

# Let's see how this works.
flocking()
flocking(400,100,.99)
flocking(400,100,.99, plot.fixed=T)

Generate GIS integrated animation:Gun deaths animation - R

# I wanted to briefly revisit a previous post in order to update the graphics.

# The code used in the previous post has also become more useful since the database of guns deaths has tripled since the original post (yay!).

# If you would like to generate your own graph, make sure to use the r code attached to the original post rather than the original post.

# We can see a steady decline in gun deaths.  I am not sure what is driving this exactly.  I suspect maybe a decrease in gun suicides due to the winter coming to an end is driving this though I think gun deaths will rise as murders start increasing during the summer months.

US Daily Gun Deaths R Animation - Sandy Hook

R script

# Listenning to NPR about gun deaths in the US got me thinking.

# Find the article here:
# http://www.slate.com/articles/news_and_politics/crime/2012/12/gun_death_tally_every_american_gun_death_since_newtown_sandy_hook_shooting.html

# Let's create an animation of gun deaths since Dec 14, 2012
gun.deaths <- getcsv.php="" gun-deaths="" http:="" p="" read.csv="" slate-interactives-prod.elasticbeanstalk.com="">
gun.deaths$victimID

# I first need to change the dates from days to days for later reference
  # First I will make a table of the dates wich will include a count
  deaths.tab = table(gun.deaths$date)
    # Calculate Cumulative Dead
    cum.deaths = deaths.tab
    for (i in 1:(length(cum.deaths)-1)) cum.deaths[i+1] = cum.deaths[i]+deaths.tab[i+1]

  plot(deaths.tab, main="Daily Total US Deaths by Gun", ylab="Death count")

->

  # Number of days in our data (constantly updated every time we run the code
  ndays = length(deaths.tab)

  # This complicated bit of code will force the dates which are currently factor variables into string variables
  gun.deaths$dates = t(data.frame(lapply(gun.deaths$date, as.character), stringsAsFactors=FALSE))

  # Create an empty factor to be filled
  gun.deaths$day = NA
  # This command loops through all of the days and checks if each individual entry in the data from is from that day.
  # If it is, then it assigns that day to the day entry.
  for (i in 1:ndays) gun.deaths$day[gun.deaths$dates==names(deaths.tab)[i]] = i

# We will cut the data into different age categories

# Some individuals have ages missing.  We will code these as category 0.
  gun.deaths$age[is.na(gun.deaths$age)] <- 999="" p="">
  gun.deaths$age.cat = ""
  gun.deaths$age.cat[gun.deaths$age<12 children="" p="">  gun.deaths$age.cat[gun.deaths$age>11 & gun.deaths$age<20 p="" teens="">  gun.deaths$age.cat[gun.deaths$age>=20 & gun.deaths$age<30 adults20="" p="">  gun.deaths$age.cat[gun.deaths$age>=30 & gun.deaths$age<40 adults30="" p="">  gun.deaths$age.cat[gun.deaths$age>=40 & gun.deaths$age<65 madults="" p="">  gun.deaths$age.cat[gun.deaths$age>65 & gun.deaths$age<999 nbsp="" oadults="" p="">
# Adjust the latitude and logitude variables to account for a rescaling of the graph later
# as well as some noise which will help identify when there are multiple deaths in the same city.
  nll = length(gun.deaths$lng)
  gun.deaths$x = ((gun.deaths$lng+125)/(60))*ndays+rnorm(nll)*.006
  # For the graph that will be produced the 20 year olds have the highest likelihood of death.
  # Thus they will provide the y upper limit.
    ymax = ceiling(sum(gun.deaths$age.cat=="adults20")/50)*50
  gun.deaths$y = ((gun.deaths$lat-24)/(27))*ymax+rnorm(nll)*.06

# Generate subsets of the data.
  children         = gun.deaths[gun.deaths$age.cat == "children",]
  teens            = gun.deaths[gun.deaths$age.cat == "teens",]
  adults20         = gun.deaths[gun.deaths$age.cat == "adults20",]
  adults30         = gun.deaths[gun.deaths$age.cat == "adults30",]
  madults          = gun.deaths[gun.deaths$age.cat == "madults",]
  oadults          = gun.deaths[gun.deaths$age.cat == "oadults",]

# This will count cumulative deaths by data subset
children.tab = table(children$date)
  cum.children = children.tab
  for (i in 1:(length(cum.children)-1)) cum.children[i+1] = cum.children[i]+children.tab[i+1]

teens.tab = table(teens$date)
  cum.teens = teens.tab
  for (i in 1:(length(cum.teens)-1)) cum.teens[i+1] = cum.teens[i]+teens.tab[i+1]

adults20.tab = table(adults20$date)
  cum.adults20 = adults20.tab
  for (i in 1:(length(cum.adults20)-1)) cum.adults20[i+1] = cum.adults20[i]+adults20.tab[i+1]

adults30.tab = table(adults30$date)
  cum.adults30 = adults30.tab
  for (i in 1:(length(cum.adults30)-1)) cum.adults30[i+1] = cum.adults30[i]+adults30.tab[i+1]

madults.tab = table(madults$date)
  cum.madults = madults.tab
  for (i in 1:(length(cum.madults)-1)) cum.madults[i+1] = cum.madults[i]+madults.tab[i+1]

oadults.tab = table(oadults$date)
  cum.oadults = oadults.tab
  for (i in 1:(length(cum.oadults)-1)) cum.oadults[i+1] = cum.oadults[i]+oadults.tab[i+1]

# This counts the total deaths
cum.total = cum.adults20+cum.adults30+cum.teens+cum.children+cum.madults+cum.oadults

# In order to get a map of the US we will need to install the library maps
  library(maps)

# Rescale the US map to fit in our data
  map.us = map('state', plot=F)
  map.us$x = ((map('state', plot=F)$x+125)/(60))*ndays
  map.us$y = ((map('state', plot=F)$y-24)/(27)*ymax)
999>65>40>30>20>12>->
# Static Plot - Final image
i=ndays
dev.new(width=15, height=10)
plot(cum.adults20, type="n", ylim=c(0,ymax),
     ylab="Cumulative Deaths by Age Group" ,
     main="US gun Deaths Since Dec 14, 2012",  cex.main=1.5,
     xlab=paste(names(deaths.tab)[i],"day",toString(i),
           "   ", toString(deaths.tab[i]), "Killed",
           "  ", toString(cum.deaths[i]), "Dead"))
   
  grid(lwd = 2) # grid only in y-direction

  # Draw the US state map
  lines(map.us, type="l")

  # Place dots for every death
  lines(adults20$x, adults20$y, type="p", col="blue",pch=16, cex=.5)
  lines(adults30$x, adults30$y, type="p", col="purple",pch=16, cex=.5)
  lines(madults$x, madults$y, type="p", col="orange",pch=16, cex=.5)
  lines(oadults$x, oadults$y, type="p", col=colors()[641],pch=16, cex=.5)
  lines(teens$x, teens$y, type="p", col="red",pch=16, cex=.5)
  lines(children$x, children$y, type="p", col=colors()[258],pch=16, cex=.5)

  lines(cum.teens,    type="l", col=gray(.9), lwd=2)
  lines(cum.children, type="l", col=gray(.9), lwd=2)
  lines(cum.adults20, type="l", col=gray(.9), lwd=2)
  lines(cum.adults30, type="l", col=gray(.9), lwd=2)
  lines(cum.madults,  type="l", col=gray(.9), lwd=2)
  lines(cum.oadults,  type="l", col=gray(.9), lwd=2)

  lines(cum.teens,    type="l", col="red", lwd=1, cex=.5)
  lines(cum.children, type="l", col=colors()[258], lwd=1, cex=.5)
  lines(cum.adults20, type="l", col="blue", lwd=1, cex=.5)
  lines(cum.adults30, type="l", col="purple", lwd=1, cex=.5)
  lines(cum.madults,  type="l", col="orange", lwd=1, cex=.5)
  lines(cum.oadults,  type="l", col=colors()[641], lwd=1, cex=.5)

  lines(c(ndays,ndays), c(-15,ymax-27), lwd=2, lty=2)

  legend("topright", "Today", cex=1.5, bty="n")

legend(0, ymax+20, expression(Children, Teens, "Adults 20s", "Adults 30s",
   "Adults 40-65", "Adults 65+"), lty = 1:1, col = c(colors()[258], "red",
   "blue", "purple", "orange", colors()[641]),  adj = c(0, 0.6), ncol = 6,
   cex=.75, lwd=2
   )

# Sequential Graphic - Animation

# Animation package must be installed
library(animation)

deaths.animation = function() {

for (i in c(1:ndays, rep(ndays,10))) {
  # Adding the rep(ndays,10)) causes the animation to wait for the equivalent of 10 days after ending.

plot(cum.adults20, type="n", ylim=c(0,ymax),
     ylab="Cumulative Deaths by Age Group" ,
     main="US gun Deaths Since Dec 14, 2012", cex.main=2,
     xlab=paste(names(deaths.tab)[i],"day",toString(i), "   ", toString(deaths.tab[i]), "Killed", "  ", toString(cum.deaths[i]), "Dead"))
   
  # Draw the US state map
  grid(lwd = 2) # grid only in y-direction

  # Place dots for every death
  lines(map.us, type="l")

  t.adults20 = adults20[adults20$day  t.adults30 = adults30[adults30$day  t.madults  = madults[madults$day  t.oadults = oadults[oadults$day  t.teens = teens[teens$day  t.children = children[children$day
  lines(t.adults20$x, t.adults20$y, type="p", col="blue",pch=16, cex=1.5)
  lines(t.adults30$x, t.adults30$y, type="p", col="purple",pch=16, cex=1.5)
  lines(t.madults$x, t.madults$y, type="p", col="orange",pch=16, cex=1.5)
  lines(t.oadults$x, t.oadults$y, type="p", col=colors()[641],pch=16, cex=1.5)
  lines(t.teens$x, t.teens$y, type="p", col="red",pch=16, cex=1.5)
  lines(t.children$x, t.children$y, type="p", col=colors()[258],pch=16, cex=1.5)

  p.adults20 = adults20[adults20$day==i,]
  p.adults30 = adults30[adults30$day==i,]
  p.madults  = madults[madults$day==i,]
  p.oadults = oadults[oadults$day==i,]
  p.teens = teens[teens$day==i,]
  p.children = children[children$day==i,]

  lines(p.adults20$x, p.adults20$y, type="p", col="blue",pch=1, cex=5)
  lines(p.adults30$x, p.adults30$y, type="p", col="purple",pch=1, cex=5)
  lines(p.madults$x, p.madults$y, type="p", col="orange",pch=1, cex=5)
  lines(p.oadults$x, p.oadults$y, type="p", col=colors()[641],pch=1, cex=5)
  lines(p.teens$x, p.teens$y, type="p", col="red",pch=1, cex=5)
  lines(p.children$x, p.children$y, type="p", col=colors()[258],pch=1, cex=5)

  lines(cum.teens[1:i],    type="b", col=gray(.9), lwd=7)
  lines(cum.children[1:i], type="b", col=gray(.9), lwd=7)
  lines(cum.adults20[1:i], type="b", col=gray(.9), lwd=7)
  lines(cum.adults30[1:i], type="b", col=gray(.9), lwd=7)
  lines(cum.madults[1:i],  type="b", col=gray(.9), lwd=7)
  lines(cum.oadults[1:i],  type="b", col=gray(.9), lwd=7)

  lines(cum.teens[1:i],    type="b", col="red", lwd=2)
  lines(cum.children[1:i], type="b", col=colors()[258], lwd=2)
  lines(cum.adults20[1:i], type="b", col="blue", lwd=2)
  lines(cum.adults30[1:i], type="b", col="purple", lwd=2)
  lines(cum.madults[1:i],  type="b", col="orange", lwd=2)
  lines(cum.oadults[1:i],  type="b", col=colors()[641], lwd=2)


  lines(c(ndays,ndays), c(-15,ymax-27), lwd=2, lty=2)

  legend("topright", "Today", bty="n", cex=2.5)

legend(0, ymax+20, expression(Children, Teens, "Adults 20s", "Adults 30s",
   "Adults 40-65", "Adults 65+"), lty = 1:1, col = c(colors()[258], "red",
   "blue", "purple", "orange", colors()[641]),  adj = c(0, 0.6), ncol = 6,
   cex=1.5, lwd=1
   )
 
  ani.pause()
}
}
# deaths.animation()

ani.options(ani.width=1200, ani.height=900)
saveGIF(deaths.animation())
# In order for the saveGIF function to work you need install Image Magic Display (http://www.imagemagick.org/script/index.php)

# Map Updated Jan-31-2012

(OLD MAP BELOW)