This exercise involves the exploratory analysis of point patterns. The goal is to learn how to visualize and do some preliminary analysis on sets of locations (points). Many of the skills learned here will be essential for hypothesis tests involving spatial point patterns. We will use a number of libraries today including:
library(spatstat)
library(sp)
library(rgdal)
library(maptools)
library(ggmap)
Read the csv file “deaths only.csv” into R using read.csv():
id <- read.csv("Deaths only.csv")
summary(id)
The data is a little sloppy, lets recode the types so that “CRIMINAL EVENT” and “Criminal Event” are not separate categories. To to do this we'll use the function sub() which substitutes one text string for another text string in the specified object:
id$Type <- sub(id$Type, pattern = "CRIMINAL EVENT", replacement = "Criminal Event")
id$Type <- sub(id$Type, pattern = "criminal event", replacement = "Criminal Event")
id$Type <- sub(id$Type, pattern = "EXPLOSIVE HAZARD", replacement = "Explosive Hazard")
In order to be sure the old categories are removed from the list of possible event types we have to run:
id$Type <- factor(id$Type)
levels(id$Type)
## [1] "Criminal Event" "Enemy Action" "Explosive Hazard"
## [4] "Friendly Action" "Friendly Fire" "Non-Combat Event"
## [7] "Other" "Suspicious Incident" "Threat Report"
The spatstat library is for the analysis and manipulation of spatial points. In a spatstat point file the attributes of each location are stored in “marks”. Each point on the map can have zero, one, or many associated marks. A file with too many marks is cumbersome, select a subset of attributes to serve as potential marks:
id.marks.simple <- id[, 4]
id.marks.full <- id[, c(4, 7, 9:18)]
If you looked carefully at summary(id) you noticed that latitude and longitude are attributes in the table, we'll use these to map out the events:
id.spp <- SpatialPointsDataFrame(coords = id[, 19:20], proj4string = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"),
data = id.marks.full)
The above text creates a SpatialPointsDataFrame which is a type of map that stores point objects, if you're a GIS person its roughly equivalent to a point shapefile. I'm guessing at the projection of the original data.
To make maps we'll use a very new package called ggmaps that allows us to send a query to online mapping services like Google maps. It returns a georeferenced image of the specified type and scale. We can then plot on top of that image. The code below uses a package called ggplot2 which we'll look at more carefully later in the semester. The first step is to get the base map.
map <- get_map(location = "Iraq", maptype = "hybrid", zoom = 6) #get map
ggmap(map, extent = "panel")
I think its cool how ggmaps was able to translate the word “Iraq” into a map. Now we can overlay the WikiLeaks Iraq war data:
ggmap(map, extent = "panel") + geom_point(data = id, aes(x = id$Longitude, y = id$Latitude),
alpha = 0.1, color = "red") + ggtitle("WikiLeaks Iraq Data \n Incidents involving Death(s)")
Alter the map the so that the size of the dots are linked to the total number of deaths in each events (note the use of size inside of aes()):
ggmap(map, extent = "panel") + geom_point(data = id, aes(x = id$Longitude, y = id$Latitude,
size = Total.deaths), alpha = 0.3, color = "red") + ggtitle("WikiLeaks Iraq Data \n Incidents involving Death(s)")
Map a subset of events, for example, those involving detentions:
ggmap(map, extent = "panel") + geom_point(data = id[id$Enemy.detained > 0, ],
aes(x = Longitude, y = Latitude, size = Enemy.detained), alpha = 0.5, color = "green") +
ggtitle("WikiLeaks Iraq Data \n Incidents involving Dententions(s)")
Again, we'll discuss the use of ggplot2 more later in the semester. Quickly, the code above plots 1 layer which is the Google map and then overlays a set of points, the data for the points is the data.frame id. We then map specific columns in id to visual variables, such as size and position using aes(). I'm a big fan of ggplot2's functionality for visualization (even though it kinda sucks for mapping).
Changing gears a bit we'll begin using the package spatstat do some analysis. Spatstat is a R library designed in Australia by Adrian Baddeley. The logic of spatstat is a little odd from a geographic perspective, but once you understand how it works it is a really useful and robust library.
Many of the analytic tools in spatstat are based on distance. Working with decimal degrees is a pain, they're difficult to interpret as distance units because 1 degree is a different distance on different parts of the globe. Re-project the data to a coordinate system that uses meters instead of degrees of latitude/longitude. Knowing how many deaths occurred per decimal degree isn't that useful. My knowledge of Iraqi coordinate systems is slim, I've found one that seems to work, but I'm not 100% sure it's correct…
## Convert to a different coordinate system. So that coordinates are in
## meters not degrees
id.spp.rp <- spTransform(id.spp, CRS("+init=epsg:3839"))
Now, we'll save the new coordinates:
id.spp.rp$Easting <- coordinates(id.spp.rp)[, 1]
id.spp.rp$Northing <- coordinates(id.spp.rp)[, 2]
I find it difficult to look at the entire country, because violence rates vary so much. I'm going to limit my analysis to Baghdad:
## Remove all events outside of Baghdad and convert back to a table
id.df.rp.bag <- id.spp.rp@data[id.spp.rp@data$Region == "MND-BAGHDAD", ]
Notice my use of @data I have pulled out a subset of the map and saved it as a non-spatial data.frame. The spatial information is preserved in the columns Easting and Northing. We can still make a “map” of the points:
plot(x = id.df.rp.bag$Easting, y = id.df.rp.bag$Northing)
There is an eastern outlier and few northern outliers, let's remove them:
id.df.rp.bag <- id.df.rp.bag[id.df.rp.bag$Northing < 4960000 & id.df.rp.bag$Easting >
9960000, ]
plot(x = id.df.rp.bag$Easting, y = id.df.rp.bag$Northing)
On much of the map the points are still very sparse, let's zoom into the region with the dense concentration of events:
id.df.rp.bag <- id.df.rp.bag[(id.df.rp.bag$Northing < 4950000 & id.df.rp.bag$Northing >
4940000) & (id.df.rp.bag$Easting > 1e+07 & id.df.rp.bag$Easting < 10010000),
]
These are the points we'll use in our analysis.
spatstatToday we will use three types of spatstat() objects:
ppp objects must have an owin.ppp object.The first step in making a ppp for spatial analysis is to define the window. Define the window using convexhull.xy() which traces the outermost points of a scatter plot. Give the function the x and y coordinates of the points, and it returns an owin object – which is necessary for making a point pattern object…
id.chull <- convexhull.xy(x = id.df.rp.bag$Easting, y = id.df.rp.bag$Northing)
It is possible (and sometimes necessary) to manually draw a window. This can be done with the following code. It is important to draw the window counter clockwise:
plot(x = id.df.rp.bag$Easting, y = id.df.rp.bag$Northing)
bdry <- locator()
id.win <- owin(poly = bdry)
Use this window to create the ppp object (ignore the warning):
id.ppp <- ppp(x = id.df.rp.bag$Easting, y = id.df.rp.bag$Northing, window = id.chull,
marks = id.df.rp.bag[, c(1, 3:12)])
## Warning: data contain duplicated points
unitname(id.ppp) <- c("meter", "meters") ##Assign names to units
If you want to change the window to one that you've drawn manually:
# To change the window of an existing ppp object
id.ppp <- id.ppp[id.win]
Now visualize the ppp object:
p1 <- plot.ppp(id.ppp)
legend(max(coords(id.ppp))[1] + 1000, mean(coords(id.ppp))[2], pch = p1, legend = names(p1),
cex = 0.5)
## Warning: mean(<data.frame>) is deprecated. Use colMeans() or sapply(*,
## mean) instead.
Notice that by default the first “mark” is use to label the points. This mark was the “Type” column in the original data, it describes the type of incident. The data are so dense that the resulting graph is impossible to read. There are alternative visualizations:
## Plot the density of points using a gaussian (normal) kernel with a
## 1000m standard deviation
plot(density(id.ppp, sigma = 1000))
This density plot places a kernel over the center of each pixel and counts up the number of points per unit area. The points further away from the center of each pixel count less toward the density of the pixel than points that are close to the center. This is the standard kernel density estimator that is used in most GIS programs.
The density command creates a pixel image or im object. These objects allow us to do cool stuff when combined with points, as we will see later…
We can make a contour map of the density:
contour(density(id.ppp, 1000), axes = F)
## NULL
A 3D, “perspective” plot, of the data:
layout(matrix(c(1, 1, 1, 2, 3, 4), 2, 3, byrow = TRUE))
persp(density(id.ppp, 1000))
## [,1] [,2] [,3] [,4]
## [1,] 2.018e-04 0.000e+00 0.000e+00 0.000e+00
## [2,] 0.000e+00 5.223e-05 -1.949e-04 1.949e-04
## [3,] 0.000e+00 1.834e+04 4.915e+03 -4.915e+03
## [4,] -2.019e+03 -2.590e+02 9.611e+02 -9.601e+02
persp(density(id.ppp, 1000), theta = 150, phi = 15) #rotate by adjusting theta
## [,1] [,2] [,3] [,4]
## [1,] -1.748e-04 -2.612e-05 9.747e-05 -9.747e-05
## [2,] 1.009e-04 -4.524e-05 1.688e-04 -1.688e-04
## [3,] -5.814e-13 1.834e+04 4.915e+03 -4.915e+03
## [4,] 1.250e+03 4.843e+02 -1.813e+03 1.814e+03
persp(density(id.ppp, 1000), theta = 120, phi = 15) #rotate by adjusting theta
## [,1] [,2] [,3] [,4]
## [1,] -1.009e-04 -4.524e-05 1.688e-04 -1.688e-04
## [2,] 1.748e-04 -2.612e-05 9.747e-05 -9.747e-05
## [3,] -1.007e-12 1.834e+04 4.915e+03 -4.915e+03
## [4,] 1.453e+02 5.810e+02 -2.174e+03 2.175e+03
persp(density(id.ppp, 1000), theta = 100, phi = 15) #rotate by adjusting theta
## [,1] [,2] [,3] [,4]
## [1,] -3.505e-05 -5.144e-05 1.920e-04 -1.920e-04
## [2,] 1.988e-04 -9.070e-06 3.385e-05 -3.385e-05
## [3,] -1.145e-12 1.834e+04 4.915e+03 -4.915e+03
## [4,] -6.322e+02 5.588e+02 -2.091e+03 2.092e+03
dev.off() #clearlayout
## null device
## 1
Alternatively, we can visualize each type of event:
plot(split(id.ppp))
plot(density(split(id.ppp)))
contour(density(split(id.ppp)), axes = F)
Notice that by default only the “Type” mark is plotted. Map different marks as follows:
plot(id.ppp, which.marks = "Total.deaths", maxsize = 3000, main = "Total Deaths")
## 0 50 100 150 200
## 0.0 828.7 1657.5 2486.2 3314.9
It takes some fiddling with maxsize to get a decent plot, by default all the circles are tiny.
There are a lot of potential ways to visualize points in R, we've just scratched the surface. One really cool thing about spatstat is that its possible to use pixel images (like the density maps) in the analysis of points. Let's make a weighted density map that reflects the total number of deaths at each location in our window. We'll use this map to do some preliminary analysis of the points. To start with lets simplify the ppp object:
id.ppp.simp <- unmark(id.ppp) #remove all data from the points
deaths <- marks(id.ppp)["Total.deaths"] #create a list of # of deaths
id.ppp.simp.death <- id.ppp.simp %mark% deaths #make a point pattern (ppp) with a single variable
The density function has a lot of options, for example we can weight each point by the total deaths using the “weight” option, we can specify the type of kernel (using kernel), and the size of the pixels (using eps):
id.death.im <- density(id.ppp.simp.death, sigma = 1000, weights = marks(id.ppp.simp.death),
kernel = "rectangular", eps = 100) #eps sets pixel size in map units
plot(id.death.im, main = "Weighted Density Map of \nWikiLeaks Deaths in Baghdad")
I don't really like the default color scheme for the density maps but changing the colors is a pain. To change the colors we first have to define a function called a color ramp. The color ramp will return a specified number of colors between the two user selected colors. To see a list of possible colors type colors().
# We can change the colors by making a custom color ramp
grRd <- colorRampPalette(c("grey60", "red"))
The results of the above code is a function that returns the specified number of colors on a color ramp between grey and red. For example, grRd(7) returns seven colors with grey on the low end and red on the high end. To use this function we have to make a colourmap that assigns each of the seven colors to a range of values.
co <- colourmap(grRd(5), range = c(3.37e-06, 0.00036)) #range from summary(id.death.im)
plot(id.death.im, main = "Weighted Density Map of \nWikiLeaks Deaths in Baghdad",
col = co)
It is possible to do map algebra like operations with pixel images (like the density map). Its also possible to export the images to a table, such that the each pixel becomes a cell in the table.
The death density image can be used to divide the study region into areas of high, medium, and low “death density”. The function tess is used to create a special type of spatstat object called a tessellation. Tessellations can be used to group points.
death.map.breaks <- quantile(id.death.im, probs = (0:3)/3) #define density cut points
death.map.cuts <- cut(id.death.im, breaks = death.map.breaks, labels = c("Low DD",
"Medium DD", "High DD")) #cut the image
death.map.quartile <- tess(image = death.map.cuts)
plot(death.map.quartile)
We can also use the tessellation to count the number of events in each category on the map:
qCount.im <- quadratcount(X = id.ppp, tess = death.map.quartile)
## Warning: Tessellation does not contain all the points of X
plot(qCount.im)
We could also divide the window up into a \( n \) by \( n \) grid and count the number of events in each grid cell:
qCount <- quadratcount(X = id.ppp, nx = 4, ny = 4) #returns the count of points in each quadrant
plot(qCount)
We can also use the density map to look at how the intensity of events (number of points per unit area) relates to the density of death image:
plot(rhohat(unmark(id.ppp), id.death.im))
We'll discuss rhohat more in lecture. This plot basically shows you the density of points as a function of the values on the weighted death density map. The yellow (High DD) regions on the map contain not only more points but the points are higher density than parts of the death density of image that have lower numbers of deaths.
We could also split the data frame up into a series of point patterns, one for each type of event:
id.ppp.types <- split(id.ppp, which.marks = "Type")
id.ppp.types <- unmark(id.ppp.types)
The object id.ppp.types now holds a bunch of unmarked point patterns, each type of event is in its own ppp object. For example, we could look at only the criminal events:
id.ppp.types$"Criminal Event"
There are so many criminal events that its difficult to make heads or tails of the map. Let's try to filter the significant clusters (high density areas) from the background rate (noise). We'll do this using an approach developed by Byers and Raferty (1998) called Nearest Neighbor Cleaning. We'll discuss the theory in class, the basic idea is to estimate the parameters of two point generating processes (one for noise and one for signal). Once the parameters of each process are estimated points are assigned to one of the two processes.
nnc <- nnclean(X = id.ppp.types$"Criminal Event", k = 20)
## Iteration 1 logLik = -361421.031345153 p = 0.7109
## Iteration 2 logLik = -358561.899117613 p = 0.6307
## Iteration 3 logLik = -356876.844782523 p = 0.5624
## Iteration 4 logLik = -356003.613549312 p = 0.4991
## Iteration 5 logLik = -355766.598100418 p = 0.4594
## Estimated parameters:
## p [cluster] = 0.45944
## lambda [cluster] = 7.1591e-05
## lambda [noise] = 1.1937e-05
plot(split(nnc))
The map above shows intense clusters of criminal activity against the background “noise.” In the next exercise we'll do hypothesis tests to understand the processes that generated the observed point patterns.