Start Submission

Reading: Spatial Patterns and Grave Goods Differences at the Cemetery of Vedrovice (Czech Republic): ...


A- A+
Alt. Display

Case Study

Spatial Patterns and Grave Goods Differences at the Cemetery of Vedrovice (Czech Republic): A Resampling Approach to Identity Markers in the Early Neolithic


Petr Pajdla

Department of Archaeology and Museology, Faculty of Arts, Masaryk University, CZ
X close


The spatial patterns and grave good associations at the cemetery of Vedrovice – Široká u Lesa (Czech Republic) are explored using computationally intensive formal methods. Resampling approaches to both intra-site spatial arrangement and the distribution of grave goods are explored in relation to layers of identity represented by different buried body attributes – sex, age, and geographical origin. Polished stone adzes, grinding tools and ceramic bowls are non-randomly distributed grave goods. The investigation of differences among the non-randomly distributed artifacts and attributes of burials indicates that the most significant are dichotomies between sexes, accentuated by both origin and age. Spatial analysis shows clear trends toward clustering between burials, and the emerging clusters are further explored for similarities in burial attributes.

How to Cite: Pajdla, P., 2023. Spatial Patterns and Grave Goods Differences at the Cemetery of Vedrovice (Czech Republic): A Resampling Approach to Identity Markers in the Early Neolithic. Journal of Computer Applications in Archaeology, 6(1), pp.1–15. DOI:
  Published on 15 Feb 2023
 Accepted on 02 Jan 2023            Submitted on 24 Aug 2021


The Cemetery of Vedrovice is one of the earliest and most intensively studied Linearbandkeramik (LBK) cemeteries in the Czech Republic. Recent articles use evidence from the cemetery in a broader discussion of post-marital residence patterns (Hrnčíř, Vondrovský and Kvĕtina 2020), sexual inequalities, and gender in the Neolithic (Augereau 2022; Bickle 2020; Masclans Latorre, Bickle and Hamon 2020). Focus of this article is an exploration of the relationship between three proposed layers of identity, biological sex and age of the body determined from skeletal remains by anthropological methods and the origin of the body based on the strontium isotope signal. The question is whether the structure observed in these three categories is mirrored in the grave good associations and spatial organization within the cemetery. The analysis is set into a broader discussion of gender relationships and identity in Neolithic communities.

Robb and Harris (2017: 128) noted a lack of debate in archaeology on gender for the European Neolithic compared to other periods. The authors present gender in the Neolithic as being of a very different nature than gender in later prehistory. They argue that it is not clearly expressed and biological differences between sexes are not transformed into binary cultural categories (Robb and Harris 2017: 140–142). Considering the Neolithic burial evidence across Europe, authors find that it ‘shows very little gender distinction; when it does, it is rarely a clear-cut binary categorization but a matter of degree and overlap, sometimes polythetically combined with factors such as age and localness’ (Robb and Harris 2017: 138). In response, Bickle (2020: 14) argues that ‘biological sex may not have defined gender during the LBK’ although the ‘sexed bodies did matter’.

I focus mainly on investigating non-random patterns in artifact co-occurrences and spatial arrangement in the cemetery in relation to the proposed layers of Neolithic identity, i.e., sex, age, and localness of the body. My question is how grave goods relate to the categories mentioned above and whether there is a consistent relationship between these categories creating the identity and structuring spatial arrangement of the cemetery.

The chosen approach is quantitative, embracing resampling techniques and computer-intensive methods in general. Resampling techniques were introduced in mortuary studies by Manly (1996) but the application is still limited (Sosna, Galeta and Sládek 2008: 343). As noted by Sosna, Galeta, and Sládek (2008), the most important benefit of resampling is its ability to distinguish random and non-random patterns. This is crucial for mortuary analysis because the studied burial grounds, including the one presented here, usually consist of a small number of burials, and identifying which observed patterns are caused by mere chance and which are truly present is problematic.

1 Material and Methods

The Vedrovice site is situated in south Moravia (Czech Republic, see Figure 1). The earliest known finds and archaeological explorations are from the 19th century, but excavations were conducted in the second half of the 20th century. There are two known burial areas, Za Dvorem and the larger Široká u Lesa, and a settlement area. For more details on the site, see Masclans Latorre, Bickle and Hamon (2020) or the original publication by Podborský (2002). The material analyzed here comes from the Široká u lesa cemetery in Vedrovice. Burials at this cemetery are radiocarbon dated to a period of five to six generations from 5300 to 5200 BC (Pettitt and Hedges 2008).

Topographic map of Central Europe with the location of the Vedrovice site marked
Figure 1 

The location of Vedrovice – Široká u Lesa cemetery indicated on a map of the Czech Republic and neighboring countries.

Originally, this burial ground consisted of at least 96 LBK burials and one (later) Lengyel culture cremation (Ondruš 2002, 21). Many of the LBK burials were at least partially disturbed by later agricultural activities, especially deep plowing. Some were disrupted by Lengyel culture settlement pits in the area of the cemetery, and 11 identified burial pits were probably excavated by amateur archaeologists during the 19th century (see the plan of the cemetery, Figure 2A). Any destroyed or heavily disturbed burials were excluded from the analysis of grave goods co-occurrence due to possible missing artifacts and/or bodies. This leaves a sample of 75 well-preserved burials in the analysis of grave-goods co-occurrence. In the analysis of spatial relationships, I include disturbed burials and empty burial pits (excavated in earlier periods and identified as LBK) because their locations, and not the contents, are of interest. Spatial analysis is thus performed on the whole sample of 96 LBK burials.

Plans of the cemetery arranged in one figure as small multiples
Figure 2 

