#######################
### SETUP WORKSPACE ###
#######################

library(spatstat)
library(plyr)
library(lattice)

## Clear the workspace
rm(list=ls())

##############################################################
### NOTE: This line should be changed to point to the folder
###       where you have unzipped the contents of the
###       replication archive.
##############################################################

## Set working directory
setwd("C:/warren_jpr2015_rep")

#######################################
#This file will allow exact replication
#of the results reported in Table 1
#using the same randomly generated 
#quadrature points
#######################################
load("ppm_data.RData")
#######################################

#Alternatively, load new versions of
#data created through included scripts

# ## Execute load_data.R 
# source('load_data.R')
# 
# ## Execute krig_data.R 
# source('krig_data.R')
# 
# ## Execute ppm_data.R 
# source('ppm_data.R')
#
# ## Clear workspace
# rm(list=ls())
#
# ## Load newly created data from working directory
# load('ppm_data_new.RData')


###############################################
### Recode zeros for variables to-be logged ###
### (add minimum observed non-zero value    ###
###  of each variable, to all observations) ###
###############################################

logvars <- c("elev",
             "popdens",
             "gdp_ppp",
             "ged.pre5surv",
             "ged.pre2005")

for (v in logvars) {
  varname <- paste0(v, ".im")
  var <- get(varname)
  minval <- min(var)
  step <- min(var[eval.im(var > 0)])
  if (minval < 0) {
    assign(varname, eval.im(var + abs(minval) + step))
  }
  if (minval == 0) {
    assign(varname, eval.im(var + step))
  }
}
rm(logvars, v, var, varname, minval, step)


####################################
### Estimate point process model ###
####################################

### Define function to fit PPM model
ppm_estimate <- function(
  dv_version,           #temporal window for DV
  obs_rad,              #radius for observation window
  strauss_rad = 900000, #radius for Strauss process
  incl_ivs = FALSE,     #include media IVs
  incl_mult = FALSE,    #include rad*cell interaction
  incl_cont = FALSE,    #include control variables
  incl_lag = FALSE,     #include lagged DV
  incl_diff = FALSE,    #include change variables (Uganda)
  uganda = FALSE)       #estimate Uganda model
{

  #Select DV for given version and radius
  depvar <- get(paste0("ged.", dv_version, ".r", obs_rad))
  
  #Select quadrature scheme for given radius
  quad.pnts <- get(paste0("quad.r", obs_rad))
  
  #If estimating longitudinal Uganda model, 
  #use Uganda pre-post observation window
  if (uganda == TRUE) {    
    depvar <- depvar[dhs.uganda.owin]
    quad.pnts <- quad.ug[dhs.uganda.owin]
    strauss_rad <- 100000
  }
  
  #Define inter-point interaction function
  inter <- StraussHard(r=strauss_rad) 
  
  ##Build list of model terms
  
  #Control variables
  controls <- c("elec",  
                "log(elev)", 
                "log(gdp_ppp)", 
                "log(popdens)",
                "urban",
                "shrpower")
  
  #Media variables
  ivs <-  c("rad", "cell")
  
  #Multiplicative interaction term
  mult <- "rad*cell"
  
  #Longitudinal changes (Uganda only)
  diff <- c('rad_diff', 'cell_diff')
  
  varnames <- as.character()
  
  if (incl_ivs) {
    varnames <- c(varnames, ivs)
  }
  
  if (incl_mult) {
    varnames <- c(varnames, mult)
  }
  
  if (incl_diff) {
    varnames <- c(varnames, diff)
  }
  
  if (incl_cont) {
    varnames <- c(varnames, controls)
  }
  
  #Include lagged DV corresponding to DV version
  if (incl_lag) {  
    
    #If using the post-2005 DV window,
    #use violent events from the period 2000 - 2004
    #(see ppp_data.R for construction details)
    if (dv == "post2005") {
      varnames <- c(varnames,  "log(ged.pre2005)")
    }
    
    #If using the post-survey DV window,
    #use violent events from the five years 
    #preceding each survey
    #(see ppp_data.R for construction details)
    if (dv == "postsurv") {
      varnames <- c(varnames,  "log(ged.pre5surv)")
    }
  }
  
  #Count events within each country
  c_cnts <- count(depvar$marks$cowcode)
  
  #Include fixed effects for countries with more than 1 event
  if (uganda == FALSE) {
    cowcodes <- c_cnts$x[c_cnts$freq > 1]
    for (c in cowcodes) {
      varnames <- c(varnames, paste0("as.factor(cowcode == ", c, ")"))
    }
  }
  
  #Create model formula
  varnames.fmla <- as.formula(paste("~ ", paste(varnames, collapse= "+")))
  
  #Remove covariate info from DV points
  depvar <- unmark(depvar)
  
  #Create list of covariate surfaces
  vardata.list <- list()
  vars <- all.vars(varnames.fmla)
  for (var in vars) {
    vardata.list[[var]] <- get(paste0(var,'.im'))
  }  
  
  #Estimate model
  m <- ppm(quadscheme(depvar, quad.pnts), 
           varnames.fmla, 
           interaction=inter,  
           rbord=1,
           covariates=vardata.list,
           gcontrol=list(maxit=100))
}
###################################### 
##End ppm_estimate()

