Transcript Lecture 25
Spatial Statistics
Point Patterns
Spatial Statistics
• Increasing sophistication of GIS
allows archaeologists to apply a
variety of spatial statistics to their
data
– Predictive Modeling
– Intra-site Spatial Analysis
Predictive Modeling 1
• Goal is to predict where sites will be
located
• Usually involves two samples:
known sites, surveyed areas where
sites have not been found
• For each group, we collect data:
slope, aspect, distance to water,
soil, vegetation zone, etc
Predictive Modeling 2
• Must convert nominal scales to
dichotomies or an interval scale of
some kind
• Analysis by logistic regression or
discriminant functions
• For new areas we compute the
probability that a site will be found
Predictive Modeling 3
• Problems:
– Usually purely inductive
– Goal is management not anthropology
– Independent variables are those
gathered for other reasons
Intra-Site Spatial Analyses
• Nearest Neighbor
– Can use on any point plot data (sites,
houses, artifacts)
– Find distance to nearest neighbor for
each item
– Mean nearest neighbor compared to
expected value (random distribution)
Nearest Neighbor
• If observed mean distance is
significantly less than expected, the
points are clustered
• If the mean distance is significantly
more than expected, the points are
evenly spread
• But problems with borders
Clustered
Random
Regular
12
12
12
10
10
10
8
8
8
6
6
6
4
4
4
2
2
2
0
0
0
0
2
4
6
8
Mean Distance = 0.52
Expected Dist = 0.78
Probability = 0.002 *
R = Mean/Exp = 0.67
10
0
2
4
6
8
Mean Distance = .74
Expected Dist = 0.63
Probability = 0.091
R = Mean/Exp = 1.18
10
0
2
4
6
8
Mean Distance = 1.04
Expected Dist = 0.81
Probability = 0.008 *
R = Mean/Exp = 1.29
10
Point Patterns in R
• Package spatstat
• Create a ppp object (point process)
• Plotting and analytical tools are
extensive
# Fix duplicate coordinates by adding 1-2 mm to each
load("C:/Users/Carlson/Documents/Courses/Anth642/R/Data/BTF3a.RData")
Win3a <- data.frame(x=c(982,982,983,983,984.5,985,985,987,987,986.2,
985,985,984.5,984,983.5,983,982.7,982.5),
y=c(1015.5,1021,1021,1022,1022,1021.3,1018,1018,1017.6,1017,
1017,1016.9,1016.8,1016.6,1016.3,1016,1015.6,1015.5))
# coordinates must be counterclockwise
Win3a <- Win3a[order(as.numeric(rownames(Win3a)), decreasing=TRUE),]
library(spatstat)
BTF3a.p <- ppp(BTF3a$East, BTF3a$North, window=owin(poly=Win3a,
unitname=c("meter", "meters")), marks=BTF3a$Type)
summary(BTF3a.p)
plot(BTF3a.p, main="Bifacial Thinning Flakes", cex=.75, chars=16,
cols=c("red", "blue", "green"))
legend("topright", c("BTF", "CBT", "NCBT"), pch=16,
col=c("red", "blue", "green"))
plot(split(BTF3a.p), main="Bifacial Thinning Flakes")
mtext("BTF = Fragments CBT = Cortex NCBT = No Cortex", side=1)
Marked planar point pattern: 230 points
Average intensity 12.8 points per square meter
Multitype:
frequency proportion intensity
BTF
87
0.378
4.86
CBT
49
0.213
2.74
NCBT
94
0.409
5.25
Window: polygonal boundary
single connected closed polygon with 18 vertices
enclosing rectangle: [982, 987]x[1015.5,
1022]meters
Window area = 17.91 square meters
Unit of length: 1 meter
plot(982, 982, xlim=c(982, 987), ylim=c(1015, 1022),
main="Bifacial Thinning Flakes", xlab="", ylab="",
axes=FALSE,
asp=1, type="n")
contour(density(BTF3a.p, adjust=.5), add=TRUE)
polygon(Win3a)
points(BTF3a.p, pch=20, cex=.75)
plot(density(BTF3a.p, adjust=.5), main="Bifacial
Thinning Flakes")
polygon(Win3a, lwd=2)
points(BTF3a.p, pch=20, cex=.75)
plot(density(BTF3a.p, adjust=.5), main="Bifacial
Thinning Flakes")
polygon(Win3a, lwd=2)
points(BTF3a.p, pch=20, cex=.75)
windows(10, 5)
V <- split(BTF3a.p)
A <- lapply(V, density, adjust=.5)
plot(as.listof(A), main="Bifacial Thinning Flakes")
Tab <- quadratcount(BTF3a.p, xbreaks=982:987, ybreaks=1015:1022)
quadrat.test(Tab)
# Warning: Some expected counts are small; chi^2 approximation may
#
be inaccurate
#
Chi-squared test of CSR using quadrat counts
# data:
# X-squared = 274.8859, df = 19, p-value < 2.2e-16
# Quadrats: 20 tiles (levels of a pixel image)
E <- kstest(BTF3a.p, "x")
plot(E)
N <- kstest(BTF3a.p, "y")
plot(N)
E
N
Spatial Kolmogorov-Smirnov test of CSR
data:
covariate 'x' evaluated at points of 'BTF3a.p'
and transformed to uniform distribution under
CSRI
D = 0.1101, p-value = 0.007611
alternative hypothesis: two-sided
Spatial Kolmogorov-Smirnov test of CSR
data:
covariate 'y' evaluated at points of 'BTF3a.p'
and transformed to uniform distribution under
CSRI
D = 0.2891, p-value < 2.2e-16
alternative hypothesis: two-sided
Gest(BTF3a.p)
plot(Gest(BTF3a.p), main="Nearest Neighbor Function G")
# Above poisson is clustered, below is regular
plot(envelope(BTF3a.p, Gest, nsim=39, rank=1), main="Nearest
Neighbor Envelope")
# The test is constructed by choosing a fixed value of r, and
rejecting the null
# hypothesis if the observed function value lies outside the
envelope at this
# value of r. This test has exact significance level alpha = 2 *
nrank/(1 + nsim).
plot(envelope(BTF3a.p, Gest, nsim=19, rank=1, global=TRUE),
main="Global Nearest Neighbor Envelope")
# The estimated K function for the data transgresses these limits
if and only if
# the D-value for the data exceeds Dmax. Under H0 this occurs with
probability
# 1/(M + 1). Thus, a test of size 5% is obtained by taking M = 19.
plot(alltypes(BTF3a.p, "G"))
plot(alltypes(BTF3a.p, "Gdot"))