
Oliver Nakoinz, Daniel Knitter
MOSAIC Summer School 2016
Interacting partners
Interacting systems
An approach and a set of methods that helps you to be explicit about the processes that caused the spatial distribution of your points (e.g. ceramic finds, settlements, graveyards, ...) [from pattern to process]
It uses the simplest possible form of spatial data: points/events in an area/region/space
Random point pattern
independent from space ...
Structured point patterns
are influenced by:
Some terminology: in case we have an random point pattern, Complete Spatial Randomness or an Independent Random Process prevails. The condition are:
Some terminology: in case we have an random point pattern, Complete Spatial Randomness or an Independent Random Process prevails. The condition are:
CSR means that the process is random not the resulting point pattern!
In point pattern analyses we test against CSR and different forms of specific processes to learn more about our pattern.
A structured point pattern violates CSR conditions:
First-order effects influence the probability of events being in any position of the region --> to trace such influences we investigate the intensity function (~ density) of the points
In case second-order effects are present, points are not independent from one another --> to trace such influences we investigate the distance distributions of the points
Events are independent
Events are not independent
Simple measures: mean, standard deviation, intensity (~ density)
download.file(
url = "https://raw.githubusercontent.com/dakni/mhbil/master/data/meg_dw.csv",
destfile = "2data/meg_dw.csv")
meg_dw <- read.table(file = "2data/meg_dw.csv",
header = TRUE,
sep = ";")
library(spatstat)
meg_pp <- ppp(x = meg_dw$x, y = meg_dw$y,
window = owin(xrange = c(min(meg_dw$x),
max(meg_dw$x)
),
yrange = c(min(meg_dw$y),
max(meg_dw$y)
),
unitname = c("meter", "meters")
)
)
plot(meg_pp)
mc <- cbind(sum(meg_pp$x/meg_pp$n),
sum(meg_pp$y/meg_pp$n)
)
stdist <- sqrt(sum((meg_pp$x-mean(meg_pp$x))^2 +
(meg_pp$y-mean(meg_pp$y))^2) /
meg_pp$n
)
plot(meg_pp)
points(mc)
library(plotrix)
draw.circle(x = mc[1], y = mc[2], radius = stdist, border = "red")
Global intensity = Number of points per area
## A = a * b
area.sqm <- diff(meg_pp$window$xrange) * diff(meg_pp$window$yrange)
area.sqkm <- area.sqm*10^-6
# area <- area/1000000
area.sqkm
## [1] 431.3529
## calculate intensity
intensity <- meg_pp$n/area.sqkm
intensity
## [1] 0.6189827
Local intensity
qc.meg <- quadratcount(X = meg_pp)
plot(qc.meg);
points(meg_pp, pch = 20, cex = .5, col = rgb(.2,.2,.2,.5))
Local intensity
Does the quadratcount indicates CSR?
Local intensity
Does the quadratcount indicates CSR?
To check we use a \(\chi^2\) test approach (remember: relation between observed (i.e. empirical) and expected (i.e. theoretical, here CSR) amounts of points in quadrants)
qt_meg <- quadrat.test(meg_pp)
qt_meg
##
## Chi-squared test of CSR using quadrat counts
## Pearson X2 statistic
##
## data: meg_pp
## X2 = 274.48, df = 24, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
## Quadrats: 5 by 5 grid of tiles
Local intensity
Meaning of bottom values:
plot(qt_meg)
Kernel density estimation
Kernel density estimation
meg_dens <- density(x = meg_pp,
sigma = 2500)
plot(raster(meg_dens))
Kernel density estimation
Kernel density estimation
We assume that the intensity of the point process is a function of the covariate (Z):
\[\lambda(u) = \rho(Z(U))\]
\(\lambda(u)\) can be regarded as location selection function --> it causes an inhomogeneous probability of points to be located in same areas of the study region
Calculation in R
is straightforward
cov_meg <- rhohat(object = "YOUR POINT PATTERN",
covariate = "YOUR COVARIATE RASTER",
bw = 100
)
pred_cov_meg <- predict(cov_meg)
cov_compare <- meg_dens - pred_cov_meg
Kernel density estimation
Kernel density estimation
\[R = \frac{observed~\overline{d_{min}}}{expected~\overline{d_{min}}}\]
\[R = \frac{\overline{d_{min}}}{\frac{1}{2\sqrt{\lambda}}}\]
meg_nn <- nndist(meg_pp)
mean(meg_nn)
## [1] 358.8447
nnE <- 1/(2*sqrt((meg_pp$n/area.sqm)))
nnE
## [1] 635.5222
R.meg <- mean(meg_nn)/nnE
R.meg
## [1] 0.5646453
hist(meg_nn)
abline(v=mean(meg_nn))
abline(v=median(meg_nn), lty=2)
cumulative frequency distribution of the nearest-neighbor distances
\[G(d) = \frac{\#(d_{min}(S_{i}) < d)}{n}\]
The function tells us what fraction of all n-n distances is less than d
meg_g <- Gest(meg_pp)
plot(meg_g)
cumulative frequency distribution of the nearest-neighbor distances of arbitrary events to known events
\[F(d) = \frac{\#(d_{min}(p_i,S) < d)}{m}\]
meg_f <- Fest(meg_pp)
plot(meg_f)
hmm...this can be advanced...in the Workshop :)
~ cumulative frequency distribution of all points within a certain radius
\[K(d) = \frac{\sum\limits_{i=1}^{n}\#(S \in C(s_i,d))}{n\lambda}\]
meg_k <- Kest(meg_pp)
plot(meg_k)
Definition
A networks are objects, in which elements (vertices) are connected by edges.
Networks are models, mapping certain facets of the real world.
Network theory has roots in geography and in social sciences
## Error in plot(n1): object 'n1' not found
library("igraph")
n1 <- graph( edges=c(1,5, 2,4, 1,3, 2,5,
3,5, 1,2, 3,4), n=6, directed=T )
plot(n1)
n1
## IGRAPH D--- 6 7 --
## + edges:
## [1] 1->5 2->4 1->3 2->5 3->5 1->2 3->4
E(n1)
## + 7/7 edges:
## [1] 1->5 2->4 1->3 2->5 3->5 1->2 3->4
V(n1)
## + 6/6 vertices:
## [1] 1 2 3 4 5 6
get.adjacency(n1)
## 6 x 6 sparse Matrix of class "dgCMatrix"
##
## [1,] . 1 1 . 1 .
## [2,] . . . 1 1 .
## [3,] . . . 1 1 .
## [4,] . . . . . .
## [5,] . . . . . .
## [6,] . . . . . .
n1 <- set_vertex_attr(n1, "label",
value =c("p1", "p2", "p3",
"Ppoint4", "Punkt 5", "6"))
plot(n1)
n2 <- make_full_graph(22)
plot(n2)
n3 <- make_tree(22, children = 3,
mode = "undirected")
plot(n3)
Delaunay graph
Construction rules for some neighbourhood graphs
The connections represent the liklyness of interaction
spdep
library("spdep")
wd <- "/home/fon/daten/analyse/mosaic"
setwd(wd)
set.seed(1242)
co.weapons <- read.csv("2data/
shkr-weapons.csv", header=TRUE,
sep=";")[sample(1:220,10),1:2]
##
## PLEASE NOTE: The components "delsgs" and "summary" of the
## object returned by deldir() are now DATA FRAMES rather than
## matrices (as they were prior to release 0.0-18).
## See help("deldir").
##
## PLEASE NOTE: The process that deldir() uses for determining
## duplicated points has changed from that used in version
## 0.0-9 of this package (and previously). See help("deldir").
coords <- as.matrix(coordinates
(co.weapons))
ids <- row.names(as.data.frame
(co.weapons))
wts <- co.weapons[,1]; wts[] <- 1
fs_nb_del <- tri2nb(co.weapons,
row.names=ids)
del <- nb2lines(fs_nb_del,
wts=wts, coords=coords,
proj4string = CRS(as.character(crs1)))
plot(del)
library(RANN)
fs_nb_soi <- graph2nb(soi.graph(fs_nb_del,
coords), row.names=ids)
soi <- nb2lines(fs_nb_soi, wts=wts,
coords=coords, proj4string =
CRS(as.character(crs1)))
plot(soi)
fs_nb_gabriel <- graph2nb(gabrielneigh
(coords), row.names=ids)
gabriel <- nb2lines(fs_nb_gabriel,
wts=wts, coords=coords,
proj4string = CRS(as.character(crs1)))
plot(gabriel)
fs_nb_relative <- graph2nb(
relativeneigh(coords),
row.names=ids)
relative <- nb2lines(fs_nb_relative,
wts=wts, coords=coords,
proj4string = CRS(as.character(crs1)))
plot(relative)
spdep-graph
to igraph-graph
n4nb <- nb2mat(fs_nb_del,
style="B", zero.policy=TRUE)
n4 <- graph.adjacency(n4nb,
mode="undirected")
plot(n4)
Centrality maps the structural importance of a node/edge in a network.
degree(n4)
## [1] 5 3 6 5 5 3 6 3 5 3
closeness(n4)
## [1] 0.07142857 0.06666667 0.08333333 0.07692308 0.07692308 0.05882353
## [7] 0.08333333 0.06250000 0.07692308 0.06666667
betweenness(n4)
## [1] 3.0000000 0.6666667 5.1666667 2.5000000 3.9166667 0.0000000 5.1666667
## [8] 0.0000000 3.9166667 0.6666667
edge_betweenness(n4)
## [1] 3.666667 2.833333 2.833333 2.000000 3.666667 4.166667 2.500000
## [8] 3.750000 2.666667 3.500000 4.083333 2.500000 1.833333 2.833333
## [15] 3.750000 1.833333 3.083333 4.083333 3.666667 3.083333 3.500000
## [22] 4.166667
spdep-graph
to igraph-graph
ceb <- cluster_edge_
betweenness(n4)
dendPlot(ceb, mode="hclust")
spdep-graph
to igraph-graph
plot(ceb, n4)
ABM comprises
The process
Actors can:
Monday, 5th of September
Tuesday, 6th of September
Wednesday, 7th of September
Thursday, 8th of September