Purpose:
Explore the geography of the Iraq conflict as represented by the WikiLeaks data
Describe and explain that patterns in event types
Research Questions:
1. For the 9 event types, describe the pattern in relation to CSR
2. Are the events distributed in relationship to one another?
3. Is it possible to predict the number of events of one type based upon the density surface of another?
####################################
####################################
#################################### Part 0: Load libraries, refine AOI,
#################################### and prepare data for PPA:
# install.packages('rgdal') install.packages('ggmap')
library(spatstat)
## Loading required package: mgcv
## This is mgcv 1.7-22. For overview type 'help("mgcv-package")'.
## Loading required package: deldir
## deldir 0.0-21
## spatstat 1.30-0 Type 'help(spatstat)' for an overview of spatstat
## 'latest.news()' for news on latest version 'licence.polygons()' for
## licence information on polygon calculations
library(sp)
library(rgdal)
## rgdal: version: 0.8.1, (SVN revision 415) Geospatial Data Abstraction
## Library extensions to R successfully loaded Loaded GDAL runtime: GDAL
## 1.9.2, released 2012/10/08 Path to GDAL shared files: C:/Program
## Files/R/R-2.15.2/library/rgdal/gdal GDAL does not use iconv for recoding
## strings. Loaded PROJ.4 runtime: Rel. 4.7.1, 23 September 2009,
## [PJ_VERSION: 470] Path to PROJ.4 shared files: C:/Program
## Files/R/R-2.15.2/library/rgdal/proj
library(maptools)
## Loading required package: foreign
## Loading required package: grid
## Loading required package: lattice
## Checking rgeos availability: TRUE
library(ggmap)
## Warning: package 'ggmap' was built under R version 2.15.3
## Loading required package: ggplot2
deaths <- read.csv("E:/Quant/Assignments/Lab04/Deaths only.csv") #Load the data
names(deaths)
## [1] "Report.key" "Date.and.time"
## [3] "YEAR" "Type"
## [5] "Category" "Title"
## [7] "Region" "Attack.on"
## [9] "Coalition.forces.wounded" "Coalition.forces.killed"
## [11] "Iraq.forces.wounded" "Iraq.forces.killed"
## [13] "Civilian.wia" "Civilian.kia"
## [15] "Enemy.wia" "Enemy.kia"
## [17] "Enemy.detained" "Total.deaths"
## [19] "Latitude" "Longitude"
# Create a spatial points data frame and project it to meters so distances
# are reasonable: First, we must clean up the data a bit:
deaths$Type <- sub(deaths$Type, pattern = "CRIMINAL EVENT", replacement = "Criminal Event")
deaths$Type <- sub(deaths$Type, pattern = "criminal event", replacement = "Criminal Event")
deaths$Type <- sub(deaths$Type, pattern = "EXPLOSIVE HAZARD", replacement = "Explosive Hazard")
deaths$Type <- factor(deaths$Type)
levels(deaths$Type) #Check that only 9 classes of Type exist
## [1] "Criminal Event" "Enemy Action" "Explosive Hazard"
## [4] "Friendly Action" "Friendly Fire" "Non-Combat Event"
## [7] "Other" "Suspicious Incident" "Threat Report"
# Refine AOI: plot using ggmaps Plot a 'locator map'
map.Iraq <- get_map(location = "Iraq", maptype = "hybrid", zoom = 6) # Get map
ggmap(map.Iraq, extent = "panel") + geom_point(data = deaths, aes(x = deaths$Longitude,
y = deaths$Latitude, size = Total.deaths), alpha = 0.3, color = "red") +
ggtitle("WikiLeaks Iraq Data:\nIncidents involving Death(s)")
# Plot a refined AOI
map.Baghdad <- get_map(location = "Baghdad", maptype = "hybrid", zoom = 12)
ggmap(map.Baghdad, extent = "panel") + geom_point(data = deaths, aes(x = deaths$Longitude,
y = deaths$Latitude, size = Total.deaths), alpha = 0.3, color = "red") +
ggtitle("Incidents involving Death(s)\nin Baghdad near AOI")
## Warning: Removed 35539 rows containing missing values (geom_point).
# Note that this map doesn't show the exact AOI we are using. See the
# attached map on the lab writeup instead.
# Subset the data:
deaths_AOI <- deaths[(deaths$Longitude < 44.465 & deaths$Longitude > 44.33) &
(deaths$Latitude > 33.235 & deaths$Latitude < 33.33), ]
deaths_marks <- deaths_AOI[, c(4, 7, 8, 18)]
names(deaths_marks)
## [1] "Type" "Region" "Attack.on" "Total.deaths"
# Now, create the Spatial points data frame
deaths_df <- SpatialPointsDataFrame(coords = deaths_AOI[, 19:20], proj4string = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"),
data = deaths_marks)
# Project the data:
deaths_df <- spTransform(deaths_df, CRS("+init=epsg:3839")) #Convert coords to meters
deaths_df$Easting <- coordinates(deaths_df)[, 2]
deaths_df$Northing <- coordinates(deaths_df)[, 1]
names(deaths_df)
## [1] "Type" "Region" "Attack.on" "Total.deaths"
## [5] "Easting" "Northing"
plot(x = deaths_df$Easting, y = deaths_df$Northing, asp = 1) #asp keeps the aspect ration at 1.0 (no distortion)
####################################
####################################
#################################### Question 1: For the 9 event types,
#################################### describe the distribution in relation
#################################### to CSR
# Start making a ppp object:
convexHull <- convexhull.xy(x = deaths_df$Easting, y = deaths_df$Northing)
plot(convexHull)
deaths_ppp <- ppp(x = deaths_df$Easting, y = deaths_df$Northing, window = convexHull,
marks = deaths_marks)
## Warning: data contain duplicated points
# Map it:
plot(deaths_ppp)
## Criminal Event Enemy Action Explosive Hazard
## 1 2 3
## Friendly Action Friendly Fire Non-Combat Event
## 4 5 6
## Other Suspicious Incident Threat Report
## 7 8 9
plot(density(deaths_ppp, sigma = 500))
plot(split(deaths_ppp)) #Create point maps
plot(density(split(deaths_ppp), sigma = 500)) #Kernel density map
# Weighted density map:
deaths_simp <- unmark(deaths_ppp) #remove all data from the points
deaths_mks <- marks(deaths_ppp)["Total.deaths"] #create a list of # of deaths
death_simp_ppp <- deaths_simp %mark% deaths_mks #make a point pattern (ppp) with a single variable
plot(density(death_simp_ppp, sigma = 500, weights = marks(death_simp_ppp)))
# Do some hypothesis testing using chi-squared:
qTest <- quadrat.test(deaths_ppp, nx = 5, ny = 5)
# plot(qTest)
qTest
##
## Chi-squared test of CSR using quadrat counts
##
## data: deaths_ppp
## X-squared = 5930, df = 24, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
## Quadrats: 25 tiles (levels of a pixel image)
# Look at clustering using the K-function
e_ppp <- ppp(x = deaths_df$Easting, y = deaths_df$Northing, window = convexHull) #The K-function doesn't like to have the marks...
## Warning: data contain duplicated points
plot(envelope(e_ppp, Kest, nsim = 99, correction = "border"))
## Generating 99 simulations of CSR ...
## 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.
##
## Done.
## lty col key label
## obs 1 1 obs K[obs](r)
## theo 2 2 theo K[theo](r)
## hi 1 8 hi K[hi](r)
## lo 1 8 lo K[lo](r)
## meaning
## obs observed value of K(r) for data pattern
## theo theoretical value of K(r) for CSR
## hi upper pointwise envelope of K(r) from simulations
## lo lower pointwise envelope of K(r) from simulations
# Get the quadrat Test and K-Functions for all type classes using a FOR
# loop:
for (t in 1:length(levels(deaths$Type))) {
cat(levels(deaths$Type)[t], "\n")
tempdat <- deaths_df[deaths_df$Type == levels(deaths$Type)[t], ]
temp_ppp <- ppp(x = tempdat$Easting, y = tempdat$Northing, window = convexHull,
marks = tempdat@data)
temp_q <- quadrat.test(temp_ppp, nx = 5, ny = 5)
cat("Q:", temp_q$statistic, "\n")
cat("P-value:", temp_q$p.value, "\n**************\n")
remove(temp_ppp)
remove(temp_q)
temp_ppp <- ppp(x = tempdat$Easting, y = tempdat$Northing, window = convexHull)
plot(envelope(temp_ppp, Kest, nsim = 99, correction = "border"))
remove(temp_ppp)
remove(tempdat)
}
## Criminal Event
## Warning: data contain duplicated points
## Q: 5291
## P-value: 0
## **************
## Warning: data contain duplicated points
## Generating 99 simulations of CSR ...
## 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.
##
## Done.
## Enemy Action
## Warning: data contain duplicated points
## Q: 888.8
## P-value: 6.834e-172
## **************
## Warning: data contain duplicated points
## Generating 99 simulations of CSR ...
## 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.
##
## Done.
## Explosive Hazard
## Warning: data contain duplicated points
## Q: 379.5
## P-value: 2.438e-65
## **************
## Warning: data contain duplicated points
## Generating 99 simulations of CSR ...
## 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.
##
## Done.
## Friendly Action
## Warning: data contain duplicated points
## Warning: Some expected counts are small; chi^2 approximation may be
## inaccurate
## Q: 198.3
## P-value: 4.43e-29
## **************
## Warning: data contain duplicated points
## Generating 99 simulations of CSR ...
## 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.
##
## Done.
## Friendly Fire
## Warning: Some expected counts are small; chi^2 approximation may be
## inaccurate
## Q: 21.64
## P-value: 0.7982
## **************
## Generating 99 simulations of CSR ...
## 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.
##
## Done.
## Non-Combat Event
## Warning: Some expected counts are small; chi^2 approximation may be
## inaccurate
## Q: 33.78
## P-value: 0.1773
## **************
## Generating 99 simulations of CSR ...
## 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.
##
## Done.
## Other
## Warning: Some expected counts are small; chi^2 approximation may be
## inaccurate
## Q: 24.16
## P-value: 0.9048
## **************
## Generating 99 simulations of CSR ...
## 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.
##
## Done.
## Suspicious Incident
## Warning: Some expected counts are small; chi^2 approximation may be
## inaccurate
## Q: 21.76
## P-value: 0.8128
## **************
## Generating 99 simulations of CSR ...
## 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.
##
## Done.
## Threat Report
## Warning: Some expected counts are small; chi^2 approximation may be
## inaccurate
## Q: 20.29
## P-value: 0.6392
## **************
## Generating 99 simulations of CSR ...
## 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.
##
## Done.
# Some warnings occur. This is because there are few of some of the point
# events, but a large number of quadrats.
####################################
# Question 2: How are events distributed amongst other events
deaths_TypeOnly <- deaths_AOI[, c(4)] #Create a simplified marks object
deaths_simpMarks_ppp <- ppp(x = deaths_df$Easting, y = deaths_df$Northing, window = convexHull,
marks = deaths_TypeOnly)
## Warning: data contain duplicated points
# alltypes(deaths_simpMarks_ppp, 'Kcross', envelope = T)#Gives lots of
# errors
# Do a lot of cross-K of significant non-CSR events: Including criminal
# event, enemy action, explosive hazard, friendly action
plot(envelope(deaths_simpMarks_ppp, Kcross, i = "Criminal Event", j = "Enemy Action",
nsim = 99, correction = "border"))
## Generating 99 simulations of CSR ...
## 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.
##
## Done.
## lty col key label
## obs 1 1 obs K[list(Criminal.Event, Enemy.Action)][obs](r)
## theo 2 2 theo K[list(Criminal.Event, Enemy.Action)][theo](r)
## hi 1 8 hi K[list(Criminal.Event, Enemy.Action)][hi](r)
## lo 1 8 lo K[list(Criminal.Event, Enemy.Action)][lo](r)
## meaning
## obs observed value of Kcross["Criminal.Event", "Enemy.Action"](r) for data pattern
## theo theoretical value of Kcross["Criminal.Event", "Enemy.Action"](r) for CSR
## hi upper pointwise envelope of Kcross["Criminal.Event", "Enemy.Action"](r) from simulations
## lo lower pointwise envelope of Kcross["Criminal.Event", "Enemy.Action"](r) from simulations
plot(envelope(deaths_simpMarks_ppp, Kcross, i = "Criminal Event", j = "Explosive Hazard",
nsim = 99, correction = "border"))
## Generating 99 simulations of CSR ...
## 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.
##
## Done.
## lty col key label
## obs 1 1 obs K[list(Criminal.Event, Explosive.Hazard)][obs](r)
## theo 2 2 theo K[list(Criminal.Event, Explosive.Hazard)][theo](r)
## hi 1 8 hi K[list(Criminal.Event, Explosive.Hazard)][hi](r)
## lo 1 8 lo K[list(Criminal.Event, Explosive.Hazard)][lo](r)
## meaning
## obs observed value of Kcross["Criminal.Event", "Explosive.Hazard"](r) for data pattern
## theo theoretical value of Kcross["Criminal.Event", "Explosive.Hazard"](r) for CSR
## hi upper pointwise envelope of Kcross["Criminal.Event", "Explosive.Hazard"](r) from simulations
## lo lower pointwise envelope of Kcross["Criminal.Event", "Explosive.Hazard"](r) from simulations
plot(envelope(deaths_simpMarks_ppp, Kcross, i = "Criminal Event", j = "Friendly Action",
nsim = 99, correction = "border"))
## Generating 99 simulations of CSR ...
## 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.
##
## Done.
## lty col key label
## obs 1 1 obs K[list(Criminal.Event, Friendly.Action)][obs](r)
## theo 2 2 theo K[list(Criminal.Event, Friendly.Action)][theo](r)
## hi 1 8 hi K[list(Criminal.Event, Friendly.Action)][hi](r)
## lo 1 8 lo K[list(Criminal.Event, Friendly.Action)][lo](r)
## meaning
## obs observed value of Kcross["Criminal.Event", "Friendly.Action"](r) for data pattern
## theo theoretical value of Kcross["Criminal.Event", "Friendly.Action"](r) for CSR
## hi upper pointwise envelope of Kcross["Criminal.Event", "Friendly.Action"](r) from simulations
## lo lower pointwise envelope of Kcross["Criminal.Event", "Friendly.Action"](r) from simulations
plot(envelope(deaths_simpMarks_ppp, Kcross, i = "Enemy Action", j = "Criminal Event",
nsim = 99, correction = "border"))
## Generating 99 simulations of CSR ...
## 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.
##
## Done.
## lty col key label
## obs 1 1 obs K[list(Enemy.Action, Criminal.Event)][obs](r)
## theo 2 2 theo K[list(Enemy.Action, Criminal.Event)][theo](r)
## hi 1 8 hi K[list(Enemy.Action, Criminal.Event)][hi](r)
## lo 1 8 lo K[list(Enemy.Action, Criminal.Event)][lo](r)
## meaning
## obs observed value of Kcross["Enemy.Action", "Criminal.Event"](r) for data pattern
## theo theoretical value of Kcross["Enemy.Action", "Criminal.Event"](r) for CSR
## hi upper pointwise envelope of Kcross["Enemy.Action", "Criminal.Event"](r) from simulations
## lo lower pointwise envelope of Kcross["Enemy.Action", "Criminal.Event"](r) from simulations
plot(envelope(deaths_simpMarks_ppp, Kcross, i = "Enemy Action", j = "Explosive Hazard",
nsim = 99, correction = "border"))
## Generating 99 simulations of CSR ...
## 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.
##
## Done.
## lty col key label
## obs 1 1 obs K[list(Enemy.Action, Explosive.Hazard)][obs](r)
## theo 2 2 theo K[list(Enemy.Action, Explosive.Hazard)][theo](r)
## hi 1 8 hi K[list(Enemy.Action, Explosive.Hazard)][hi](r)
## lo 1 8 lo K[list(Enemy.Action, Explosive.Hazard)][lo](r)
## meaning
## obs observed value of Kcross["Enemy.Action", "Explosive.Hazard"](r) for data pattern
## theo theoretical value of Kcross["Enemy.Action", "Explosive.Hazard"](r) for CSR
## hi upper pointwise envelope of Kcross["Enemy.Action", "Explosive.Hazard"](r) from simulations
## lo lower pointwise envelope of Kcross["Enemy.Action", "Explosive.Hazard"](r) from simulations
plot(envelope(deaths_simpMarks_ppp, Kcross, i = "Enemy Action", j = "Friendly Action",
nsim = 99, correction = "border"))
## Generating 99 simulations of CSR ...
## 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.
##
## Done.
## lty col key label
## obs 1 1 obs K[list(Enemy.Action, Friendly.Action)][obs](r)
## theo 2 2 theo K[list(Enemy.Action, Friendly.Action)][theo](r)
## hi 1 8 hi K[list(Enemy.Action, Friendly.Action)][hi](r)
## lo 1 8 lo K[list(Enemy.Action, Friendly.Action)][lo](r)
## meaning
## obs observed value of Kcross["Enemy.Action", "Friendly.Action"](r) for data pattern
## theo theoretical value of Kcross["Enemy.Action", "Friendly.Action"](r) for CSR
## hi upper pointwise envelope of Kcross["Enemy.Action", "Friendly.Action"](r) from simulations
## lo lower pointwise envelope of Kcross["Enemy.Action", "Friendly.Action"](r) from simulations
plot(envelope(deaths_simpMarks_ppp, Kcross, i = "Explosive Hazard", j = "Criminal Event",
nsim = 99, correction = "border"))
## Generating 99 simulations of CSR ...
## 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.
##
## Done.
## lty col key label
## obs 1 1 obs K[list(Explosive.Hazard, Criminal.Event)][obs](r)
## theo 2 2 theo K[list(Explosive.Hazard, Criminal.Event)][theo](r)
## hi 1 8 hi K[list(Explosive.Hazard, Criminal.Event)][hi](r)
## lo 1 8 lo K[list(Explosive.Hazard, Criminal.Event)][lo](r)
## meaning
## obs observed value of Kcross["Explosive.Hazard", "Criminal.Event"](r) for data pattern
## theo theoretical value of Kcross["Explosive.Hazard", "Criminal.Event"](r) for CSR
## hi upper pointwise envelope of Kcross["Explosive.Hazard", "Criminal.Event"](r) from simulations
## lo lower pointwise envelope of Kcross["Explosive.Hazard", "Criminal.Event"](r) from simulations
plot(envelope(deaths_simpMarks_ppp, Kcross, i = "Explosive Hazard", j = "Enemy Action",
nsim = 99, correction = "border"))
## Generating 99 simulations of CSR ...
## 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.
##
## Done.
## lty col key label
## obs 1 1 obs K[list(Explosive.Hazard, Enemy.Action)][obs](r)
## theo 2 2 theo K[list(Explosive.Hazard, Enemy.Action)][theo](r)
## hi 1 8 hi K[list(Explosive.Hazard, Enemy.Action)][hi](r)
## lo 1 8 lo K[list(Explosive.Hazard, Enemy.Action)][lo](r)
## meaning
## obs observed value of Kcross["Explosive.Hazard", "Enemy.Action"](r) for data pattern
## theo theoretical value of Kcross["Explosive.Hazard", "Enemy.Action"](r) for CSR
## hi upper pointwise envelope of Kcross["Explosive.Hazard", "Enemy.Action"](r) from simulations
## lo lower pointwise envelope of Kcross["Explosive.Hazard", "Enemy.Action"](r) from simulations
plot(envelope(deaths_simpMarks_ppp, Kcross, i = "Explosive Hazard", j = "Friendly Action",
nsim = 99, correction = "border"))
## Generating 99 simulations of CSR ...
## 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.
##
## Done.
## lty col key label
## obs 1 1 obs K[list(Explosive.Hazard, Friendly.Action)][obs](r)
## theo 2 2 theo K[list(Explosive.Hazard, Friendly.Action)][theo](r)
## hi 1 8 hi K[list(Explosive.Hazard, Friendly.Action)][hi](r)
## lo 1 8 lo K[list(Explosive.Hazard, Friendly.Action)][lo](r)
## meaning
## obs observed value of Kcross["Explosive.Hazard", "Friendly.Action"](r) for data pattern
## theo theoretical value of Kcross["Explosive.Hazard", "Friendly.Action"](r) for CSR
## hi upper pointwise envelope of Kcross["Explosive.Hazard", "Friendly.Action"](r) from simulations
## lo lower pointwise envelope of Kcross["Explosive.Hazard", "Friendly.Action"](r) from simulations
plot(envelope(deaths_simpMarks_ppp, Kcross, i = "Friendly Action", j = "Criminal Event",
nsim = 99, correction = "border"))
## Generating 99 simulations of CSR ...
## 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.
##
## Done.
## lty col key label
## obs 1 1 obs K[list(Friendly.Action, Criminal.Event)][obs](r)
## theo 2 2 theo K[list(Friendly.Action, Criminal.Event)][theo](r)
## hi 1 8 hi K[list(Friendly.Action, Criminal.Event)][hi](r)
## lo 1 8 lo K[list(Friendly.Action, Criminal.Event)][lo](r)
## meaning
## obs observed value of Kcross["Friendly.Action", "Criminal.Event"](r) for data pattern
## theo theoretical value of Kcross["Friendly.Action", "Criminal.Event"](r) for CSR
## hi upper pointwise envelope of Kcross["Friendly.Action", "Criminal.Event"](r) from simulations
## lo lower pointwise envelope of Kcross["Friendly.Action", "Criminal.Event"](r) from simulations
plot(envelope(deaths_simpMarks_ppp, Kcross, i = "Friendly Action", j = "Enemy Action",
nsim = 99, correction = "border"))
## Generating 99 simulations of CSR ...
## 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.
##
## Done.
## lty col key label
## obs 1 1 obs K[list(Friendly.Action, Enemy.Action)][obs](r)
## theo 2 2 theo K[list(Friendly.Action, Enemy.Action)][theo](r)
## hi 1 8 hi K[list(Friendly.Action, Enemy.Action)][hi](r)
## lo 1 8 lo K[list(Friendly.Action, Enemy.Action)][lo](r)
## meaning
## obs observed value of Kcross["Friendly.Action", "Enemy.Action"](r) for data pattern
## theo theoretical value of Kcross["Friendly.Action", "Enemy.Action"](r) for CSR
## hi upper pointwise envelope of Kcross["Friendly.Action", "Enemy.Action"](r) from simulations
## lo lower pointwise envelope of Kcross["Friendly.Action", "Enemy.Action"](r) from simulations
plot(envelope(deaths_simpMarks_ppp, Kcross, i = "Friendly Action", j = "Explosive Hazard",
nsim = 99, correction = "border"))
## Generating 99 simulations of CSR ...
## 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.
##
## Done.
## lty col key label
## obs 1 1 obs K[list(Friendly.Action, Explosive.Hazard)][obs](r)
## theo 2 2 theo K[list(Friendly.Action, Explosive.Hazard)][theo](r)
## hi 1 8 hi K[list(Friendly.Action, Explosive.Hazard)][hi](r)
## lo 1 8 lo K[list(Friendly.Action, Explosive.Hazard)][lo](r)
## meaning
## obs observed value of Kcross["Friendly.Action", "Explosive.Hazard"](r) for data pattern
## theo theoretical value of Kcross["Friendly.Action", "Explosive.Hazard"](r) for CSR
## hi upper pointwise envelope of Kcross["Friendly.Action", "Explosive.Hazard"](r) from simulations
## lo lower pointwise envelope of Kcross["Friendly.Action", "Explosive.Hazard"](r) from simulations
####################################
# Question 3 - Predict the number of events of one type based upon the
# density surface of another
# Explore parameters for predicting number of Criminal Events using
# Difference K-Function: Criminal Event and Enemy Action
KCrimEvent <- Kest(deaths_ppp[deaths_marks$Type == "Criminal Event", ], correction = "best")
KEnAct <- Kest(deaths_ppp[deaths_marks$Type == "Enemy Action", ], correction = "best")
KdiffCEEA <- KCrimEvent$iso - KEnAct$iso # Difference in K Functions
plot(KdiffCEEA, type = "o", ylim = c(-1e+07, 1e+07))
abline(h = 0)
median(KdiffCEEA)/median(KCrimEvent$iso)
## [1] 0.2827
# Criminal Event and Explosive Hazard
KExploHaz <- Kest(deaths_ppp[deaths_marks$Type == "Explosive Hazard", ], correction = "best")
KdiffCEEH <- KCrimEvent$iso - KExploHaz$iso
plot(KdiffCEEH, type = "o", ylim = c(-1e+07, 1e+07))
abline(h = 0)
median(KdiffCEEH)/median(KCrimEvent$iso)
## [1] 0.4569
# Criminal Event and Friendly Action
KFriendAct <- Kest(deaths_ppp[deaths_marks$Type == "Friendly Action", ], correction = "best")
KdiffCEFA <- KCrimEvent$iso - KFriendAct$iso
plot(KdiffCEFA, type = "o", ylim = c(-1e+07, 1e+07))
abline(h = 0)
median(KdiffCEFA)/median(KCrimEvent$iso)
## [1] 0.152
# Enemy Action and Explosive Hazard
KdiffEAEH <- KEnAct$iso - KExploHaz$iso
plot(KdiffEAEH, type = "o", ylim = c(-1e+07, 1e+07))
abline(h = 0)
median(KdiffEAEH)/median(KEnAct$iso)
## [1] 0.2433
# Enemy Action and Friendly Action
KdiffEAFA <- KEnAct$iso - KFriendAct$iso
plot(KdiffEAFA, type = "o", ylim = c(-1e+07, 1e+07))
abline(h = 0)
median(KdiffEAFA)/median(KEnAct$iso)
## [1] -0.1993
# Explosive Hazard and Friendly Action
KdiffEHFA <- KExploHaz$iso - KFriendAct$iso
plot(KdiffEHFA, type = "o", ylim = c(-1e+07, 1e+07))
abline(h = 0)
median(KdiffEHFA)/median(KExploHaz$iso)
## [1] -0.5872
# Develop a point process model for Criminal Event using density surfaces
# from other types of events Two Approaches: From difference in K
# functions,it appears as though Friendly Action might best predict the
# number of Criminal Events From the Cross K functions it appears , Enemy
# Action might best predict Criminal Events
# Create ppp objects and/or density surfaces to be used in a ppm to
# predict lambda intensity
# Criminal Event - This will be our Independent Variable (no density
# surface needed)
deaths_CrimEvent <- deaths_df[deaths_df$Type == "Criminal Event", ]
deaths_ppp_CrimEvent <- ppp(x = deaths_CrimEvent$Easting, y = deaths_CrimEvent$Northing,
window = convexHull)
## Warning: data contain duplicated points
# Enemy Action
deaths_EnAct <- deaths_df[deaths_df$Type == "Enemy Action", ]
deaths_ppp_EnAct <- ppp(x = deaths_EnAct$Easting, y = deaths_EnAct$Northing,
window = convexHull)
## Warning: data contain duplicated points
density_EA <- density(deaths_ppp_EnAct)
# Friendly Action
deaths_FriendAct <- deaths_df[deaths_df$Type == "Friendly Action", ]
deaths_ppp_FriendAct <- ppp(x = deaths_FriendAct$Easting, y = deaths_FriendAct$Northing,
window = convexHull)
## Warning: data contain duplicated points
density_FA <- density(deaths_ppp_FriendAct)
# Use point process model to simulate Criminal Event points based on a
# raster image of density surfaces of other event types Criminal Event ~
# Enemy Action
ppm.CEEA <- ppm(deaths_ppp_CrimEvent, ~density_EA, covariates = list(density_EA = density_EA))
## Warning: Some tiles with zero area contain points
## Warning: Some weights are zero
## Warning: Values of the covariate 'density_EA' were NA or undefined at
## 0.01% (1 out of 18395) of the quadrature points. Occurred while executing:
## ppm(deaths_ppp_CrimEvent, ~density_EA, covariates = list(density_EA =
## density_EA))
ptsim_CEEA1 <- simulate(ppm.CEEA) #Simulates CSR using ppm model
ptsim_CEEA2 <- simulate(ppm.CEEA)
plot(ptsim_CEEA1, main = "1st Simulation") #Plot the simulated CSR with a variable lambda intensity
plot(ptsim_CEEA2, main = "2nd Simulation")
plot(predict.ppm(ppm.CEEA, type = "lambda"), main = "Predicted Criminal Event Intensity\nbased on Enemy Action") #predicted intensity
ppm.CEEA.k <- envelope(ppm.CEEA, Kest, nsim = 99, correction = "border", verbose = F)
plot(ppm.CEEA.k, main = "Observed Criminal Events against\nModel based on Enemy Action")
## lty col key label
## obs 1 1 obs K[obs](r)
## mmean 2 2 mmean bar(K)(r)
## hi 1 8 hi K[hi](r)
## lo 1 8 lo K[lo](r)
## meaning
## obs observed value of K(r) for data pattern
## mmean sample mean of K(r) from simulations
## hi upper pointwise envelope of K(r) from simulations
## lo lower pointwise envelope of K(r) from simulations
# If the D.V. density plot is good at predicting criminal events, we would
# expect to now see a CSR process Because the line is above the CSR
# envelope, we cannot determine that Enemy Action is a good predictor of
# Criminal Events
# Criminal Event ~ Friendly Action
ppm.CEFA <- ppm(deaths_ppp_CrimEvent, ~density_FA, covariates = list(density_FA = density_FA))
## Warning: Some tiles with zero area contain points
## Warning: Some weights are zero
## Warning: Values of the covariate 'density_FA' were NA or undefined at
## 0.01% (1 out of 18395) of the quadrature points. Occurred while executing:
## ppm(deaths_ppp_CrimEvent, ~density_FA, covariates = list(density_FA =
## density_FA))
ptsim_CEFA1 <- simulate(ppm.CEFA)
ptsim_CEFA2 <- simulate(ppm.CEFA)
plot(ptsim_CEFA1, main = "1st Simulation")
plot(ptsim_CEFA2, main = "2nd Simulation")
plot(predict.ppm(ppm.CEFA, type = "lambda"), main = "Predicted Criminal Event Intensity\nbased on Friendly Action") #predicted intensity
ppm.CEFA.k <- envelope(ppm.CEFA, Kest, nsim = 99, correction = "border", verbose = F)
plot(ppm.CEFA.k, main = "Observed Criminal Events against\nModel based on Friendly Action")
## lty col key label
## obs 1 1 obs K[obs](r)
## mmean 2 2 mmean bar(K)(r)
## hi 1 8 hi K[hi](r)
## lo 1 8 lo K[lo](r)
## meaning
## obs observed value of K(r) for data pattern
## mmean sample mean of K(r) from simulations
## hi upper pointwise envelope of K(r) from simulations
## lo lower pointwise envelope of K(r) from simulations
# If the D.V. density plot is good at predicting criminal events, we would
# expect to now see a CSR process Because the line is above the CSR
# envelope, we cannot determine that Friendly Action is a good predictor
# of Criminal Events
# Sadly, we were not able to show that it was possible to predict criminal
# events at this time.
~Fin