X Wang - Lab 4 Appendix

GEOG 5023: Quantitative Methods in Geography

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 complete spatial randomness (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?

Question 0 - Load libraries, refine AOI, and prepare data for PPA.

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.31-1 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)

## Warning: package 'sp' was built under R version 2.15.3

library(rgdal)

## Warning: package 'rgdal' was built under R version 2.15.3

## rgdal: version: 0.8-6, (SVN revision Unversioned directory) Geospatial
## Data Abstraction Library extensions to R successfully loaded Loaded GDAL
## runtime: GDAL 1.9.0, released 2011/12/29 Path to GDAL shared files:
## /Library/Frameworks/R.framework/Versions/2.15/Resources/library/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480] Path to
## PROJ.4 shared files:
## /Library/Frameworks/R.framework/Versions/2.15/Resources/library/rgdal/proj

library(maptools)

## Loading required package: foreign

## Loading required package: grid

## Loading required package: lattice

## Checking rgeos availability: FALSE Note: when rgeos is not available,
## polygon geometry computations in maptools depend on gpclib, which has a
## restricted licence. It is disabled by default; to enable gpclib, type
## gpclibPermit()

library(ggmap)

## Loading required package: ggplot2

## Load data
deaths <- read.csv("/Users/xiwang/Dropbox/GEOG 5023 - offline/Labs/Lab 4/Deaths only.csv")
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"

str(deaths$Type)  # Check structure to ensure conversion back to factor type

##  Factor w/ 9 levels "Criminal Event",..: 6 1 3 2 4 3 3 3 3 3 ...

# Look at map of Iraq and overlay WikiLeaks Iraq War data
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)")

# Explore and subset the data to define 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 = 0.12 * Total.deaths), alpha = 0.3, color = "red") + 
    ggtitle("Incidents involving Death(s)\nin Baghdad near AOI")  # Note that this map doesn't show the exact AOI we are using.  See the attached map on the lab writeup instead

## Warning: Removed 35539 rows containing missing values (geom_point).

# Subset the data
deaths_AOI <- deaths[(deaths$Longitude < 44.465 & deaths$Longitude > 44.33) & 
    (deaths$Latitude > 33.235 & deaths$Latitude < 33.33), ]

AOI coordinates specify southern Baghdad around the Tigris River centered on the Two Stories and Al Jadriyah Bridges. The lower panel shows the AOI.

deaths_marks <- deaths_AOI[, c(4, 7, 8, 18)]
names(deaths_marks)

## [1] "Type"         "Region"       "Attack.on"    "Total.deaths"

str(deaths_marks)

## 'data.frame':    5928 obs. of  4 variables:
##  $ Type        : Factor w/ 9 levels "Criminal Event",..: 2 1 1 2 2 1 1 2 3 1 ...
##  $ Region      : Factor w/ 8 levels "","MND-BAGHDAD",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ Attack.on   : Factor w/ 3 levels "ENEMY","FRIEND",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Total.deaths: int  1 1 1 1 1 1 2 2 1 2 ...

deaths_marks$Total.deaths <- factor(deaths_marks$Total.deaths)  # Make sure all elements are factor type

# 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 ratio at 1.0 (no distortion)

Question 1 - For the 9 event types, describe the pattern in relation to CSR.

# Start making a point pattern (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

# Point and density maps
plot(deaths_ppp)  # Point map

##      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))  # Kernel density map

plot(split(deaths_ppp))

plot(density(split(deaths_ppp), sigma = 500))

# 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 ppp object with a single variable
plot(density(death_simp_ppp, sigma = 500, weights = marks(death_simp_ppp)))

## Error: error in evaluating the argument 'x' in selecting a method for
## function 'plot': Error in Summary.factor(c(1L, 1L, 2L, 2L, 1L), na.rm =
## FALSE) : sum not meaningful for factors Calls: density ... pixellate.ppp
## -> tapply -> lapply -> Summary.factor

# 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 - Are the events distributed in relationship to one another?

deaths_TypeOnly <- deaths_AOI[, c(4)]
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 errors
# because marks have to be converted to a factor string

# 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 - Is it possible to predict the number of events of one type based upon the density surface of another?

Explore parameters for predicting number of Criminal Events.

