Thursday, January 31, 2013

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:

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

# 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 = table(gun.deaths$date)
    # Calculate Cumulative Dead
    cum.deaths =
    for (i in 1:(length(cum.deaths)-1)) cum.deaths[i+1] = cum.deaths[i][i+1]

  plot(, 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(

  # 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([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[$age)] <- 999="" p="">
  gun.deaths$ = ""
  gun.deaths$[gun.deaths$age<12 children="" p="">  gun.deaths$[gun.deaths$age>11 & gun.deaths$age<20 p="" teens="">  gun.deaths$[gun.deaths$age>=20 & gun.deaths$age<30 adults20="" p="">  gun.deaths$[gun.deaths$age>=30 & gun.deaths$age<40 adults30="" p="">  gun.deaths$[gun.deaths$age>=40 & gun.deaths$age<65 madults="" p="">  gun.deaths$[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$"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$ == "children",]
  teens            = gun.deaths[gun.deaths$ == "teens",]
  adults20         = gun.deaths[gun.deaths$ == "adults20",]
  adults30         = gun.deaths[gun.deaths$ == "adults30",]
  madults          = gun.deaths[gun.deaths$ == "madults",]
  oadults          = gun.deaths[gun.deaths$ == "oadults",]

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

# This counts the total deaths = 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

# Rescale the US map to fit in our data = map('state', plot=F)$x = ((map('state', plot=F)$x+125)/(60))*ndays$y = ((map('state', plot=F)$y-24)/(27)*ymax)

# Static Plot - Final image
i=ndays, 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,
           "   ", toString([i]), "Killed",
           "  ", toString(cum.deaths[i]), "Dead"))
  grid(lwd = 2) # grid only in y-direction

  # Draw the US state map
  lines(, 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

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([i],"day",toString(i), "   ", toString([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(, 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
# deaths.animation()

ani.options(ani.width=1200, ani.height=900)
# In order for the saveGIF function to work you need install Image Magic Display (

# Map Updated Jan-31-2012


No comments :

Post a Comment