Plan of the cemetery. (A) Preserved burials in white, disturbed burials in gray. Burial numbers as in the original publication by Ondruš 2002, 16. Approximate extent of the excavation (a window) in gray. (B) Distribution of burials based on localness. Panels for the local and non-local categories with kernel density estimate. (C) Distribution of age groups. Panels for juveniles (juv.), adults (ad.), and mature bodies (mat.) with KDE. (D) Distribution of sex groups. Panels for non-adults (n. a.), female-sexed bodies (F), and male-sexed bodies (M) with KDE for each of the categories.

Most of the data were obtained from the original publication (Podborský 2002), and the determination of the sex and age of the bodies were corrected based on the anthropological analysis by Dočkalová (2008) and Bramanti (2008). Carbon, nitrogen, and strontium isotope ratios were collected in Whittle et al. (2013) and Richards et al. (2008). Further details on the use and function of the artifacts were gathered from Masclans Latorre, Bickle and Hamon (2020). The method to derive the localness of various bodies based on the Sr isotopes was adapted from Hrnčíř, Vondrovský and Kvĕtina (2020).

Each burial was initially characterized by more than 20 variables recording the number of various artifacts present in the graves. Some of these original variables are omitted in the analysis due to their rare occurrence; for instance, a spondylus arm ring is present in only a single burial, rendering it useless in establishing relationships between pairs of burials. Only grave goods occurring in more than 5% of the analyzed burials are used in the analysis. This step reduces the number of input variables to 14. These include variables noting quantities of various types of pottery vessels in the graves, i.e., bottles, globular pots, and bowls, counts of chipped stone tools made of either local or non-local raw materials, and the presence of polished and beveled stone adzes, grinding slabs, pebbles, bone tools, and different kinds of personal ornaments. The personal ornaments include beads made of spondylus shells, spondylus pendants, and oval-shaped buckles. Last but not least, the presence of ocher pigment in the grave is also noted. Quantities of various grave goods are simplified into a binary matrix indicating the presence or absence of a given trait. Effects of grave goods counts are removed because some numerous artifact categories, e.g. spondylus beads, would distort the results.

Another group of variables is derived from retrieved human skeletal remains. These include biological sex, age simplified into age groups, and an indicator variable as to whether the strontium signal in the remains is of local origin. The variable representing the recorded sex is simplified into the following categories: female-sexed body (abbreviated F), male-sexed body (M), non-adult (n. a.), and undetermined sex (ind., see the spatial distribution in Figure 2D). The age of the body is also simplified into three categories: juveniles (juv.) comprised of the groups Infans and Juvenis, adults (ad., all bodies aged as Adultus), and mature bodies (mat., anthropological categories Adultus-Maturus, Maturus, and Senilis) and bodies of undetermined age (ind., see Figure 2C). The last category is the indicator variable for the localness of the buried person. It derives from the local strontium range which is defined as two standard deviations from the mean of Sr ranges of bodies with an age category determined as Infans (method adapted from Hrnčíř, Vondrovský and Kvĕtina 2020). Bodies falling within this range are labeled as local, outside the range as non-local; where data on Sr values is missing, the origin is not determined (ind., see Figure 2B). This approach is based on the premise that children were born locally and their Sr signal thus gives approximate boundaries for local Sr range. As mentioned above, I consider these three variables – sex, age, and origin (localness) – to be indicative of an individual’s identity in the Neolithic.

1.1 Randomness of Artifact Combinations

Non-randomness of various artifact combinations in burials is determined to relate artifact categories to the sex, age, and origin of the buried bodies. This is achieved by conducting a permutation test on the presence/absence data matrix. This method based on generalized Monte Carlo tests was introduced into archaeology by Manly (1996) and adapted by Sosna, Galeta, and Sládek (2008). My approach is based on both of these works. The original artifact co-occurrence matrix, i.e. how many times each combination of artifact categories is co-present, was permuted n = 9,999 times with constant row and column totals. With this procedure, I simulate virtual burial grounds where both (1) the overall distribution of artifact types across the whole burial ground, and (2) the count of artifact types in individual burials are identical to the Vedrovice – Široká u Lesa site.

To examine which associations of artifacts are infrequent or even not present and which tend to appear together more often than it would be expected by random chance, the co-occurrences of observed and simulated data are summarized using two statistics suggested by Manly (1996: 475). The S statistic is a measure of overall deviation between the expected and observed number of artifact occurrences, while the v statistic is a measure of the difference between co-occurrences of individual artifacts with other artifacts.

This step introduces a multiple testing problem. There are 14 tests to determine whether any of the variables (types of artifacts) co-occur by chance (null hypothesis). At the 5% alpha level, 5% of the tests will be significant by definition. Thus, one of the significant tests (precisely 0.7 tests) will be a false positive (type I error). Common approaches to control for the multiple hypothesis testing problem, for example Bonferroni-type correction, minimize the type I error rate (false positive rate) to 5%. That means null hypothesis is rejected even though the artifact type is distributed randomly. Type II error rate increases by this approach (false negative rate), i.e. a test fails to reject the null hypothesis that the distribution of an artifact follows a random pattern even though the artifact is distributed non-randomly. In the case of non-random co-occurence of artifact types I find it more acceptable to risk that false positives are determined as non-randomly co-occuring than not to identify non-randomly distributed artifact types (cf. Perneger 1998). Thus, I do not correct for the multiple hypothesis tests problem.

To assess whether differences between artifacts that are non-randomly distributed reflect in relationships between groups of people based on different sex, age, and origin assignments, I calculate mean intergroup Euclidean distance (ED) in multidimensional space settings defined by counts of non-randomly distributed artifacts. The observed mean ED for each of the groups is compared to the mean ED distribution derived by random reshuffling of the group category labels between the burials 9,999 times with the constant probabilities of appearance of the given category label. For each of the categories, the deviation of the mean ED from its randomly generated distribution is expressed as the fraction of random values of ED larger than the one observed. The deviation is deemed significant if less than 5% of random values are smaller or larger (both tails of the distribution are taken into account) than the observed value. The process is a variant of the one described for resampling in the case of artifact combinations, and a similar approach is taken by Sosna, Galeta, and Sládek (2008: 346–347).