### TABLE 1: Main Results ###

### Set maximum distance from survey location (kilometers)
### (must be one of distances included in radius_list)

maxdist <- 40

### Set DV temporal window 
### (must be either "post2005" or "postsurv")
##
##    This corresponds to the two versions of the
##    dependent variable described in the paper.
##       "post2005" - temporal windows range from 2005 to 2010 
##       "postsurv" - temporal windows follow survey dates

dv <- "post2005"


###############
### MODEL 1 ###
###############
model1 <- ppm_estimate(dv, maxdist, 
                       incl_cont=TRUE,
                       incl_lag=TRUE)
coef(summary(model1))
AIC(model1)


###############
### MODEL 2 ###
###############
model2 <- ppm_estimate(dv, maxdist, 
                       incl_cont=TRUE, 
                       incl_lag=TRUE,
                       incl_ivs=TRUE)
coef(summary(model2))
AIC(model2)


###############
### MODEL 3 ###
###############
model3 <- ppm_estimate(dv, maxdist, 
                       incl_cont=TRUE, 
                       incl_lag=TRUE,
                       incl_ivs=TRUE, 
                       incl_mult=TRUE)
coef(summary(model3))
AIC(model3)


###############
### MODEL 4 ###
###############
dv <- "postsurv"
model4 <- ppm_estimate(dv, maxdist, 
                       incl_cont=TRUE, 
                       incl_lag=TRUE,
                       incl_ivs=TRUE, 
                       incl_mult=TRUE)
coef(summary(model4))
AIC(model4)



##################################################
### Robustness checks with observation windows ###
### constructed using alternative maximum      ###
### distances from survey locations            ###
##################################################

maxdist <- 30
model4a <- ppm_estimate(dv, maxdist, 
                       incl_cont=TRUE, 
                       incl_lag=TRUE,
                       incl_ivs=TRUE, 
                       incl_mult=TRUE)
coef(summary(model4a))
AIC(model4a)

maxdist <- 50
model4b <- ppm_estimate(dv, maxdist, 
                        incl_cont=TRUE, 
                        incl_lag=TRUE,
                        incl_ivs=TRUE, 
                        incl_mult=TRUE)
coef(summary(model4b))
AIC(model4b)



### TABLE 2: Uganda ###

###############
### MODEL 5 ###
###############
dv <- "postsurv"
maxdist <- 30
model5 <- ppm_estimate(dv, maxdist, 
                       incl_ivs=TRUE, 
                       uganda=TRUE)
coef(summary(model5))
AIC(model5)


###############
### MODEL 6 ###
###############
dv <- "postsurv"
maxdist <- 30
model6 <- ppm_estimate(dv, maxdist, 
                       incl_ivs=TRUE, 
                       incl_diff=TRUE, 
                       uganda=TRUE)
coef(summary(model6))
AIC(model6)



# NOTE: Due to a small error in the construction of the Uganda observation window, 
# not discovered before going to print, the results generated here for Models 5 and 6 
# differ slightly from those reported in the paper.  However, they remain substantively
# equivalent, with coefficients pointing in the same direction, and similar levels of 
# statistical significance.  For replication purposes, the new results are printed below, 
# which should match the output on screen.