# The “cross-K” function returns the number of points of type b within r
# units of a type-a point. The difference in K functions examines
# whetheter two types of events (a and b) have similar spatial
# distributions, so the difference should be zero at all scales. Look at
# difference in K functions between all event types that do not apper to
# come from a CSR process: Criminal Event, Enemy Action, Explosive Hazard,
# Friendly Action

# 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)  # 0.2827: The K function for Enemy Action is 28% smaller than the K function for Criminal Event at the median, meaning the former is less clustered than the latter at this scale.

## [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)  # 0.4569: The K function for Explosive Hazard is 46% smaller than the K function for Criminal Event at the median, meaning the former is substantially less clustered than the latter at this scale.

## [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)  # 0.1519: The K function for Friendly Action is 15% smaller than the K function for Criminal Event at the median, meaning the former is a bit less clustered than the latter at that scale.

## [1] 0.1519

# 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)  # 0.2433: The K function for Explosive Hazard is 24% smaller than the K function for Enemy Action at the median, meaning the former is less clustered than the latter at that scale.

## [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)  #-0.1995: The K function for Friendly Action is 20% larger than the K function for Enemy Action at the median, meaning the former is more clustered than the latter at that scale.

## [1] -0.1995

# 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)  #0.5873: The K function for Friendly Action is 59% larger than the K function for Explosive Hazard at the median, meaning the former is substantially more clustered than the latter at that scale.

## [1] -0.5873

Develop a point process model for Criminal Event using density surfaces from other types of events.

# Test two approaches: 1) From the difference in K functions,it appears as
# though Friendly Action might best predict the number of Criminal Events.
# 2) From the 'cross-K'' functions, it appears Enemy Action might best
# predict Criminal Events.

# Subset significant event types as separate ppp objects and/or find
# density surfaces to find predicted lambda density of Criminal Event.

# Criminal Event - this is 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)

# Explosive Hazard
deaths_ExploHaz <- deaths_df[deaths_df$Type == "Explosive Hazard", ]
deaths_ppp_ExploHaz <- ppp(x = deaths_ExploHaz$Easting, y = deaths_ExploHaz$Northing, 
    window = convexHull)

## Warning: data contain duplicated points

density_EH <- density(deaths_ppp_ExploHaz)

# 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
plot(deaths_ppp_CrimEvent, main = "Observed Criminal Events")  # Observed events (for reference)

# 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(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

# The results show significant clustering up to about 1800 m. Then the
# observed k-function enters the grey envelope. This means that at some
# scales, the observed data is a plausible realization of the
# inhomogeneous intensity surface we created, however, at scales less than
# 1800 m, it is too clustered. This simulation can be viewed as a failed
# attempt to explain the observed distribution of Criminal Event.

# Criminal Event ~ Explosive Hazard
ppm.CEEH <- ppm(deaths_ppp_CrimEvent, ~density_EH, covariates = list(density_EH = density_EH))

## Warning: Some tiles with zero area contain points

## Warning: Some weights are zero

## Warning: Values of the covariate 'density_EH' were NA or undefined at
## 0.01% (1 out of 18395) of the quadrature points. Occurred while executing:
## ppm(deaths_ppp_CrimEvent, ~density_EH, covariates = list(density_EH =
## density_EH))

ptsim_CEEH1 <- simulate(ppm.CEEH)
ptsim_CEEH2 <- simulate(ppm.CEEH)
plot(ptsim_CEEH1, main = "1st Simulation")

plot(ptsim_CEEH2, main = "2nd Simulation")

plot(predict.ppm(ppm.CEEH, type = "lambda"), main = "Predicted Criminal Event Intensity\nbased on Explosive Hazard")  #predicted intensity

ppm.CEEH.k <- envelope(ppm.CEEH, Kest, nsim = 99, correction = "border", verbose = F)
plot(ppm.CEEH.k, main = "Observed Criminal Events against\nModel based on Explosive Hazard")

##       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

# At scales less than 2100 m, the observed data is too clustered to be a
# plausible realization of the inhomogeneous intensity surface we created.
# This simulation is a failed attempt to explain the observed distribution
# of Criminal Event.

# 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

# At scales less than 1900 m, the observed data is too clustered to be a
# plausible realization of the inhomogeneous intensity surface we created.
# This simulation is a failed attempt to explain the observed distribution
# of Criminal Event.

Sadly, we were not able to show that it was possible to predict criminal events at this time. :(