The capture of 'El Mochomo' and homicides in Jalisco

mochomo-jalisco.R

########################################################
# Author: Diego Valle-Jones
# Website: www.diegovalle.net
# Date Created: Tue May 29 17:03:32 2012
# Email: diegovalle at gmail.com
# Purpose: The capture of 'El Mochomo' increased homicides in Jalisco contrary
# to statements by the Mexican government
# Copyright (c) Diego Valle-Jones. All rights reserved

##homicide data for Jalisco
hom.jal <-structure(list(week = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 
13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 
29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 
45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 
77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 
93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 
107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 
120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 
133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 
146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 
159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 
172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 
185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 
198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 
211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 
224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 
237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 
250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260, 261, 262, 
263, 264, 265, 266, 267, 268, 269, 270, 271, 272, 273, 274, 275, 
276, 277, 278, 279, 280, 281, 282, 283, 284, 285, 286, 287, 288, 
289, 290, 291, 292, 293, 294, 295, 296, 297, 298, 299, 300, 301, 
302, 303, 304, 305, 306, 307, 308, 309, 310, 311, 312, 313, 314, 
315, 316, 317, 318, 319, 320, 321, 322, 323, 324, 325, 326, 327, 
328, 329, 330, 331, 332, 333, 334, 335, 336, 337, 338, 339, 340, 
341, 342, 343, 344, 345, 346, 347, 348, 349, 350, 351, 352, 353, 
354, 355, 356, 357, 358, 359, 360, 361, 362, 363, 364), V1 = c(8, 
13, 4, 13, 11, 5, 7, 6, 9, 14, 7, 8, 6, 8, 10, 12, 3, 5, 4, 6, 
11, 6, 10, 6, 14, 7, 0, 4, 4, 16, 4, 4, 11, 10, 6, 6, 7, 7, 6, 
11, 17, 10, 14, 7, 6, 13, 6, 1, 9, 3, 7, 15, 1, 9, 5, 5, 10, 
12, 7, 6, 4, 8, 9, 4, 8, 9, 14, 7, 18, 9, 4, 11, 10, 8, 8, 8, 
11, 3, 8, 5, 5, 4, 9, 6, 14, 7, 10, 12, 12, 7, 9, 5, 14, 10, 
9, 5, 4, 18, 8, 14, 8, 16, 5, 12, 12, 19, 12, 8, 11, 10, 7, 7, 
7, 7, 10, 5, 14, 13, 10, 10, 8, 12, 9, 11, 10, 4, 11, 5, 13, 
6, 11, 3, 11, 7, 6, 9, 8, 8, 14, 5, 7, 11, 8, 3, 10, 5, 8, 9, 
17, 7, 9, 15, 7, 14, 6, 4, 16, 6, 5, 4, 16, 10, 4, 4, 11, 2, 
5, 3, 6, 11, 6, 4, 11, 7, 10, 14, 5, 11, 13, 6, 8, 8, 8, 7, 8, 
4, 11, 15, 8, 7, 13, 11, 5, 8, 10, 11, 12, 5, 11, 6, 5, 9, 14, 
10, 9, 10, 9, 11, 14, 10, 7, 6, 6, 11, 6, 6, 7, 7, 17, 11, 7, 
6, 8, 9, 9, 16, 13, 12, 10, 6, 17, 10, 16, 13, 4, 14, 8, 14, 
16, 10, 5, 16, 9, 4, 11, 13, 7, 13, 18, 5, 14, 16, 6, 8, 8, 10, 
12, 11, 11, 20, 15, 18, 11, 10, 17, 5, 25, 4, 7, 15, 19, 9, 10, 
5, 8, 16, 9, 18, 10, 5, 12, 17, 8, 22, 11, 11, 9, 10, 12, 16, 
8, 20, 7, 11, 24, 19, 9, 12, 12, 12, 14, 20, 16, 18, 8, 24, 14, 
11, 7, 22, 17, 14, 14, 8, 9, 15, 11, 16, 16, 11, 13, 16, 15, 
21, 12, 8, 23, 12, 31, 17, 15, 17, 25, 23, 24, 14, 14, 25, 16, 
21, 27, 25, 20, 18, 18, 28, 22, 32, 19, 23, 23, 42, 26, 35, 11, 
40, 32, 20, 30, 30, 24, 29, 12, 19), date = structure(c(12422, 
12429, 12436, 12443, 12450, 12457, 12464, 12471, 12478, 12485, 
12492, 12499, 12506, 12513, 12520, 12527, 12534, 12541, 12548, 
12555, 12562, 12569, 12576, 12583, 12590, 12597, 12604, 12611, 
12618, 12625, 12632, 12639, 12646, 12653, 12660, 12667, 12674, 
12681, 12688, 12695, 12702, 12709, 12716, 12723, 12730, 12737, 
12744, 12751, 12758, 12765, 12772, 12779, 12786, 12793, 12800, 
12807, 12814, 12821, 12828, 12835, 12842, 12849, 12856, 12863, 
12870, 12877, 12884, 12891, 12898, 12905, 12912, 12919, 12926, 
12933, 12940, 12947, 12954, 12961, 12968, 12975, 12982, 12989, 
12996, 13003, 13010, 13017, 13024, 13031, 13038, 13045, 13052, 
13059, 13066, 13073, 13080, 13087, 13094, 13101, 13108, 13115, 
13122, 13129, 13136, 13143, 13150, 13157, 13164, 13171, 13178, 
13185, 13192, 13199, 13206, 13213, 13220, 13227, 13234, 13241, 
13248, 13255, 13262, 13269, 13276, 13283, 13290, 13297, 13304, 
13311, 13318, 13325, 13332, 13339, 13346, 13353, 13360, 13367, 
13374, 13381, 13388, 13395, 13402, 13409, 13416, 13423, 13430, 
13437, 13444, 13451, 13458, 13465, 13472, 13479, 13486, 13493, 
13500, 13507, 13514, 13521, 13528, 13535, 13542, 13549, 13556, 
13563, 13570, 13577, 13584, 13591, 13598, 13605, 13612, 13619, 
13626, 13633, 13640, 13647, 13654, 13661, 13668, 13675, 13682, 
13689, 13696, 13703, 13710, 13717, 13724, 13731, 13738, 13745, 
13752, 13759, 13766, 13773, 13780, 13787, 13794, 13801, 13808, 
13815, 13822, 13829, 13836, 13843, 13850, 13857, 13864, 13871, 
13878, 13885, 13892, 13899, 13906, 13913, 13920, 13927, 13934, 
13941, 13948, 13955, 13962, 13969, 13976, 13983, 13990, 13997, 
14004, 14011, 14018, 14025, 14032, 14039, 14046, 14053, 14060, 
14067, 14074, 14081, 14088, 14095, 14102, 14109, 14116, 14123, 
14130, 14137, 14144, 14151, 14158, 14165, 14172, 14179, 14186, 
14193, 14200, 14207, 14214, 14221, 14228, 14235, 14242, 14249, 
14256, 14263, 14270, 14277, 14284, 14291, 14298, 14305, 14312, 
14319, 14326, 14333, 14340, 14347, 14354, 14361, 14368, 14375, 
14382, 14389, 14396, 14403, 14410, 14417, 14424, 14431, 14438, 
14445, 14452, 14459, 14466, 14473, 14480, 14487, 14494, 14501, 
14508, 14515, 14522, 14529, 14536, 14543, 14550, 14557, 14564, 
14571, 14578, 14585, 14592, 14599, 14606, 14613, 14620, 14627, 
14634, 14641, 14648, 14655, 14662, 14669, 14676, 14683, 14690, 
14697, 14704, 14711, 14718, 14725, 14732, 14739, 14746, 14753, 
14760, 14767, 14774, 14781, 14788, 14795, 14802, 14809, 14816, 
14823, 14830, 14837, 14844, 14851, 14858, 14865, 14872, 14879, 
14886, 14893, 14900, 14907, 14914, 14921, 14928, 14935, 14942, 
14949, 14956, 14963), class = "Date")), .Names = c("week", "V1", 
"date"), row.names = 2:365, class = "data.frame")
hom.jal$V1 <- log1p(hom.jal$V1)