1.2 Spatial Analysis

Spatial patterning in prehistoric burial grounds is often examined to find any internal structure, arrangement, or possible grouping. I approach the task as a point pattern analysis problem given that a burial pit (an event in the point-pattern analysis terms) is simplified as its centroid location in two-dimensional space. A polygon of a simplified excavation extent is used as the study area (a window).

1.2.1 Spatial randomness functions

The overall structure of the cemetery is examined using methods utilizing distances between neighboring burials. These methods are based on distances between pairs and threesomes of neighbors and characterize the point pattern based on second- and third-order properties (Nakoinz and Knitter 2016: 130–135; O’Sullivan and Unwin 2010: 107). These are in contrast to the assumption underlying first-order effects that any point location is independent of other point locations and different parameters influence its location (O’Sullivan and Unwin 2010: 126). In the spatial analysis of settlement patterns, individual locations of settlements can depend on several parameters, e.g., elevation, soil type, etc. In the case of a grave placement in a burial ground, it is presumed that the choice of the grave location is either random or influenced predominantly by other (previous) grave locations and is largely independent of outside parameters. In other words, first-order effects are of little to no importance from my point of view. The second-order properties and third-order properties are thus considered major factors in the location of the graves.

The analysis utilizes the G, F, Ripley’s K (Ripley 1976; Sayer and Wienhold 2013) and T (Schladitz and Baddeley 2000) spatial randomness functions. These functions are based on nearest neighbor distances and are commonly used to characterize point patterns in archaeology (Nakoinz and Knitter 2016: 135–144). The G function is a metric of how close burials (i.e. events) are and summarizes the cumulative frequency distribution of the nearest neighbor distances. Because of its dependence on the intensity of burial locations in the study area, it is subject to changes in its size (window size). The F function, an empty-space function, is based on simulating a point pattern and calculating the distances between these random locations and the observed locations. It is a measure of how far observed burials are from any arbitrary location (O’Sullivan and Unwin 2010: 132–135). In contrast to the G function, Ripley’s K function is resilient to changes in the window size and works with all the proximal locations up to a certain threshold, not only the nearest neighbors (Nakoinz and Knitter 2016: 138). Negre, Muñoz and Barceló (2018) find it the most useful in most archaeology cases. Last but not least, the T function is a third-order parallel to the K function (Schladitz and Baddeley 2000).

The given statistics for the independent random process are estimated using Monte Carlo simulation with 999 iterations. Although a smaller number of permutations should be sufficient to achieve credible results (Baddeley, Rubak and Turner 2015: 384), I follow the approach taken by Sayer and Wienhold (2013: 78). The ten most extreme realizations of the simulation are discarded as outliers, so the significance level is 0.01 (cf. Baddeley, Rubak and Turner 2015: 410). The remaining realizations of given point processes are used to construct envelopes around function curves derived from observed data.

To explore the frequency of neighboring burials falling into the same categories, I explore several methods. A modification of the K function for a multitype point pattern does not take into account all the neighboring individuals in the point pattern but pairs of observations of the same value within the given category are explored. For instance, in the case of localness, bodies designated as being of a local origin and non-local bodies are considered as two separate point patterns. In the same way that the K function is indicative of the distances at which clustering or repulsive processes take place, the multitype K function shows these distances for the given pairs of values within a given categorical variable. The envelopes are constructed using the same procedure as described previously.

1.2.2 Neighborhood graphs

To further explore the neighborhood structure, I construct two neighborhood graphs (Figure 3). Delaunay triangulation (Figure 3A) and its subgraph, a Gabriel graph (Figure 3B, Gabriel and Sokal 1969) with burial locations as nodes and edges connecting perceived pairs of neighboring burials is used. In the Gabriel graph, an edge is created between two burials (nodes) if a circle with a diameter equal to the distance between the nodes does not contain any other points. Similarly, in the Delaunay triangulation, edges between three nodes (a triangle) lying on a circumference of a circle exist if there are no other points present in the circumcircle. The Delaunay triangulation is regarded here as the maximum number of neighbors per individual burial case scenario, where, especially at the rims of the cemetery, the edges connecting the neighboring burials tend to be quite long. The Gabriel graph is the opposite with only the closest pairs of burials forming neighbors.

Plans of the cemetery, one with a network based on Delaunay triangulation and the other on Gabriel graph
Figure 3 

Neighborhood based on graph structures – (A) Delaunay triangulation and (B) Gabriel graph.

For each burial, the mean number of neighbors of a given category defined by sex, age and localness allocations is determined and recorded. To check how common the observed value is, a permutation test is performed where the individual categories of sex, age, and origin are randomly assigned while preserving the overall proportion of given categories at the cemetery and the neighborhood structure determined by the graph. The observed values of neighbor categories are compared to distributions derived from 999 permutations of the simulated data.

1.2.3 Percolation analysis

Another approach considered for establishing neighborhood relationships between graves is constructing buffer zones around each of the burials. To avoid determinism in designating a certain radius of the buffer zone around each of the individuals, I utilize percolation analysis, a data-centric method used to identify clusters of points based on spatial separation measured by Euclidean distance (Maddison and Schmidt 2020: 269). Although the introduction of percolation analysis in archaeology is relatively recent, its effectiveness and usefulness in studies with varied scales is well-proven (Maddison and Schmidt 2020). Percolation analysis is based on a city clustering algorithm (Rozenfeld et al. 2008) and its generalization which employs Euclidean distance (Rozenfeld et al. 2011). Points falling within a given distance threshold (radius) are added to a cluster. This process is repeated for a set of thresholds while recording maximum cluster sizes and possibly other statistics. This makes it possible to identify distances at which clusters exhibit rapid growth (absorbing numerous nodes) and distances at which the clusters remain relatively stable. For more details, refer to Maddison and Schmidt (2020) and the associated R package (Schmidt and Maddison 2020).

