Generating artificial fragmentation graphs
The frag.simul.process function generates a pair of
spatial units containing fragmented objects with connection
relationships within and between these units. The next command creates
two spatial units populated with 20 initial objects (corresponding to
the “connected components” of a graph) which are fragmented into 50
pieces.
simul.g <- frag.simul.process(n.components=20, vertices=50)
This illustrates the simplest use of the
frag.simul.process function, which has several other
parameters to control the features of the simulation.
The number of initial spatial units is a crucial parameter, set using
the initial.layers parameter with “1” or “2”. This
parameter determines the method used to construct the graph and,
accordingly, the underlying formation process hypothesis.
If initial.layers is “1”, the fragmentation process is
simulated assuming that all the objects were originally buried in a
single spatial unit. The two clusters observed at the end of the process
are due to fragmentation and displacement.
- A single spatial unit is populated with the initial objects,
- the fragmentation process is applied,
- spatial units are assigned to the fragments,
- some fragments are moved as determined by the value of the
disturbance parameter.
If initial.layers is “2”, it assumes that the objects
were buried in two different spatial units, which were later partially
mixed due to fragmentation and displacement:
- two spatial units are populated with the initial objects
(components),
- the fragmentation process is applied,
- disturbance is applied.
The vertices and edges parameters are
related: at least one of them must be set, or both (only if
initial.layers is set to 1). Note that using both
parameters at the same time increases the constraints and reduces the
number of possible solutions to generate the graph. When there is no
solution, an error occurs and a message suggests how to change the
parameters.
The balance argument determines the number of fragments
in the smaller spatial unit (before the application of
the disturbance process). The components.balance also
determines the contents of the two spatial units by affecting the
distribution of the initial objects (components). Note that this
argument is used only when initial.layers is set to 2.
The aggreg.factor parameter affects the distribution of
the sizes of the components: this distribution tends to be more unequal
when aggreg.factor has values close to 1.
By default, fragments from two spatial units can be disturbed and
moved to another other spatial unit. However, the
asymmetric.transport.from can be used to move fragments
from only one given spatial unit.
Finally, the planar argument determines if the generated
graph has to be planar or not (a graph is planar when it can be drawn on
a plane, without edges crossing). This function requires to install the
optional RBGL package.
An example of a complete configuration of the function is:
frag.simul.process(initial.layers=1,
n.components=20,
vertices=50,
edges=40,
balance=.4,
components.balance=.4,
disturbance=.1,
aggreg.factor=0,
planar=FALSE,
asymmetric.transport.from="1")
An additional function is intended to simulate the failure of an
observer to determine the relationships between fragments. The
frag.observer.failure function takes a fragmentation graph
and randomly removes a given proportion of edges.
frag.observer.failure(abu.g12, likelihood=0.2)
Testing hypotheses
The versatile frag.simul.process function can generate
fragmentation graphs under multiple hypotheses about the initial
conditions (number of initial objects, number of initial spatial units,
etc.). Testing measurements on observed empirical data against
measurements made under these hypotheses can determine the most likely
initial conditions and fragmentation process.
Here, this is illustrated by comparing measurements from Liang Abu
layers 1 and 2 with measurements from simulated data under two
hypotheses about the number of initial spatial units (e.g., layers),
using the initial.layers parameter with two values, namely
one or two initial layers.
A fragmentation graph is generated for each
initial.layers value, using the parameters observed in the
Liang Abu layers 1 and 2 fragmentation graph. Setting the simulator is
made easier by using the frag.get.parameters function,
which takes a graph and computes a series of parameters that are
returned as a list.
params <- frag.get.parameters(abu.g12, layer.attr="layer")
params
#> $n.components
#> [1] 28
#>
#> $vertices
#> [1] 72
#>
#> $edges
#> [1] 52
#>
#> $balance
#> [1] 0.32
#>
#> $components.balance
#> [1] 0.29
#>
#> $disturbance
#> [1] 0.04
#>
#> $aggreg.factor
#> [1] 0.7
#>
#> $planar
#> [1] TRUE
#>
#> $edge.weights.sum
#> [1] 210.1143
#>
#> $edge.weights.median
#> [1] 1.849523
#>
#> $edge.weights.median.abs.dev.
#> [1] 0.7554912
# for H2:
test.2layers.g <- frag.simul.process(initial.layers=2,
n.components=params$n.components,
vertices=params$vertices,
disturbance=params$disturbance,
aggreg.factor=params$aggreg.factor,
planar=params$planar)
# for H1:
test.1layer.g <- frag.simul.process(initial.layers=1,
n.components=params$n.components,
vertices=params$vertices,
disturbance=params$disturbance,
aggreg.factor=params$aggreg.factor,
planar=params$planar)
Let us now generate not only one graph, but a large number of graphs
to statistically compare measurements in the empirical and simulated
graphs. The frag.simul.process function is set for the “two
initial layers” hypothesis and embedded into an ad hoc
function:
run.test2 <- function(x){
frag.simul.process(initial.layers=2, # note the different value
n.components=params$n.components,
vertices=params$vertices,
disturbance=params$disturbance,
aggreg.factor=params$aggreg.factor,
planar=params$planar)
}
The function is then executed a sufficient number of times:
test2.results <- lapply(1:100, run.test2)
The empirical values observed for Liang Abu layers 1 and 2 (red line)
can now be compared to the values measured in the simulated graph
generated under the hypothesis of two initial layers. This shows, for
example, that the empirical admixture value is slightly lower than the
simulated admixture values:
edges.res <- sapply(test2.results,
function(g) frag.get.parameters(g, "layer")$edges)
plot(stats::density(edges.res), main="Edges")
abline(v=params$edges, col="red")
Similarly, the empirical admixture value is lower than the simulated
admixture values:
admix.res <- sapply(test2.results,
function(g) frag.layers.admixture(g, "layer"))
plot(stats::density(admix.res), main="Admixture")
abline(v=frag.layers.admixture(abu.g12, "layer"), col="red")
Two functions (frag.simul.compare and
frag.simul.summarise) facilitate the execution of the
analytical process described above on the initial number of spatial
units. The frag.simul.compare function takes an observed
fragmentation graph, generates two series of simulated graphs
corresponding to two hypotheses on the number of initial spatial units
(H1 for one initial spatial unit and H2 for two initial spatial units),
and returns a data frame of measurements made on each series (including
the edge count, weights sum, balance value, disturbance value, admixture
value, and cohesion values of the two spatial units).
compare.res <- frag.simul.compare(abu.g12, layer.attr="layer",
iter=30, summarise=FALSE)
head(compare.res$h1.data)
#> edges balance components.balance disturbance admixture cohesion1 cohesion2
#> 1 53 0.28 0.26 0.04 0.0299 0.3798 0.5903
#> 2 52 0.27 0.27 0.04 0.0178 0.1896 0.7926
#> 3 58 0.29 0.30 0.04 0.0463 0.1783 0.7754
#> 4 60 0.32 0.32 0.04 0.0064 0.2858 0.7078
#> 5 55 0.28 0.22 0.04 0.0850 0.4211 0.4939
#> 6 53 0.32 0.33 0.04 0.0425 0.2450 0.7125
#> edge.weights.sum edge.weights.median edge.weights.median.abs.dev.
#> 1 207.1647 2.370370 1.171437
#> 2 202.6179 2.384503 1.553962
#> 3 250.2895 3.069665 2.528381
#> 4 350.6027 6.972864 3.803087
#> 5 193.3982 2.355198 1.421427
#> 6 187.1661 2.362922 1.444101
For each of these parameters, the frag.simul.summarise
function facilitates the comparison between empirical observed values
and simulated values generated for H1 and H2.
frag.simul.summarise(abu.g12, layer.attr="layer",
compare.res$h1.data,
compare.res$h2.data)
#> H1 != H2? p.value Obs. value/H1 Obs. value/H2
#> edges FALSE 0.16 lower lower
#> balance TRUE 0.03 within within
#> components.balance TRUE 0 lower within
#> disturbance <NA> <NA> <NA> <NA>
#> admixture FALSE 0.26 lower lower
#> cohesion1 TRUE 0 higher within
#> cohesion2 TRUE 0 lower within
#> edge.weights.sum FALSE 0.27 within within
#> edge.weights.median FALSE 0.15 lower lower
#> edge.weights.median.abs.dev. FALSE 0.45 lower lower
This function returns a data frame with four columns, containing, for
each parameter studied:
- whether the series of H1 values are statistically different to the
H2 series (Boolean),
- the p-value of the Wilcoxon test (numerical),
- whether the observed value is “within”, “higher”, or “lower” to the
interquartile range of values for H1,
- whether the observed value is “within”, “higher”, or “lower” to the
interquartile range of values for H2.
Note that the frag.simul.compare function can optionally
be set to execute and return the results of the
frag.simul.summarise function.