library(xtable)
library(ggplot2)
library(dlm)
library(grid)
library(gridExtra)
library(strucchange)
library(reshape)
library(directlabels)
library(scales)
library(testthat)
library(Cairo)
options(xtable.type = 'html')


#Add text at the bottom of a plot
addSource <- function(plot, text = "Data Source: Mexican Mortality Database SINAIS/INEGI") {
  plot <- arrangeGrob(plot, 
              sub = textGrob(text,
                x = 0, hjust = -0.1, vjust=0.1,
                gp = gpar(fontface = "italic", fontsize = 9,
                  col = "gray50")))
  return(plot)
}

#Clean up the data for the structural change model
hom.jal$date <- as.Date(hom.jal$date)
jalisco <- hom.jal
jalisco$dhom <- c(NA, diff(jalisco$V1))
jalisco$lag <- c(NA, jalisco$V1[-359])

##AR(1) model
##h=20 == five month min span
bp.ri <- breakpoints(V1 ~ lag, h = 20, data = jalisco) 
confint(bp.ri)

#Week number of the capture of Mochomo and kidnapping of Nacho's son
nacho.kidnapping <- which(hom.jal$date == as.Date("2010-04-05"))
mochomo.capture <- which(hom.jal$date == as.Date("2008-01-28"))


#subset December 2010 because it's under-counted by around 20%
hom.jal <- subset(hom.jal, date < as.Date("2010-11-20"))
ts.jalisco <- ts(hom.jal$V1, start = c(2004,1), freq = 52)


##Create dummy variables indicating whether 'el mochomo' had been captured or whether the kidnapping had taken place 
hom.jal$group <- as.factor(c(rep("peace", mochomo.capture),   
  rep("mochomo", length(mochomo.capture:nacho.kidnapping)-1),
  rep("kidnapping", length(nacho.kidnapping:length(ts.jalisco))-1)))