Percolation analysis identifies threshold distances at which maximum cluster size is growing, and clusters formed at these distances are further explored. For each cluster, the number of individual bodies of a given category are recorded. The same procedure is repeated for 999 resamples with the given categories assigned at random with constant probabilities. The observed numbers of occurrences are compared against the simulated distribution to check how extreme it is. The observed value is regarded as occurring non-randomly if it is larger than 95% of simulated values (right tail of the distribution) or smaller than 5% of simulated values (left tail of the distribution). The first case is evidence for overrepresentation, i.e., there are more bodies of a given category than expected under a random process. In the latter case, the given category is underrepresented, i.e., there are fewer bodies of a given category in the given cluster than would be expected under the assumption of a random process taking place.

Data is analyzed in an R environment for statistical computing and graphics (R Core Team 2022) using functions of several different packages; the most frequently used are packages from the tidyverse family (Wickham et al. 2019). Packages sf (Pebesma 2018), spatstat (Baddeley, Rubak and Turner 2015), spdep (Bivand and Wong 2018), and concaveman (Gombin, Vaidyanathan and Agafonkin 2020) are used to perform various spatial analysis tasks; igraph (Csardi and Nepusz 2006) to construct and analyze networks; vegan (Oksanen et al. 2020) and randomizr (Coppock 2022) in the randomization experiments.

2 Results

2.1 Artifact Combinations

Counting the co-occurrences of analyzed variables on simulated data matrices and comparing the distributions of the statistics proposed by Manly (1996) with the values calculated for the original data set establishes whether the overall pattern of artifact distribution in the graves is non-random and further identifies artifacts that deviate from the random pattern. The overall deviation from expected co-occurrences, i.e., the S statistic (Manly 1996: 475), is significantly large (S = 1.88, p = 0.03) at the 0.05 alpha level. This suggests that there is evidence that the overall distribution of artifacts in different graves is not random (see Figure 4).

Probability density distribution of the simulated and observed S statistic
Figure 4 

Probability density distribution of simulated and observed S statistics. Note the value of the S statistic and a p-value in the top right corner.

Examining distributions of 14 studied artifact types, I find that at the 0.05 significance level only adzes (v = 2.93, p = 0.048), grinding tools (v = 1.54, p = 0.03) and ceramic bowls (v = 2.65, p = 0.049) are distributed non-randomly (see Figure 5). At the 0.05 alpha level, one of these is by definition a false positive. Close to the 0.05 alpha level are distributions of bone tools (v = 1.82, p = 0.054). Among the grave goods found to be randomly distributed in the graves across the burial ground are ocher pigment, pebbles, both local and non-local lithics, globular pots, bottles and special types of pottery vessel shapes, spondylus beads, buckles, pendants and already mentioned bone tools.

Small multiple plot of probability density distributions of simulated and observed v statistics for different grave goods
Figure 5 

Probability density distribution of simulated and observed v statistics for different grave goods. Note the value of the v statistic and associated p-value in the top right corner.

The relationships between artifact types frequently occurring and missing together in graves are visualized in Figure 6 with clusters established using the fast-greedy algorithm (Newman and Girvan 2004) and the strength of the given association is expressed by the thickness of the edge connecting the nodes. Artifacts most frequently present with other artifacts (see Figure 6A) are bowls and adzes both with five links and grinding tools, spondylus beads, and non-local lithics all with 4 links. The association is strongest between spondylus beads and (a) ceramic bottles and (b) pigment; and adzes and (a) non-local lithics and (b) bone tools. Artifacts frequently absent from other artifacts (see Figure 6B) in the graves are local lithic tools, globular pots, and spondylus buckles with seven, six, and five links respectively. The strongest antagonist relationships are between globular pots and (a) bowls, (b) adzes and (c) bone tools; and local lithics and (a) pebbles, (b) bottles, (c) bowls, and (d) ochre pigment.

Networks of grave goods frequently present together and missing from each other in the graves
Figure 6 

Networks of variables (A) frequently present together, and (B) excluding each other. 1 – bottle, 2 – bowl, 3 – globular pot, 4 – polished stone tool, 5 – local lithics, 6 – non-local lithics, 7 – pebble, 8 – grinding slab, 9 – bone tool, 10 – U-shaped pendant, 11 – spondylus beads, 12 – O-shaped buckle. Images adapted from Podborský (2002).

Artifacts non-randomly distributed across the cemetery, i.e. adzes, grinding stone tools, and ceramic bowls are further used to check whether significant differences exist between their occurrences in different groups defined based on origin, sex, and age. This is done by calculating the mean Euclidean distance between the groups in three-dimensional space settings defined by counts of the non-randomly distributed artifacts and assessment whether the observed ED occurs by chance or not using the Monte-Carlo simulation.

In the case of localness, all of the relationships between groups of local, non-local, and undetermined individuals are well explained by random processes. In other words, the non-random distribution of adzes, grinding stone tools and ceramic bowls are not reflected in the structure given by the origin of the bodies. Similar results apply to the age categories except for pairs of juvenile bodies. In this case, the observed mean ED is smaller than 95.9% of simulated values (see Table 1), i.e. in the space defined by non-randomly occurring artifacts, the juvenile bodies are closer than would be expected. Contrary to origin and age categories, most of the mean intergroup ED observations in the sex category do not fit the expected values. In general, this means that the space defined by adzes, grinding tools, and bowls is structured by sex. The mean ED between pairs of female-sexed bodies is smaller than 99.9% of simulated values. Similarly, the mean ED is small for female-sexed and non-adult bodies. On the other hand, only 1% of simulated values is larger than the mean ED observed for pairs of male-sexed bodies. The distances between female-sexed and male-sexed bodies (only 0.1% of simulated values is larger) and male-sexed and non-adult bodies (1.2%) are also larger than expected. Bodies of undetermined sex category have larger distances than expected with female-sexed and non-adult bodies (see Table 1 for the percentages).