# > coef(summary(model5))
# Estimate       S.E.      CI95.lo     CI95.hi Ztest       Zval
# (Intercept) -21.49753216 0.30053652 -22.08657292 -20.9084914   *** -71.530515
# rad          -3.87188150 0.68315954  -5.21084960  -2.5329134   ***  -5.667609
# cell          5.37495170 0.98159967   3.45105169   7.2988517   ***   5.475706
# Interaction   0.08183887 0.01484029   0.05275244   0.1109253   ***   5.514642
# > AIC(model5)
# [1] 2431.179

# > coef(summary(model6))
# Estimate       S.E.      CI95.lo     CI95.hi Ztest        Zval
# (Intercept) -21.42295485 0.30886535 -22.02831981 -20.8175899   *** -69.3601752
# rad          -2.86650000 1.35426246  -5.52080564  -0.2121944     *  -2.1166503
# cell          4.92714773 1.24990994   2.47736927   7.3769262   ***   3.9420022
# rad_diff      1.75985049 2.62913978  -3.39316879   6.9128698         0.6693636
# cell_diff    -1.62278728 2.07408109  -5.68791152   2.4423370        -0.7824126
# Interaction   0.07906789 0.02528933   0.02950171   0.1286341    **   3.1265317
# > AIC(model6)
# [1] 2434.554


#################################
### Generate Predicted Values ###
#################################

#Select key indepdent variables
xvar.name <- 'rad'
yvar.name <- 'cell'

#Find range of xvar values seen within obsevation window
xvar <- get(paste0(xvar.name, ".im"))
xvar <- na.omit(xvar[dhs.all.owin.r40])
xvar <- xvar[xvar >= 0]
xvar.min <- round(min(xvar), 3)
xvar.max <- round(max(xvar), 3)

#Break range into 100 steps
xvar.vals <- do.breaks(c(xvar.min, xvar.max), 100)

#Find range of yvar values seen within obsevation window
yvar <- get(paste0(yvar.name, ".im"))
yvar <- na.omit(yvar[dhs.all.owin.r40])
yvar <- yvar[yvar >= 0]
yvar.min <- round(min(yvar), 3)
yvar.max <- round(max(yvar), 3)

#Break range into 100 steps
yvar.vals <- do.breaks(c(yvar.min, yvar.max), 100)

#Generate matrix of prediction values
predvals <- expand.grid(xvar.name = xvar.vals,
                        yvar.name = yvar.vals)

names(predvals) <- c(xvar.name, 
                     yvar.name)

predvals$x <- predvals[[xvar.name]]
predvals$y <- predvals[[yvar.name]]

#Set values of other independent variables to their means
#(excluding key IVs and interaction terms)
varnames.text <- all.vars(model3$trend)
varnames.means <-   varnames.text[which(
                     ! varnames.text  == xvar.name &
                     ! varnames.text  == yvar.name  )]

for (vname in varnames.means) {
  vname.im <- paste0(vname, ".im")
  var <- get(vname.im)
  var <- var[dhs.all.owin.r40]
  predvals[[vname]] <- mean(var, na.rm=TRUE)
} 

#Set categorical variables to most frequently observed value
calc_mode <- function(x) { 
  ux <- unique(x) 
  ux[which.max(tabulate(match(x, ux)))] 
}
predvals$cowcode <- calc_mode(na.omit(as.vector(cowcode.im$v)))


#Calculate predicted outcomes (conditional intensity function)
predvals$pred.cif <- predict.ppm(model3, 
                       locations=predvals, 
                       covariates=predvals,
                       type="cif")

#Shift scale from square meters to square kilometers
predvals$z <- predvals$pred.cif * 10^6

#Plot Predicted Values
library(RColorBrewer)
library(rasterVis)

#Convert conditional intensity function to raster
pred_raster <- rasterFromXYZ(cbind(predvals$x, predvals$y, predvals$z))

levelplot(pred_raster, xvar=xvar, yvar=yvar, zscaleLog=TRUE, margin=FALSE,
          par.settings=rasterTheme(region=brewer.pal(9, 'YlOrRd')),
          xlab="Radio Penetration",
          ylab="Cellular Penetration",
          main="Predicted Intensity of Violent Events")