x <- model.matrix(~-1 + group, hom.jal)

#Seasonal effects for the dlm
mod1 <- dlmModTrig(s = 12, dV = 5.1118, dW = 0) + dlmModPoly(1, dV = 0, dW = 81307e-3)
smoothTem1 <- dlmSmooth(ts.jalisco, mod1)
plot(ts(smoothTem1$s[2 : 13, c(1, 3, 5, 7, 9, 11)],
        names = paste("S", 1 : 6, sep = "_")),
     oma.multi = c(2, 0, 1, 0), pch = 16, nc = 1,
     yax.flip = TRUE, type = 
     'o', xlab = "", main = "")


#A braver man than me would add dlmModSeas(52)
modJaliscoReg <- dlmModReg(x) + dlmModTrig(s = 52, q = 1, dW = 0)#+dlmModSeas(52, dV=0)

buildFun <- function(theta) {
  V(modJaliscoReg) <- exp(theta[1])
  diag(W(modJaliscoReg))[1] <- exp(theta[2])
  return(modJaliscoReg)
}

#obtain MLE values 
fit <- dlmMLE(ts.jalisco, parm = rep(0, 2), build = buildFun)
#the component convergence is set to zero if and only if
#the algorithm converged successfully
expect_that(fit$convergence, equals(0))
#compute the model with the values obtained by MLE
modJaliscoReg <- buildFun(fit$par)


modSmooth <- dlmSmooth(ts.jalisco, mod = modJaliscoReg)

#Diagnostics of the dlm
res <- residuals(dlmFilter(ts.jalisco, mod = modJaliscoReg), sd=FALSE)
qqnorm(res)
abline(a=0,b=1)
tsdiag(dlmFilter(ts.jalisco, mod = modJaliscoReg))

#Error bands for the smoothed estimate
sesmooth <- dlmSvd2var(modSmooth$U.S,modSmooth$D.S)[-1]
width <- qnorm(.95) * sqrt(sapply(sesmooth,diag))

#Add the kalman smooth and the error bands to the original data frame
hom.jal$dlm <- ts(modSmooth$s[-1, 1] + (modSmooth$s[-1, 2] * x[,1]) + ((modSmooth$s[-1, 3] * x[,2]) + (modSmooth$s[-1, 4] * x[,3])), start = 2004)[1:359]
hom.jal$low <- hom.jal$dlm - as.data.frame(t(width))$V5
hom.jal$high <- hom.jal$dlm + as.data.frame(t(width))$V6

CairoSVG("jalisco.svg", h = 5, w = 8)
p <- ggplot(hom.jal, aes(date, expm1(V1))) +
  geom_point(size=1.2) +
  geom_line(aes(date, expm1(dlm)), color = "#3e6dfc") +
  scale_y_log10() +
  ylab("number of homicides") +
  xlab("date") +
  geom_ribbon(aes(ymin=expm1(low), ymax=expm1(high)), alpha = .2) +
  opts(title = "Weekly number of homicides in Jalisco,\nwith a smoothed estimate of homicides using a state space model") +
  theme_bw()  
  ## geom_segment(aes(yend = 20, xend =as.Date("2008-01-20"), 
  ##                               y=30, x=as.Date("2008-01-20")),
  ##              arrow=arrow(length=unit(0.25,"cm"))) +
  ## geom_segment(aes(yend = 25, xend =as.Date("2010-04-05"), 
  ##                               y=35, x=as.Date("2010-03-05")),
  ##              arrow=arrow(length=unit(0.25,"cm")))
addSource(p)
dev.off()

#Perception of Insecurity table
jalisco <- c(38,39,41,50)
colima <- c(19,10,27,48)
nayarit <- c(19,21,27,52)
ensi <- data.frame(Jalisco = jalisco,
Colima = colima,
Nayarit = nayarit)
ensi <- as.data.frame(t(ensi))
names(ensi) <- c("2005", "2008","2009", "2010")
xtable(ensi, digits = 0)
ensi <- melt(ensi)
ensi$year <- rep(c(2005, 2008,2009, 2010), each = 3)
ensi$state <- rep(c("Jalisco", "Colima", "Nayarit"), n = 3)
ensi$value <- ensi$value / 100

p <- ggplot(ensi, aes(year, value, group = state, color = state)) +
  geom_line() +
  geom_point(size = 2.5, alpha = .6) +
  scale_y_continuous(labels = percent, limits = c(0, .55)) +
  xlim(2005, 2010.7) +
  opts(title = "Percentage of the population that considers the municipality\nor borough they live in insecure, by state (2005, 2008-2010)") +
  xlab("year") +
  ylab("percent who feel insecure") +
  theme_bw()
p <- direct.label(p, "last.bumpup")
p <- addSource(p, text = "Data Source: ENSI-3, ENSI-5, ENSI-6, and ENSI-7")
ggsave("insecure.png", plot = p, dpi = 100, h = 5, w = 8)