Table 1

Summary table for observed and simulated mean Euclidean distances of categories within the sex and age variables in a space defined by the counts of non-randomly distributed artifacts, i.e. adzes, grinding tools and bowls. For abbreviations, see the labels in Figure 2.


Age juv. – juv. 0.32 95.9

Sex F – F 0.29 99.9

Sex F – ind. 2.05 1.9

Sex F – M 1.03 0.1

Sex F – n. a. 0.37 99.9

Sex M – M 1.11 1.0

Sex M – n. a. 1.00 1.2

Sex n. a. – ind. 1.93 2.57

2.2 Spatial Patterning and the Neighborhood Structure

The results of functions based on nearest neighbors are reported in Figure 7. The curve for function G (Figure 7A) falls within the boundaries of the envelope, i.e., the gray shaded area. This means that there is no evidence of significant deviations from the independent random process. Function F (Figure 7B), on the other hand, shows that there are fewer neighboring points present at a radius of around 5 to 10 meters than would be expected under the assumption of the independent random process. This means that at this distance, there are empty spaces present in the point pattern and a clustering process is taking place. Similarly, the K and T functions (Figure 7C–D) show clear deviations from the expected curves and both exceed the boundaries of the given envelopes at distances around 5 meters. This means that the number of neighbors at this and larger distances are significantly higher than expected under complete spatial randomness. The overall point pattern thus shows clear signs of the clustering process taking place.

This suggests that the locations for new burials were not chosen at random but were influenced by the positions of previous burials. The pattern begins to cluster at a distance of around 5 meters, i.e., beginning at this distance, the probability of finding a neighboring burial is much higher than it would be if the burials were placed randomly in space. The community burying their dead at the Široká u lesa cemetery in Vedrovice did so in an organized manner taking into account previous burials at the same place, meaning their positions were known, probably marked somehow on the surface. This would imply the existence of non-perishable marks on the surface of the cemetery and/or their maintenance from time to time.

The results of a percolation analysis reveal a more detailed view of the clustering processes taking place; see Figures 8 and 9. The distances at which the relative stability of clusters emerges are labeled in Figure 8. The curve reveals a steady growth of clusters up to a radius of 4.1 meters. At the small distances (see Figure 9, especially panels A and B) clusters are formed mostly by pairs of neighboring burials. With growing percolation distance, larger clusters emerge (Figure 9, clusters 1 and 2 on panels C to D). Between 4.1 and 4.6 meters, there is a jump (Figure 8) toward a single large cluster (Figure 9, cluster 1 on panels E and F). This point of cluster growth is revealed by the point pattern analysis as well where there is a tendency of clustering visible in the case of F, K, and T functions (Figure 7) at the distance around 5 meters.

Curves for G, F, K and T as functions of a radius around burials
Figure 7 

G, F, K, and T as functions of a radius r. Labels in panel D, dashed line Ttheo – theoretical function under an assumption of complete spatial randomness; solid line Tobs – estimate for a function as observed in the given data set; gray area – an envelope of the curve of the theoretical function.

Curve showing a proportional size of a largest cluster as a function of distance
Figure 8 

Proportional size of the largest cluster at any given percolation distance. Most eminent threshold distances at which the cluster growth stabilizes are labeled.

Small multiple with plans of the cemetery for different percolation distance thresholds with given clusters
Figure 9 

Individual panels A to F show clusters formed at given percolation threshold distances (see Figure 8). Clusters are highlighted in white, only clusters composed of at least 4 burials are labeled with numbers 1 to 6 in descending order by the size of a cluster.

2.2.1 Spatial randomness functions

The results of the multitype K function show deviations from random processes for many of the combinations in the sex, age, and origin categories. The distances at which these deviations start occurring are in most cases rather large, well over 10 meters, and considering the fact that at the 9.2 meters radius around individual points the percolation analysis showed the existence of a single large cluster encompassing almost all the burials (Figure 9), these observations correspond to the deviations observed in the analysis of a single K function.

2.2.2 Neighborhood graphs

In the case of neighborhood structure derived from Delaunay triangulation (maximum number of neighbors per a single burial scenario, Figure 3A) and Gabriel graph (minimum number of neighbors scenario, Figure 3B), most of the observed relationships between categories defined by bodies of different age, sex and origin are not significantly different from those derived by randomization experiments. Given the neighborhood construction methods, most of the observed pairs of neighboring burials are thus sufficiently explained by chance. There are several exceptions: in the Gabriel graph, there is a relationship in the age category between pairs of mature bodies and in the localness category between bodies of local and non-local origin. In both of the graphs, the Delaunay triangulation and the Gabriel graph, there is a relationship between female-sexed and male-sexed bodies.

2.2.3 Percolation analysis

