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 51 0.33 0.36 0.04 0.0096 0.1743 0.8161
#> 2 56 0.32 0.40 0.04 0.0073 0.1111 0.8816
#> 3 60 0.30 0.35 0.04 0.0061 0.1471 0.8468
#> 4 56 0.34 0.32 0.03 0.0102 0.4473 0.5425
#> 5 54 0.35 0.36 0.04 0.0086 0.2587 0.7327
#> 6 53 0.29 0.28 0.04 0.0164 0.3457 0.6379
#> edge.weights.sum edge.weights.median edge.weights.median.abs.dev.
#> 1 206.0636 2.362922 1.455769
#> 2 290.4578 6.376574 5.908693
#> 3 365.1255 6.410656 6.666308
#> 4 254.2262 5.792639 4.294083
#> 5 238.8570 5.643078 3.140116
#> 6 225.6562 3.551193 3.217499
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.25 lower within
#> balance FALSE 0.46 within within
#> components.balance TRUE 0 within within
#> disturbance FALSE 0.58 within within
#> admixture TRUE 0.01 within lower
#> cohesion1 TRUE 0.01 higher within
#> cohesion2 TRUE 0.01 within within
#> edge.weights.sum TRUE 0 lower within
#> edge.weights.median TRUE 0 lower lower
#> edge.weights.median.abs.dev. TRUE 0.01 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.