Clusters formed at different percolation distances show several deviations from the expected number of occurrences of given categories at different scales. The results for observed counts larger than 95% of simulated values are summarized in Table 2. Results for the lower tail, i.e. where the observed value falls under 5% of the simulated values, are due to the length of the table deposited together with the code in an online repository (see section Reproducibility). Across different radii, there is an overrepresentation of individuals with undetermined age or localness assignments present in the eastern part of the cemetery. This includes cluster 5 at the distance of 2.8 meters (Figure 9B), clusters 3 at the distance of 3.4 and 5.6 meters (Figure 9C and F), and cluster 2 at the distance of 6.3 and 8.8 meters where there are more bodies of undetermined age present than we would expect under the assumption of a random process taking place. Under the random process we assume that bodies with undetermined sex, age and origin categories are distributed randomly across the area of the cemetery. Similarly, in cluster 4 at the distance threshold of 4.1 meters (Figure 9D), cluster 3 at 5.6 meters (Figure 9F), and cluster 2 at 6.3 and 8.8 meters there are more bodies of undetermined origin present. On the other hand, in the eastern part of the cemetery, there are fewer bodies determined as having non-local origin than would be expected. These occur in clusters 4 (radius 4.1 meters, Figure 9D), 3 (radii 4.6 and 5.6 meters, Figure 9, panels E and F respectively), and 2 (radii 6.3 and 8.8 meters). Furthermore, there are fewer female-sexed bodies in clusters 3 (radius 5.6, Figure 9F) and 2 (radii 6.3 and 8.8 meters). These observations are strong evidence of a bias in the eastern part of the cemetery caused by the concentration of disturbed burials in this part of the cemetery (Figure 2A).

Table 2

Summary table for observed and simulated occurrences of categories within the sex, age and origin variables. Only observed values that are in the top 5% of the simulated distribution are shown (i.e. the upper tail). For abbreviations, see the labels in Figure 2.


2.1 1 Age ad. 4 2.5

2.8 5 Age ind. 3 2.2

3.4 3 Age ind. 3 4.1

3.4 5 Age juv. 4 0.3

4.1 2 Age juv. 9 3.8

4.1 2 Sex n. a. 8 4.9

4.1 4 Origin ind. 5 4.3

4.6 4 Sex F 4 2.1

4.6 4 Age mat. 3 4.0

5.6 3 Origin ind. 6 3.9

5.6 3 Age ind. 4 2.4

5.6 4 Sex F 6 0.2

5.6 4 Age mat. 4 3.0

6.3 2 Origin ind. 6 4.4

6.3 2 Age ind. 4 3.3

6.3 3 Sex F 6 0.5

6.3 3 Age mat. 4 1.8

8.8 2 Origin ind. 6 4.7

8.8 2 Age ind. 4 2.9

The group of burials in the southwestern part of the cemetery (cluster 3 at 3.4 meters radius and clusters 2 at 4.6 and 5.6 meters radii, Figure 9, panels D to F) is characterized by the underrepresentation of undetermined bodies in the age and sex categories. This further highlights the fact that bodies of determined and undetermined sex, age, and origin categories are not distributed uniformly across the space of the cemetery. Another group consisting of fewer bodies of determined age and sex groups is situated in the south of the cemetery, it is cluster 2 at the 3.4 meters percolation radius (Figure 9C).

The group of burials in the northern part of the cemetery is characterized by more female-sexed individuals assigned to the Maturus age group at the percolation distances of 4.6, 5.6 (clusters 4, Figure 9, panels E and F respectively), and 6.3 meters. Cluster 2 at the percolation distance of 4.1 meters in the southeastern part of the cemetery (Figure 9D) is characterized by an overrepresentation of juvenile (non-adult in case of the sex determination) bodies. Male-sexed bodies are, on the other hand, underrepresented in this cluster.

In the central part of the cemetery at the percolation distance of 2.1 meters in cluster 1, there is an overrepresentation of adult bodies (Figure 9A). In the same part of the cemetery at the distance of 2.8 meters in cluster 2 (Figure 9B) bodies of non-local origin are underrepresented. Furthermore, at the 3.4 meters radius in cluster 1 (Figure 9C), there are fewer juvenile (non-adult) bodies than we would expect under the assumption of a random process. On the other hand, in cluster 5 at the same percolation distance juvenile bodies are overrepresented.

Conclusions and Discussion

The goal was to identify non-random differences between burials both in terms of the artifact contents of the graves and position of the burials within the cemetery. The buried bodies were divided into categories based on three variables commonly connected to identity in the Neolithic period. These categories are sex (male-sexed, female-sexed, non-adult, and undetermined bodies), age (juvenile, adult, mature, and undetermined bodies), and origin (local, non-local and bodies of undetermined origin). The observations are based on the extensive use of resampling approaches that allowed distinction between patterns occurring solely by chance and patterns significantly deviating from randomness. These methods prove essential in analyzing small samples like the one presented here.

Establishing whether artifacts in burials occur randomly or not enables identification of the artifact types that are somehow different from the rest of the assemblage. Although the deposition of an artifact in a burial is, in most cases, a deliberate practice, from an analytical point of view it is important to distinguish between artifacts that account for the noise, i.e. are deposited randomly across the burials, and the signal, i.e. a very specific choice is made when depositing the artifact within the given burial. In the case of Vedrovice – Široká u Lesa cemetery the distributions of ceramic bowls, adzes and grinding stone tools are found to follow non-random patterns. In the next step, I identified what associations between the non-randomly distributed artifact types and groups of bodies defined by different sex, age and origin categories were occurring at random. This was achieved by measuring mean intergroup Euclidean distances in the space defined by the counts of non-randomly occurring artifact types. Sex of the body is the most structuring variable with significant deviations from values expected under an assumption of a random processes taking place. The relationships in the origin category (i.e. localness of the body) were all well explained by random processes. In the age category, only the pairs of juvenile bodies deviated from random processes. This means that sex classification of the body, although sometimes being accentuated by age category, structures the artifact assemblage in the grave the most.

The spatial organization of the cemetery was explored using point pattern analysis techniques. The methods based on the cumulative frequency distribution of nearest neighbors (G function) did not show deviations from complete spatial random processes. This is not surprising because of its strong dependence on the intensity of burials across the studied area, i.e. its size. On the other hand, the generalization for empty space points (F function) and the methods including all the distances between burials in the given area (K and T functions), indicated that clustering processes are taking place at the 5 m distance threshold. The K function for multitype point patterns was also used to identify clustering processes between burials of different sexes, age categories, and origins of the body. This is not useful in most cases because the distances at which clustering processes are suggested imply the existence of a single cluster across the cemetery and a finer scale is of much greater interest for the analysis of the intra-site relationships. The nearest-neighbor methods are thus a great tool in analyzing the overall structure of the cemetery, but on a small scale, they do not prove useful for discovering relationships between studied layers of identity.

Thus, to identify whether close burials are more similar in terms of the categories outlined above, two neighborhood graphs – Delaunay triangulation and a Gabriel graph – were used. They are examples of high and low neighbor counts scenarios, respectively. Most of the observations from the graphs proved not significantly different from the consecutive randomization experiments. Both of the graphs reveal associations between female-sexed and male-sexed bodies. In the case of the Gabriel graph, i.e. the low number of neighbors, bodies of local and non-local origin are associated as well as bodies of mature age category themselves.

The most detailed view of the clustering processes in the cemetery and what kinds of bodies are prevalent in the given clusters is achieved by the percolation analysis. One important observation that becomes clear thanks to the percolation analysis is that there is a bias present in the eastern part of the cemetery. The bodies of undetermined age and localness categorization are present in this part of the cemetery in significantly higher amounts than would be expected if such bodies were distributed at random across the cemetery. In general, the spatial analysis proves that there is a clustering process taking place in the cemetery, and the age, sex, and origin of the buried bodies play structuring roles in the arrangement of the burials. The observed patterns are not too strict, but some are quite constant across the range of applied methods.

At first I have shown that the overall distribution of the artifacts in the graves was not random; the grave goods distributed non-randomly were identified, and subsequently only these were used in the further analysis, effectively removing possible noise created by the artifacts distributed randomly. My main premise here was that the identity of the dead is based mostly on attributes like sex, age and origin. I have successfully demonstrated that especially the sex and age of the body play substantive roles in structuring the assemblage of grave goods deposited within the burials. Some of the relationships between the explored categories and artifacts are well known at the LBK burial grounds, for instance, the association between male burials and polished stone adzes (cf. Bentley et al. 2012; Bickle 2020), and some are known specifically from the Široká u lesa cemetery in Vedrovice, e.g., the bowls being associated with male-sexed bodies (cf. Kvĕtina 2004). Furthermore, the sex, age and origin of the body as layers of identity are, to a certain degree, a structuring factor for the overall spatial organization of the cemetery.

In accordance with the analysis of Robb and Harris (2017), I have found that sex of the body mattered the most in case of structuring the artifact assemblages deposited within the burials at the cemetery of Vedrovice, but its translation into clear-cut identity or (binary) gender categories is not straightforward. In the spatial structure of the cemetery, all the examined categories, sex, age and origin of the body, play roles in the discovered clustering processes, but as well as with the artifacts, the clusters are defined polythetically. Although bodies tend to be buried in clusters, the relationships between sex, age and localness that matter for some of them are not valid for the other clusters. Perhaps there are some yet undiscovered factors like family ties that play substantive roles in structuring the LBK burial grounds spatially or in terms of the grave goods assemblages.


To enable re-use and reproducibility, a compendium containing data and code accompanying this paper is released at Zenodo data repository under DOI The original repository is located at


The research described in this paper was accomplished with support from the project ‘Lifestyle as an unintentional identity in the Neolithic’ (Project 19-16304S), financed by the Czech Science Foundation.

Competing Interests

The author has no competing interests to declare.


  1. Augereau, A. 2022. In search of the origin of inequalities: Gender study and variability of social organization in the first farmers societies of western Europe (Linearbandkeramik culture). Journal of Anthropological Archaeology, 66: 101413. DOI: 

  2. Baddeley, A, Rubak, E and Turner, R. 2015. Spatial point patterns: Methodology and applications with R. London: Chapman and Hall/CRC Press. DOI: 

  3. Bentley, RA, Bickle, P, Fibiger, L, Nowell, G, Dale, CW, Hedges, REM, Hamilton, J, Wahl, J, Francken, M, Grupe, G, Lenneis, E, Teschler-Nicola, M, Arbogast, R-M, Hofmann, D and Whittle, A. 2012. Community differentiation and kinship among Europe’s first farmers. PNAS, 109(24): 9326–9330. DOI: 

  4. Bickle, P. 2020. Thinking Gender Differently: New Approaches to Identity Difference in the Central European Neolithic. Cambridge Archaeological Journal, 30(2): 201–218. DOI: 

  5. Bivand, R and Wong, DWS. 2018. Comparing implementations of global and local indicators of spatial association. TEST, 27(3): 716–748. DOI: 

  6. Bramanti, B. 2008. Ancient DNA: Genetic Analysis of aDNA from Sixteen Skeletons of the Vedrovice Collection. Anthropologie, XLVI(2): 153–160. 

  7. Coppock, A. 2022. Randomizr: Easy-to-use tools for common forms of random assignment and sampling. 

  8. Csardi, G and Nepusz, T. 2006. The igraph software package for complex network research. InterJournal, Complex Systems, 1695. 

  9. Dočkalová, M. 2008. Anthropology of the Neolithic Population from Vedrovice (Czech Republic). Anthropologie, XLVI(2–3): 239–315. 

  10. Gabriel, KR and Sokal, RR. 1969. A New Statistical Approach to Geographic Variation Analysis. Systematic Zoology, 18(3): 259. DOI: 

  11. Gombin, J, Vaidyanathan, R and Agafonkin, V. 2020. Concaveman: A very fast 2D concave hull algorithm. 

  12. Hrnčíř, V, Vondrovský, V and Kvĕtina, P. 2020. Post-marital residence patterns in LBK: Comparison of different models. Journal of Anthropological Archaeology, 59: 101190. DOI: 

  13. Kvĕtina, P. 2004. Mocní muži a sociální identita jednotlivcu˚ – prostorová analýza pohřebištĕ LnK ve Vedrovicích. Archeologické rozhledy, LVI(2): 383–392. 

  14. Maddison, MS and Schmidt, SC. 2020. Percolation Analysis – Archaeological Applications at Widely Different Spatial Scales. Journal of Computer Applications in Archaeology, 3(1): 269–287. DOI: 

  15. Manly, BFJ. 1996. The Statistical Analysis of Artefacts in Graves: Presence and Absence Data. Journal of Archaeological Science, 23(4): 473–484. DOI: 

  16. Masclans Latorre, A, Bickle, P and Hamon, C. 2020. Sexual Inequalities in the Early Neolithic? Exploring Relationships Between Sexes/Genders at the Cemetery of Vedrovice Using Use-Wear Analysis, Diet and Mobility. Journal of Archaeological Method and Theory, 28(1): 232–273. DOI: 

  17. Nakoinz, O and Knitter, D. 2016. Modelling Human Behaviour in Landscapes. Quantitative Archaeology and Archaeological Modelling. Cham: Springer. DOI: 

  18. Negre, J, Muñoz, F and Barceló, JA. 2018. A Cost-Based Ripley’s K Function to Assess Social Strategies in Settlement Patterning. Journal of Archaeological Method and Theory, 25(3): 777–794. DOI: 

  19. Newman, MEJ and Girvan, M. 2004. Finding and evaluating community structure in networks. Physical Review, E 69(026113). DOI: 

  20. Oksanen, J, Blanchet, FG, Friendly, M, Kindt, R, Legendre, P, McGlinn, D, Minchin, PR, O’Hara, R B, Simpson, G L, Solymos, P, Stevens, MHH, Szoecs, E and Wagner, H. 2020. Vegan: Community ecology package. 

  21. Ondruš, V. 2002. Dvĕ pohřebištĕ lidu s neolitickou lineární keramikou ve Vedrovicích. In: Podborský, V (ed.) Dvĕ pohřebištĕ neolitického lidu s lineární keramikou ve Vedrovicích na Moravĕ (Zwei Gräberfelder des Neolithischen Volkes mit Linearbandkeramik in Vedrovice in Mähren). Brno: Masaryk University. pp. 9–149. 

  22. O’Sullivan, D and Unwin, D. 2010. Geographic information analysis. 2nd ed. Hoboken, NJ: John Wiley & Sons. 

  23. Pebesma, E. 2018. Simple Features for R: Standardized Support for Spatial Vector Data. The R Journal, 10(1): 439–446. DOI: 

  24. Perneger, TV. 1998 What’s wrong with Bonferroni adjustments. BMJ, 316:1236 DOI: 

  25. Pettitt, P and Hedges, R. 2008. The Age of the Vedrovice Cemetery: The AMS Radiocarbon Dating Programme. Anthropologie, XLVI(2–3): 125–134. 

  26. Podborský, V. (ed.) 2002. Dvĕ pohřebištĕ neolitického lidu s lineární keramikou ve Vedrovicích na Moravĕ. Brno: Masaryk University. 

  27. R Core Team. 2022. R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing. 

  28. Richards, MP, Montgomery, J, Nehlich, O and Grimes, V. 2008. Isotopic Analysis of Humans and Animals from Vedrovice. Anthropologie, XLVI(2–3): 185–194. 

  29. Ripley, BD. 1976. The second-order analysis of stationary point processes. Journal of Applied Probability, 13(2): 255–266. DOI: 

  30. Robb, J and Harris, OJT. 2017. Becoming gendered in European Prehistory: Was Neolithic gender fundamentally different? American Antiquity, 83(1): 128–147. DOI: 

  31. Rozenfeld, HD, Rybski, D, Andrade, JS, Batty, M, Stanley, HE and Makse, HA. 2008. Laws of population growth. Proceedings of the National Academy of Sciences, 105(48): 18702–18707. DOI: 

  32. Rozenfeld, HD, Rybski, D, Gabaix, X and Makse, HA. 2011. The Area and Population of Cities: New Insights from a Different Perspective on Cities. American Economic Review, 101(5): 2205–2225. DOI: 

  33. Sayer, D and Wienhold, M. 2013. A GIS-Investigation of Four Early Anglo-Saxon Cemeteries: Ripley’s K-function Analysis of Spatial Groupings Amongst Graves. Social Science Computer Review, 31(1): 71–89. DOI: 

  34. Schladitz, K and Baddeley, AJ. 2000. A Third Order Point Process Characteristic. Scandinavian Journal of Statistics, 27(4): 657–671. DOI: 

  35. Schmidt, SC and Maddison, MS. 2020. Percopackage. DOI: 

  36. Sosna, D, Galeta, P and Sládek, V. 2008. A resampling approach to gender relations: The Rebešovice cemetery. Journal of Archaeological Science, 35(2): 342–354. DOI: 

  37. Whittle, A, Bentley, RA, Bickle, P, Dočkalová, M, Fibiger, L, Hamilton, J, Hedges, R, Mateiciucová, I and Pavúk, J. 2013. Moravia and western Slovakia. In: Bickle, P and Whittle, A (eds.), The first farmers of central Europe: Diversity in LBK lifeways. Oxford: Oxbow Books. pp. 101–156. 

  38. Wickham, H, Averick, M, Bryan, J, Chang, W, McGowan, LD, François, R, Grolemund, G, Hayes, A, Henry, L, Hester, J, Kuhn, M, Pedersen, TL, Miller, E, Bache, SM, Müller, K, Ooms, J, Robinson, D, Seidel, DP, Spinu, V, Takahashi, K, Vaughan, D, Wilke, C, Woo, K and Yutani, H. 2019. Welcome to the Tidyverse. Journal of Open Source Software, 4(43): 1686. DOI: