(PDF) Looking at the Big Picture: Using Spatial Statistical Analyses to Study Indigenous Settlement Patterns in the North-Wester

Looking at the Big Picture: Using Spatial Statistical Analyses to Study Indigenous Settlement Patterns in the North- Western Dominican Republic RESEARCH ARTICLE EDUARDO HERRERA MALATESTA JORGE ULLOA HUNG CORINNE L. HOFMAN *Author affiliations can be found in the back matter of this article ABSTRACT CORRESPONDING AUTHOR: Eduardo Herrera Malatesta After several field campaigns between 2007 and 2018 in the northwestern region of Centre for Human Network the Dominican Republic, more than 300 archaeological sites have been registered and Evolution, Aarhus University, revisited. While several of these sites were identified through the scatter of surface Denmark material culture, others show terrain modifications in the form of anthropogenic

mounds and levelled areas. Researchers have gathered valuable information regarding these features’ functionality and construction processes through large-scale excavations in archaeological sites with anthropogenic mounds, paleoenvironmental KEYWORDS: studies and remote sensing analyses. These anthropogenic mounds represent a Caribbean archaeology; long-term process of formation and were used for a variety of purposes ranging from settlement patterns; agricultural to ritual activities. While excavations and small-scale remote sensing can anthropogenic mounds; provide a myriad of data to improve our understanding of these archaeological sites, regional spatial analysis; point a regional perspective is needed to map the relationship among archaeological sites process models with and without terrain modifications, to better understand the Indigenous cultural patterns in their regional environment. In this regard, the primary objective of this TO CITE THIS ARTICLE: paper is to explore to what extent these archaeological sites were related to the Herrera Malatesta, E, Ulloa Hung, J and Hofman, CL. 2023. environment and each other. This was achieved by correlating archaeological data Looking at the Big Picture: with a set of archaeologically recognized important environmental variables using Using Spatial Statistical advanced spatial statistics. The results provide important insights to understand the Analyses to Study Indigenous underlying pattern of archaeological sites in this region and its relationship with the Settlement Patterns in the environmental setting. North-Western Dominican Republic. Journal of Computer Applications in Archaeology, 6(1): 16–28. DOI: https://doi. org/10.5334/jcaa.83 Herrera Malatesta et al. Journal of Computer Applications in Archaeology DOI: 10.5334/jcaa.83 17 1. INTRODUCTION et al., 2020). However, owing to limited resources and time, work has been concentrated on a selection of During several field campaigns between 2007 and 2018 archaeological sites. As a result, most archaeological in the northern region of the Dominican Republic, we were sites with terrain modifications have not been extensively able to register and revisit more than 300 archaeological excavated. Based on the available data and field sites showing a wide range of cultural contexts (Hofman observations, we know that mounds usually occur on the et al., 2018; Ulloa Hung and Herrera Malatesta, 2015; De foothills and hills of the Cordillera Septentrional. Although Ruiter, 2012; Ulloa Hung, 2014; Herrera Malatesta, 2018). there are a few sites on the Yaque river valley that also While several of these sites were identified through contain mounds, these do not present recognizable the scattering of surface material culture, others show levelled areas. As the valley is a floodplain of the Yaque terrain modifications in the form of anthropogenic river, its soils are quite fertile, more so than those of the mounds and levelled areas. Their anthropogenic nature Cordillera Septentrional. This difference between sites was identified by either studying looter pits in already with terrain modifications in the hills and the valley damaged mounds, revealing their stratigraphy, or by brings up the main aim of this paper of exploring the performing test pits and/or small excavations combined relationship between the archaeological sites with and with large-scale excavations on a small number of sites without terrain modifications and the environmental (Ulloa Hung, 2014; Herrera Malatesta, 2018; Hofman setting and assessing if the environmental setting had an et al., 2020; Hofman and Hoogland, 2015). In addition, influence on the location of these terrain modifications. UAV photogrammetry was used to study the intrasite topographic structure of a sample of sites with terrain modifications (Sonnemann, Herrera Malatesta and 2. ARCHAEOLOGICAL BACKGROUND Hofman, 2016; Sonnemann, Ulloa Hung and Hofman, 2016). From the resulting digital terrain models (DTM) The recently available regional archaeological database and the excavations, we observed that artificial mounds has permitted a more comprehensive exploration of how were spatially associated with levelled areas in all Indigenous people inhabited, used, and transformed this cases. From small- and large-scale excavations on the landscape (Hofman et al., 2018). In particular, it allows levelled areas of different sites, we identified patterns the discernment of multiscalar regional patterns. Broadly of postholes. Our current hypothesis is that these areas said, the settlement pattern in the study region, i.e., the were purposely levelled out by the Indigenous people to Puerto Plata, Valverde and Montecristi provinces, shows provide even ground for settlement structures (Hofman a wide range of archaeological sites located in different and Hoogland, 2015). The soil from the areas was swept environmental and topographical settings, with a aside and together with the accumulated dirt from potential tendency towards the coast (Figure 1). Based on domestic and planting activities began to form mounds their archaeological contexts and the number of materials next to the houses. Such mounds could reach up to three distributed on the surface, these sites were grouped into meters in height and twelve meters in diameter, and size categories. The size division was done using Jenks in some cases also served as burial locations (Pagán- Optimization, which is a statistical method that attempts Jiménez et al., 2020; van Dijk, 2019; Sonnemann, Ulloa to find clusters in the data. To do this, it determines the Hung and Hofman, 2016; Hofman and Hoogland, 2015). ‘natural breaks’ in a dataset by minimizing the sum of Furthermore, we were able to perform palaeoecological the squared deviations from class means, in this case, analysis at two of the excavated sites, as well as at key the sites’ area sizes in square meters (Conolly and Lake, places along one of the major rivers, the Yaque, that has 2006, p. 142). Four size categories were defined. First, the provided fundamental information about the human small sites (≤10,000 m2), have a low diversity of material impact on the surrounding environment (Castilla-Beltrán remains, consisting of undecorated ceramics, few stone et al., 2018; Castilla-Beltrán et al., 2020). The recent tools, usually axes and/or hammerstones, and scarce identification of some of the crops that were cultivated mollusk shells. Second, the medium-sized sites (10,001– on the mounds and/or eaten by the Indigenous people 30,000 m2) tend to be more complex, presenting a wider on the site further contextualized this environmental range of material culture, consisting of undecorated and setting (Pagán-Jiménez et al., 2020). These studies decorated ceramics, a variety of stone and shell tools, highlight that the Indigenous people were active agents and more evidence of marine food resources, e.g., various in modifying their surroundings to establish a variety of types of mollusk shells and fish remains. Third, large sites settlements. (30,001–60,000 m2), have characteristics like those of Furthermore, the terrain modifications had a wide medium-sized sites but tend to have more diverse and range of functions for the Indigenous people, from small- abundant material culture and other remains. Finally, the scale agriculture, and dumping of garbage, to ritual very large sites (60,001–140,000 m2) present a material activities (Pagán-Jiménez et al., 2020; van Dijk, 2019; culture pattern like that of large sites, but their spatial Herrera Malatesta, 2018; Hofman et al., 2018; Hofman distribution covers a larger area. Herrera Malatesta et al. Journal of Computer Applications in Archaeology DOI: 10.5334/jcaa.83 18 Figure 1 Archaeological site distribution in the north-western Dominican Republic (map by Eduardo Herrera Malatesta). Archaeological sites with anthropogenic mounds and of brown and dark brown soil, ash lenses and refuse, levelled areas show an arrangement where two or more which alternate with layers of white soil removed from mounds that have either a circular, oval, or “U” shape, the surrounding area to create levelled areas for various are located around the levelled areas. In some cases, the settlement structures (Hofman, 2017). The layers of levelled area is located between mounds or in the middle dark brown soil suggest that during certain periods of of the “U” shape (Figure 2). The registered mounds on occupation, the mounds may have been used as ‘kitchen the foothills and/or hilltops of the Cordillera Septentrional gardens’, i.e. small household gardens (Pagán-Jiménez have diameters ranging from 3 m to 12 m and do not et al., 2020). Additionally, the presence of human and exceed 3 m in height (De Ruiter, 2012; Sonnemann, Ulloa animal burials in some of the mounds also indicates their Hung and Hofman, 2016; Hofman et al., 2020; Herrera possible use as ritual areas (Hofman and Hoogland, 2015; Malatesta, 2018; Hofman and Hoogland, 2015; Ulloa Hofman, Ulloa Hung and Hoogland, 2016; Hofman et al., Hung, 2014). At some sites, the topography has been 2020). Finally, during the large-scale excavations at the greatly affected by human or environmental influences, site of El Flaco, we identified a small number of mounds which complicates the identification of mounds and that were not the result of human intervention but of levelled areas. The levelled areas between the mounds the natural topography of the karstic terrain (Hofman usually have a circular shape and an area between 40 and Hoogland, 2015). Yet, we have not observed this and 200 m2. Their extent also depends on the size of at the (also extensively investigated) site of El Carril, or the settlement: investigated larger sites, featuring more in any other excavated site in the region (Hofman et and bigger mounds tend to also have larger levelled al., 2020; Herrera Malatesta, 2018). The excavations of areas (Hofman et al., 2020; Sonnemann, Ulloa Hung these sites increased our comprehension of their general and Hofman, 2016). Several of the excavated sites, e.g., chronology, covering between 800 and 1500 CE, and El Flaco, El Carril, and El Manantial, featured levelled provided a better context for the ceramic repertoire of areas, or areas in between mounds, with postholes and the three main ceramic groups present in the region, other cultural features, which confirmed the existence of i.e. Ostionoid (600–1200 CE), Meillacoid (800–1550 settlement structures (Figure 3), possibly the houses of CE), and Chicoid (900–1700 CE) (Keegan and Hofman, the Indigenous inhabitants (Hofman et al., 2018; Herrera 2017). This new evidence on the relationship between Malatesta, 2018). The mounds are composed of layers sites with terrain modifications and the environment, Herrera Malatesta et al. Journal of Computer Applications in Archaeology DOI: 10.5334/jcaa.83 19 Figure 2 Three topographically modified settlement sites from the survey region. The topography was recorded from drone imagery and a 3D model produced by photogrammetry. Each site is visualized by DEM (50% visibility), over slope (50% visibility), and over a hillshade visualization. Based on a visual analysis the mounds and levelled areas were drawn by hand (models created by Till Sonnemann). Figure 3 Stratigraphy of mound units at four sites with alternating layers and lenses of ash and soil. A) El Carril, B) El Flaco, C) El Manantial (photographs A & B by Menno Hoogland, C by Eduardo Herrera Malatesta). Herrera Malatesta et al. Journal of Computer Applications in Archaeology DOI: 10.5334/jcaa.83 20 including their topographic location strongly highlights environmental variables that have already been identified the need to study and understand the regional pattern as having some sort of relationship with the archaeological of archaeological sites with terrain modifications. distribution of sites with and without terrain modifications. The selected variables are related to the geomorphology and soil quality of the region. To facilitate locating them 3. MATERIAL AND METHODS in the Atlas, we decided to use an English translation of their original Spanish names. From the geomorphology, 3.1. ARCHAEOLOGICAL DATA the covariates are (i) alluvium, (ii) hills and plateaus, and The focal region encompasses geopolitically the modern (iii) mountainous areas. The alluvium variable was chosen provinces of Montecristi, Puerto Plata, and Valverde in the because it represents an area in the valley where some Dominican Republic and topographically the northern important sites with terrain modification have been coastal strip, Cordillera Septentrional and the Cibao valley, registered. The area itself is close to the major rivers in which includes the Yaque river (Figure 1). The database the study region, and the land on and near the alluvium is results from a compilation of regional archaeological highly arable and fertile soil. The second and third variables surveys carried out in this region since 2007 (Ulloa Hung, were selected because several sites with and without 2014; Ulloa Hung and Herrera Malatesta, 2015; Herrera terrain modifications are located in the low-relief hills of Malatesta, 2018; Herrera Malatesta and Hofman, 2019). the Cordillera Septentrional and the more mountainous Several factors, from environmental (dense thorny areas. From the quality of soils, the covariates are: (i) bushland, marshes) to topographical (steep valleys and arable soils and (ii) soils suitable for forests, pastures and mountains) hindered a systematic survey approach for mountain crops. These two variables were selected to test the great extent of the study region. Therefore, sites if the location of the archaeological sites with and without were primarily recorded using an opportunistic survey terrain modifications are being affected by the presence methodology (Orton, 2000), which involved intense of arable soils in the valley and the low relief hills and/or ethnographical work with local communities and the suitable areas for crops located in the mountains. To collaborations with local guides to locate archaeological prepare these covariates for the spatial statistical analysis, sites. Alternatively, some sites were registered by a Euclidian distance map was calculated for each variable implementing predictive modelling methods that using ESRI ArcGIS 10.7. Since the vectors created from determined areas with a high probability of archaeological the environmental variables are based on copyrighted evidence (Herrera Malatesta, 2017). For the present study, materials, though allowed to be replicated and used for the 285 of the more than 300 analyzed archaeological sites purposes of the research study, it is not possible to share were used, the remaining ones being unsuitable for a the digital data. Nonetheless, all data can be accessed in variety of reasons, e.g., the total or partial alteration of paper form in the mentioned publications. the site; dense vegetation that prevented full recording of evidence on the ground, or sites in caves or rock shelters. 3.3. METHODS Of the 285, less than 5% were previously registered by To study the spatial relationships of the Indigenous Dominican researchers (Ortega, 2005; Veloz Maggiolo settlement patterns, we focused on applying Point and Ortega, 1980). The database contains information Pattern Analysis, which is simply defined as the study regarding the location of sites, their associated material of the spatial arrangement of points in space (Riris, culture, functionality, the size of the site, and the presence/ 2020; Costanzo et al., 2021; Carrero-Pazos, Bevan and absence of mounds. Since the database is composed of Lake, 2019). When the intensity of the point pattern is surface finds, its structure is mostly in a binary format of constant, i.e., corresponding to a stationary and isotropic presence/absence of material evidence. For the spatial process within the region, it is called a Homogeneous statistical analysis, the dataset of archaeological sites was Poisson Process (HPP). Alternatively, when the intensity divided into two groups. One contains the archaeological is not spatially uniform it is called Inhomogeneous sites with terrain modifications (n = 58) and the other one, Poisson Process (IPP). If evidence indicates that a the archaeological sites without terrain modifications (n regional pattern can be the result of an IPP, it is relevant = 227). The dataset used for this paper is open access to assess whether the intensity of the events depends (Herrera Malatesta and Hofman, 2019). on the influence of spatial variables (the covariates) and, if so, it becomes necessary to quantify this dependence. 3.2. ENVIRONMENTAL DATA In spatial statistics, this is referred to as the intensity of The related environmental data consist of modern records the pattern depending on First-order effects. However, compiled in two atlases published by the Environmental when the intensity of the point pattern is the result of Ministry of the Dominican Republic (Moya Pons, 2004; the points showing independence or rejection from each MMARN, 2012). Considering the available archaeological other, then the point pattern is the result of Second- and palaeobotanical data from the sites with terrain order effects (Baddeley, Rubak and Turner, 2016; Diggle, modifications, the focus was on exploring a group of 2013). Based on this framework and the archaeological Herrera Malatesta et al. Journal of Computer Applications in Archaeology DOI: 10.5334/jcaa.83 21 background of this region, this paper aims to determine After fitting the covariates and the archaeological the impact of the regional environmental setting site distribution, the resulting models were tested by on the distribution of archaeological sites with and a goodness-of-fit for each of the six models (‘null’, without terrain modifications. For this, we defined two ‘first-’, and ‘second-order’ for the sites with and without complimentary null hypotheses: terrain modifications), comparatively employing a pair correlation function (PCF). This method provides H1 the intensity of the point pattern of an estimation of the degree of spatial clustering and archaeological sites with and without terrain regularity of the spatial pattern at multiples scales within modifications is uniform, i.e., it is stationary and a window of observations (Baddeley, Rubak and Turner, isotropic. The environmental variables (covariates) 2016). Other common methods to test second-order do not affect the spatial distribution of the points. effects are the Ripley’s K and L functions (Ripley, 1977; H2 the archaeological sites with and without Orton, 2004; Bevan and Conolly, 2006; Bivand, Pebesma terrain modifications are completely independent and Gómez-Rubio, 2008; Vanzetti et al., 2010; Bevan et of each other, and they do not show any spatial al., 2013), which analyses point patterns at different scale correlation. ranges. However, considering the characteristics of the available data and the hypothesis to be tested, PCF was a more adequate method. This is because, contrary to the K 3.4. FITTING COVARIATE SPATIAL DATA and L functions, PCF measures the intensity of pair points To statistically address these two hypotheses, the point in doughnut-shaped rings from and around each point process intensity was explored by applying parametric (Baddeley, Rubak and Turner, 2016) in a non-cumulative methods that explore the degree to which the covariates matter, allowing an improved measure for spatial affect the distribution of archaeological sites. The first association at different scales. PCF quantifies how each step to fit the covariates was to test the presence of point in the dataset is surrounded by other points, and multicollinearity using Pearson’s R, to assess if each therefore the resulting value will present the variation of variable can be used as an independent model for further site density across the study region. This allows a better analysis. This test showed that there are no collinear differentiation of patterns at different scales. effects among the five considered variables. Then, the To illustrate the critical scales where sites with and observed point pattern of archaeological sites with and without terrain modifications attract or reject one another, without terrain modifications was fitted to a Poisson point a significance envelope was computed via Monte Carlo process model, equivalent to an initial null hypothesis of simulations from the ‘null’, ‘first-’, and ‘second-order’ complete spatial randomness (CSR). For this, the Bayesian models. The envelopes were generated by simulating Information Criterion (BIC) was used as a guideline 999 point patterns for each calculation of the PCF to the for information gain per model. The stepwise model six fitted models (Baddeley and Turner, 2005). Rather selection procedure rejected the same covariate for the than expressing confidence intervals at a level of 95%, two groups of sites, i.e., soils suitable for forests, pastures these envelopes estimated acceptance intervals (or non- and mountain crops. The resulting four covariates for each rejection intervals), which are the likelihood of incorrectly group of archaeological sites, and the assumed Poisson rejecting the null hypothesis of no spatial structure process, comprise the multivariate first-order model. as a function of the number of simulated runs. For the The spatial intensity (density) of the archaeological sites envelopes generated here, the value of α corresponded to with and without terrain modifications in this region as 0.002, which is highly robust. The workflow for fitting and described by the interaction of these parameters, is a testing the models was fully implemented as R code, and function of their combined maximum logistic likelihood it is available in the supplementary information. It was (Baddeley and Turner, 2000). Next, to assess the possible based on the codes of Bevan et al. (2013) and Riris (2020), effects of interaction between the sites with and without and it draws on functions from several R packages, such terrain modifications at a landscape level, a second-order as ‘sp’, ‘rgdal’, ‘raster’, ‘MASS’, ‘maptools’, ‘dismo’, ‘DAAG’ test was incorporated into the fitted first-order model and ‘spatstat’ (Pebesma and Bivand, 2005; Maindonald to generate a combined model. This model statistically and Braun, 2010; Baddeley and Turner, 2005; Bivand, accounted for induced and internal spatial dynamics 2022; Bivand et al., 2022; Ripley et al., 2022; Hijmans, simultaneously. For this, a Strauss point process model 2021; Hijmans et al., 2022). was advanced to replace CSR as the null model and to obtain the range of observed variability in the empirical point pattern (Kelly and Ripley, 1976). Although there 4. SPATIAL POINT PROCESS are other alternatives such as Poisson, Area Interaction MODELLING or Multiscale Geyer Process (Baddeley and Turner, 2005; Riris, 2020), the selected function provided clearer results 4.1. NULL MODEL to understand the distribution of sites with and without All spatial statistical processes begin with the assumption terrain modifications in this region. that point patterns are distributed in space stochastically Herrera Malatesta et al. Journal of Computer Applications in Archaeology DOI: 10.5334/jcaa.83 22 (Illian et al., 2008), and this random pattern is called terrain modifications from 4000 m up to 10000 m. Yet, Complete Spatial Randomness (CSR). When a point pattern the results for the two groups of sites are not conclusive is defined by a stochastic process, it is said to have a and are just a baseline for further explorations. This Poisson process that affects the distribution, i.e., that the baseline assessment lacks sufficient context to properly events occur continuously and independently from each explore the implications of these results, and therefore other. For the two sets of archaeological sites used in this the analysis turns to the assessment of the fitted study, their spatial structure was illustrated by the initial covariate data and point process models. application of the PCF (Figures 4 and 5, top left). Figure 4 (top left) shows that the archaeological distribution 4.2. FIRST-ORDER MODEL of sites with terrain modifications shows a pattern of The first-order model (Figures 4 and 5, bottom left) repulsion at all spatial scales. The red dotted line (equal presents a more interesting fit over the CSR null model. to 1) marks the value of CSR, all values above this line For the sites with terrain modifications, this model represent a cluster point pattern and all below a regular (Figure 4, bottom left) controls all the signals of spatial point pattern. The grey area indicates the Monte Carlo autocorrelation at all ranges within the spatial frame. Of simulation threshold (Baddeley, Rubak and Turner, 2016). the five tested covariates, four showed highly significant For the sites without terrain modifications (Figure correlations (Table 1). Interestingly, when the covariates 5, top left) the first maximum and minimum values are included, the pattern shows two clusters of sites, outside the simulation envelope mark the range at which one in the range of 3800 m and 6000 m and another points are most strongly attracted or repelled by each at 7800 m and 9500 m, approximately. Before the first other. For the study region, there is a peak cluster of sites cluster, the pattern shows CSR; after the first and second without terrain modifications at around 500 m distance clusters, there is repulsion among the sites. The intensity from each other that decreases over space to reach a surface (log-odds of points per unit of area) from this maximum point of rejection, i.e., a regular pattern, at model (Figure 4, top right) offers a visual summary around 4000 m. This initial pattern implies that a variety of the degree to which this model explains the spatial of spatial processes may be affecting the distribution structure of sites with terrain modifications. The pattern suggesting a strong range clustering up to 2000 m and is inhomogeneous, this is spatially variable along the a weaker longer-range repulsion between sites without study region, and the predicted spatial trend is strong, Figure 4 Results of point process modelling for sites with terrain modifications. Null: The pair correlation function estimated on an assumption of complete spatial randomness. Intensity surface: the predicted intensity surface created from the first-order fit. First- order: the pair correlation function with a critical envelope conditioned on the covariate data as first-order variables. Second-order: incorporating a point interaction term in the first-order model, accounting for all spatial variability. Herrera Malatesta et al. Journal of Computer Applications in Archaeology DOI: 10.5334/jcaa.83 23 Figure 5 Results of point process modelling for sites without terrain modifications. Null: The pair correlation function estimated on an assumption of complete spatial randomness. Intensity surface: the predicted intensity surface created from the first-order fit. First- order: the pair correlation function with a critical envelope conditioned on the covariate data as first-order variables. Second-order: incorporating a point interaction term in the first-order model, accounting for all spatial variability. COVARIATES ESTIMATE STD. ERROR Z VALUE SIGNIFICANT (Intercept) –1.7504450 4.7434160 –36.902629 *** Alluvium 0,0001872681 0,00003771534 4.965303 *** Hills and Plateaus –0,0002048970 0,00004904832 –4.177451 *** Mountainous areas –0,0001380711 0,00003149600 –4.383766 *** Arable soils –0,0001164922 0,00003183083 –3.659729 *** Table 1 Fitted covariate datasets for the first-order model, sites with terrain modifications. COVARIATES ESTIMATE STD. ERROR Z VALUE SIGNIFICANT (Intercept) –1.6048930 2.0523640 –78.1972737 *** Alluvium 0,0001314754 0,00001791114 7.3404269 *** Hills and Plateaus –0,0001393845 0,00002166846 –6.4325969 *** Mountainous areas –0,00006091828 0,00001397768 –4.3582533 *** Arable soils –0,00001786111 0,00003669478 –0.4867479 *** Table 2 Fitted covariate datasets for the first-order model, sites without terrain modifications. with several zones highlighting values at the high end of correlations (Table 2). The covariate presence provides an the scale. Together with the results of the pair correlation improved fit compared to the CSR null model, controlling functions, the fitted first-order model suggests that there most of the signals of spatial autocorrelation. In this is spatial variation within this pattern. iteration, the strong pattern of short-range clustering For the sites without terrain modifications, of the and long-range repulsion remains. However, a small five tested covariates, four showed highly significant new cluster appears at a range of 4800 m and 5800 m, Herrera Malatesta et al. Journal of Computer Applications in Archaeology DOI: 10.5334/jcaa.83 24 approximately. As with the previous group of sites, the results are satisfactory as they present an explanation intensity surface (Figure 5, top right) provides a clear view of the landscape-level archaeological distribution of of the degree to which this model explains the spatial the sites with terrain modifications. Among the results structure of sites without terrain modifications. Again, of the two groups of sites, a noticeable exception not the inhomogeneous pattern shows spatial variability in accounted for by the second-order model is a short- the study region, and the strong predicted spatial trend distance statistically significant clustering of sites illustrates several zones with values at the high end of without terrain modifications in a range of 500 m and the scale. Together with the results of the pair correlation 800 m. This indicates a relevant spatial threshold for the functions, the fitted first-order model suggests that there sites without terrain modifications, which could help with is spatial variation within this pattern. the understanding of the distributional pattern of this type of site within the region. 4.3. SECOND-ORDER MODEL The second-order model was created by fitting a nonstationary Strauss process by maximizing the profile 5. DISCUSSION information entropy with respect to a series of fixed spatial interaction distances (Baddeley, Rubak and This paper has applied advanced spatial statistical Turner, 2016). For the sites with terrain modifications, methods to study a regional archaeological dataset in the analysis was permuted in 100-meter intervals from the northwestern Dominican Republic to augment the 100 m to 2500 m to produce 25 candidate models. understanding of past Indigenous regional settlement According to Akaike’s Information Criterion (AIC), this patterns. It demonstrated the analytical possibilities procedure provides the best fit for the area interaction of using spatial statistics for non-systematic regional parameter at r = 2500 m. In this case, the AIC was used survey data. As verified in other regions, applying spatial instead of the BIC, as while BIC is better at selecting the statistical methods in correlation with environmental true model than AIC, the latter is optimal for estimating and cultural variables proved to be a valuable tool to the regression function (Yang, 2005). The second model understand regional patterns (Riris, 2020; Bevan et al., also incorporates the first-order model as a trend 2013), complementing previous site-based research in the parameter. For the sites without terrain modifications, region. The performed analysis proves the importance of the analysis was permuted in 250-metre intervals from comparing archaeological data to known environmental 250 m to 6500 m to produce 25 candidate models. In variables to explore their correlation on a large scale. The this case, Akaike’s Information Criterion indicated that analyses have also proven that further tests must be the best fit for this group of sites in the area interaction done to continue looking for explanations on the cultural parameter is at r = 750 m. This combined first- and reasons behind the location of settlements and the second-order model successfully controls for both construction of mounds and levelled areas in the region. exogenous and endogenous processes within the point Nonetheless, the methods presented here provided a pattern of archaeological sites with and without terrain robust and rigorous starting point for future regional and modifications. Conditioning the simulation envelope for landscape research in the region. the pair correlation function on the combined null model The results provided two main takeaways. First, both (Figures 4 and 5, bottom right) illustrates the effect of sets of archaeological sites presented the same positive incorporating between-point interaction. and negative correlations with the same variables (Tables For the sites with terrain modifications (Figure 4, bottom 1 and 2). In archaeological terms, this could mean that, on right), the results indicate that most of the variability in the one hand, these variables do not explain the variability the spatial structure can be explained by controlling the of archaeological sites as expected. On the other hand, it second-order factors. Like with the first-order trend, the could mean that, contrary to expectations, in terms of the second-order trend affects the site patterning, as both regional pattern there is no difference between the sites models show clustering at 3800 m to 6000 m and 7800 with and without terrain modifications and the presence m to 9800 m, approximately. For the sites without terrain or absence of these environmental variables. This means modifications, a more interesting picture is presented by that the potential reasoning behind the location of the the results of the second-order trend (Figure 5, bottom sites and the eventual construction of levelled areas and right). The second-order model emphasizes the presence mounds was similar and not led by certain environmental of the second and third clusters, which are smaller in variables. The results indicated that both sets of sites have the first-order trend and almost non-existent with the a positive correlation with the alluvium variable, which null model. This represents a greater influence of the is physically located in the Cibao valley, running along second-order trend within this distribution. While other the Yaque river. This implies that, contrary to the visual analyses with a better goodness-of-fit could be done to assessment that the sites are clustered towards the coast, further untangle the relationship between the influences the statistical results suggest that the regional pattern of first and second-order effects in the site pattern, the tends to be towards the valley, the major water sources, Herrera Malatesta et al. Journal of Computer Applications in Archaeology DOI: 10.5334/jcaa.83 25 and the high-productivity soils. Along the same line, the dans-xyn-cu72. The rest of the data will be also deposited negative correlation with the other three variables does in DANS later this year. Supplementary materials of not negate the presence of archaeological sites in the this article can be found online on the Open Science same spatial areas, just that the regional pattern shows a Framework at https://doi.org/10.17605/OSF.IO/WUBA7. tendency of rejection towards those variables. Secondly, the analysis provided valuable information regarding the structures of the regional patterns ACKNOWLEDGEMENTS concerning the three models. For the two sets of archaeological sites, there were clear differences between Many thanks go out to the communities in the northern the null model and the first and second-order trends. Dominican Republic for their collaboration during our For the sites with terrain modifications, the null model research. The Ministerio de Cultura and the Museo del presented an image of a regular regional pattern with no Hombre Dominicano are acknowledged for providing insights into the decision of locating sites. For the sites us with permissions and research authorization in the without terrain modifications, the null model showed a country, as well as the Servicio Geológico Nacional and significant cluster at distances <2000 m, and significant the Ministerio de Ambiente for providing environmental regular patterns at 4000, 6–7000, and above 9000 m, and geopolitical digital data of the country. approximately. Yet, when the environmental variables were considered, patterns of clustered sites appeared at different spatial scales for both sets of archaeological FUNDING INFORMATION sites. This indicates that the regional site distribution could be influenced by the environmental setting. This research has received funding from the European Furthermore, when the pattern was fitted against the Research Council under the European Union’s Seventh second-order trend, the clusters were even more clear. Framework Programme (FP7/2007–2013)/ERC grant These results indicate that, as Riris (2020) proposed, agreement n° 319209. Supplementary funds for the regional spatial palimpsests cannot be explained by realization and publishing of this paper came from fixed-scale summary statistics under a null assumption Horizon Europe, HORIZON-MSCA-2021-PF-01-01, Marie of spatial homogeneity. They also show that the regional Skłodowska Curie Action, Grant agreement n° 101062882. archaeological pattern is affected by exogenous factors that influence the locations of points at different spatial scales, yet the cluster pattern becomes clearer and more COMPETING INTERESTS significant when the intrinsic properties of the point pattern are considered. With these results, the two null The authors have no competing interests to declare. hypotheses are rejected, as there was evidence that the environmental variables (covariates) do affect the spatial distribution of the points, and the archaeological sites AUTHOR CONTRIBUTIONS with and without terrain modifications do show a spatial correlation. EHM conceptualized and designed the study, wrote the Past Indigenous peoples took advantage of the main text, executed the analyses, and created most topographical conditions in and around their settlements figures. JUH and CLH contribute to the paper’s ideas, to remodel the landscape to their needs, and the terrain wrote specific sections, contribute to the figures, and modifications are evidence of this. They transformed provided unpublished data. the surroundings with new topographies. The results emphasize the complexity and multiscalar conditions between the archaeological sites and environmental AUTHOR AFFILIATIONS context. This gives way to a debate on how past Indigenous peoples interacted with their environment Eduardo Herrera Malatesta orcid.org/0000-0001-5265-6296 and created this complex landscape, where specific Centre for Human Network Evolution, Aarhus University, cultural features correlate with environmental conditions Denmark at different spatial scales. Jorge Ulloa Hung orcid.org/0000-0002-2680-1530 Instituto Tecnologico de Santo Domingo and Anthropology Deparment University of Miami, US DATA ACCESSIBILITY STATEMENT Corinne L. Hofman orcid.org/0000-0003-4447-5019 Leiden University, Faculty of Archaeology, The Netherlands; The The Montecristi data used in this research is open access Royal Netherlands Institute of Southeast Asian and Caribbean at the DANS EASY repository https://doi.org/10.17026/ Studies, The Netherlands Herrera Malatesta et al. Journal of Computer Applications in Archaeology DOI: 10.5334/jcaa.83 26 REFERENCES Castilla-Beltrán, A, Hooghiemstra, H, Hoogland, MLP, Pagán- Jiménez, J, van Geel, B, Field, MH, Prins, M, Donders, Baddeley, A, Rubak, E and Turner, R. 2016. Spatial point T, Herrera Malatesta, E, Ulloa Hung, J, McMichael, patterns: methodology and applications with R. Boca Raton: CH, Gosling, WD and Hofman, CL. 2018. ‘Columbus’ Chapman and Hall/CRC. DOI: https://doi.org/10.1201/ footprint in Hispaniola: A paleoenvironmental record of b19708 indigenous and colonial impacts on the landscape of Baddeley, A and Turner, R. 2000. ‘Practical Maximum the central Cibao Valley, northern Dominican Republic’. Pseudolikelihood for Spatial Point Patterns’. Australian & Anthropocene, 22: 66–80. DOI: https://doi.org/10.1016/j. New Zealand Journal of Statistics, 42(3): 283–322. DOI: ancene.2018.05.003 https://doi.org/10.1111/1467-842X.00128 Conolly, J and Lake, M. 2006. Geographical information systems Baddeley, A and Turner, R. 2005. ‘spatstat: An R Package for in archaeology. Cambridge: Cambridge University Press. Analyzing Spatial Point Patterns’. Journal of Statistical DOI: https://doi.org/10.1017/CBO9780511807459 Software, 12(6): 1–42. DOI: https://doi.org/10.18637/jss. Costanzo, S, Brandolini, F, Idriss Ahmed, H, Zerboni, A and v012.i06 Manzo, A. 2021. ‘Creating the funerary landscape of Bevan, A and Conolly, J. 2006. ‘Multiscalar approaches to Eastern Sudan’. PLOS ONE, 16(7): e0253511. DOI: https:// settlement pattern analysis’. In Lock, G and Molyneaux, BL doi.org/10.1371/journal.pone.0253511 (eds.), Confronting scale in archaeology: Issues of Theory De Ruiter, S. 2012. Mapping History: An analysis of site locations and Practice, 217–234. New York: Springer. in the northwestern Dominican Republic. RMA thesis, Leiden Bevan, A, Crema, E, Li, X and Palmisano, A. 2013. ‘Intensities, University, Leiden. interactions and uncertainties: some new approaches Diggle, PJ. 2013. Statistical Analysis of Spatial and Spatio- to archaeological distributions’. In Bevan, A and Lake, M Temporal Point Patterns. 3rd edition. Boca Raton: Chapman (eds.), Computational Approaches to Archaeological Spaces, & Hall/CRC. DOI: https://doi.org/10.1201/b15326 27–52. Walnut Creek: Left Coast Press. Dijk, K. van. 2019. Mounded Landscapes. The distribution of past Bivand, RS. 2022. ‘Combining Spatial Data’. Online vignette human activities associated with pre-colonial mounds at El originally from 120–126 of the first edition of Bivand, R. S., Carril, Dominican Republic. MA thesis, Leiden University. Pebesma, E. and Gómez-Rubio V. (2008) Applied Spatial Herrera Malatesta, E. 2017. ‘Understanding ancient patterns: Data Analysis with R, Springer-Verlag, New York. It was Predictive Modeling for field research in Northern retired from the second edition (2013) to accommodate Dominican Republic’. International Association of Caribbean material on other topics, and is made available in this form Archaeology (IACA), Sint Maarten: SIMARC Heritage Series, with the understanding of the publishers. Published online 88–97. by the Cran R Project at https://cran.r-project.org/web/ Herrera Malatesta, E. 2018. Una isla, dos mundos. Estudio packages/maptools/vignettes/combine_maptools.pdf. arqueológico sobre el paisaje indígena de Haytí y su Bivand, RS, Lewin-Koh, N, Pebesma, E, Archer, E, Baddeley, transformación al paisaje colonial de La Española (1200– A, Bearman, N, Bibiko, HJ, Brey, S, Callahan, J, Carrillo, 1550). Leiden: Sidestone Press. G, Dray, S, Forrest, D, Friendly, M, Giraudoux, P, Golicher, Herrera Malatesta, E and Hofman, CL. 2019. ‘Indigenous D, Gómez Rubio, V, Hausmann, P, Hufthammer, KO, landscape transformation on northern Haytí: An Jagger, T, Johnson, K, Lewis, M, Luque, S, MacQueen, D, archaeological and environmental database of the Niccolai, A, Perpiñán Lamigueiro, O, Plunkett, E, Rubak, Montecristi coast’. Journal of Open Archaeology Data, 7: E, Short, T, Snow, G, Stabler, B, Stokely, M and Turner, 1–5. DOI: https://doi.org/10.5334/joad.51 R. 2022. ‘Package ‘maptools’’. Published online by the Hijmans, RJ. 2021. Introduction to R. Published online by the Cran R Project at https://cran.r-project.org/web/packages/ Cran R Project at https://cran.r-project.org/doc/manuals/r- maptools/index.html. release/R-intro.pdf. Bivand, RS, Pebesma, EJ and Gómez-Rubio, V. 2008. Applied Hijmans, RJ, Phillips, S, Leathwick, J and Elith, J. 2022. spatial data analysis with R. New York: Springer. ‘Package ‘dismo’’. https://cran.r-project.org/web/packages/ Carrero-Pazos, M, Bevan, A and Lake, MW. 2019. ‘The spatial dismo/index.html. structure of Galician megalithic landscapes (NW iberia): Hofman, CL. 2017. “Nexus 1492 Fieldwork Report Dominican A case study from the Monte Penide region’. Journal of Republic, 2017”. https://doi.org/10.17026/dans-29r-9qtr, Archaeological Science, 108: 104968. DOI: https://doi. DANS Data Station Archaeology, V2. org/10.1016/j.jas.2019.05.004 Hofman, CL and Hoogland, MLP. 2015. ‘Investigaciones Castilla-Beltrán, A, Hooghiemstra, H, Hoogland, MLP, arqueológicas en los sitios El Flaco (Loma de Guayacanes) Donders, TH, Pagán-Jiménez, JR, McMichael, CNH, y La Luperona (Unijica). Informe pre-liminar’. Boletin del Rolefes, SMF, Olijhoek, T, Herrera Malatesta, E, Hung, Museo del Hombre Dominicano, 46: 61–74. JU and Hofman, CL. 2020. ‘Ecological responses to land Hofman, CL, Hoogland, MLP and Pagán-Jiménez, J. 2020. ‘En use change in the face of European colonization of Haytí la víspera de la colonización europea: Los sitios indígenas island’. Quaternary Science Reviews, 241: 106407. DOI: El Flaco, El Carril y El Cabo, República Dominicana’. Boletín https://doi.org/10.1016/j.quascirev.2020.106407 del Museo del Hombre Dominicano, XLVII(48): 133–152. Herrera Malatesta et al. Journal of Computer Applications in Archaeology DOI: 10.5334/jcaa.83 27 Hofman, CL, Ulloa Hung, J, Herrera Malatesta, E, Jean, JS, northern Dominican Republic’. Review of Palaeobotany and Sonnemann, TF and Hoogland, MLP. 2018. ‘Indigenous Palynology, 274: 104160. DOI: https://doi.org/10.1016/j. Caribbean perspectives: archaeologies and legacies of the revpalbo.2020.104160 first colonised region in the New World’. Antiquity, 92(361): Pebesma, E and Bivand, RS. 2005. ‘Classes and Methods for 200–216. DOI: https://doi.org/10.15184/aqy.2017.247 Spatial Data: the sp Package’. Published online by the Cran Hofman, CL, Ulloa Hung, J and Hoogland, MLP. 2016. ‘El R Project at https://cran.r-project.org/package=sp/sp.pdf. paisaje social indígena al momento del encuentro colonial: Ripley, BD. 1977. ‘Modelling spatial patterns’. Journal of the nuevas investigaciones en el norte de la República Royal Statistical Society B, 39: 172–212. DOI: https://doi. Dominicana’. Boletín del Museo del Hombre Dominicano, org/10.1111/j.2517-6161.1977.tb01615.x 47: 300–310. Ripley, BD, Venables, B, Bates, DM, Hornik, K, Gebhardt, Illian, J, Penttinen, A, Stoyan, H and Stoyan, D. 2008. A and Firth, D. 2022. ‘Package ‘maptools’’. Published Statistical analysis and modelling of spatial point online by the Cran R Project at https://cran.r-project.org/ patterns. New York: Wiley-Interscience. DOI: https://doi. package=maptools. org/10.1002/9780470725160 Riris, P. 2020. ‘Spatial structure among the geometric Keegan, WF and Hofman, CL. 2017. The Caribbean before earthworks of western Amazonia (Acre, Brazil)’. Journal of Columbus. New York: Oxford University Press. DOI: https:// Anthropological Archaeology, 59: 101177. DOI: https://doi. doi.org/10.1093/acprof:oso/9780190605247.001.0001 org/10.1016/j.jaa.2020.101177 Kelly, FP and Ripley, BD. 1976. ‘A Note on Strauss’s Model for Sonnemann, TF, Herrera Malatesta, E and Hofman, CL. 2016. Clustering’. Biometrika, 63(2): 357–360. DOI: https://doi. ‘Applying UAS photogrammetry to analyze spatial patterns org/10.1093/biomet/63.2.357 of indigenous settlement sites in the northern Dominican Maindonald, J and Braun, WJ. 2010. Data analysis and graphics Republic.’ In Forte, M and Campana, S (eds.), Digital using R: An example-based approach. Cambridge Series Methods and Remote Sensing in Archaeology. Archaeology in Statistical and Probabilistic Mathematics (3rd Edn.). in the Age of Sensing, 71–87. Switzerland: Springer. DOI: Cambridge: Cambridge University Press. https://doi.org/10.1007/978-3-319-40658-9_4 MMARN. 2012. Atlas de biodiversidad y recursos naturales Sonnemann, TF, Ulloa Hung, J and Hofman, CL. 2016. de la República Dominicana. Santo Domingo, República ‘Mapping indigenous settlement topography in the Dominicana: Ministerio de Medio Ambiente y Recursos Caribbean using drones’. Remote Sensing, 8(791). DOI: Naturales. https://doi.org/10.3390/rs8100791 Moya Pons, F. 2004. Atlas de los recursos naturales de la Ulloa Hung, J. 2014. Arqueología en la línea noroeste de La República Dominicana. Santo Domingo, República Española. Paisajes, cerámicas e interacciones. Santo Dominicana: Secretaría de Estado de Medio Ambiente y Domingo: Instituto Tecnológico de Santo Domingo. Recursos Naturales. Ulloa Hung, J and Herrera Malatesta, E. 2015. ‘Investigaciones Ortega, E. 2005. Compendio general arqueológico de Santo arqueológicas en el Norte de la Española, entre viejos Domingo. Santo Domingo: Academia de Ciencias de la esquemas y nuevos datos’. Boletín del Museo del Hombre República Dominicana. Dominicano, 46(XLII): 75–107. Orton, C. 2000. Sampling in archaeology. Cambridge: Vanzetti, A, Vidale, M, Gallinaro, M, Frayer, DW and Bondioli, Cambridge University Press. DOI: https://doi.org/10.1017/ L. 2010. ‘The iceman as a burial’. Antiquity, 84(325): 681– CBO9781139163996 692. DOI: https://doi.org/10.1017/S0003598X0010016X Orton, C. 2004. ‘Point pattern analysis revisited’. Archeologia Veloz Maggiolo, M and Ortega, E. 1980. ‘Nuevos hallazgos e Calcolatori, 15: 299–315. http://eprints.bice.rm.cnr.it/id/ arqueológicos en la costa norte de Santo Domingo’. eprint/920. Boletín del Museo del Hombre Dominicano, 13: 11–48. Pagán-Jiménez, JR, Ali, Z, Santiago-Marrero, CG and Hofman, Yang, Y. 2005. ‘Can the strengths of AIC and BIC be shared? CL. 2020. ‘Plantscapes of dwelling: Precolonial household A conflict between model indentification and regression mounds, phytocultural dynamics and the ensuing human estimation’. Biometrika, 92(4): 937–950. https://www.jstor. ecosystems at El Flaco and El Carril (cal. AD 990–1450), org/stable/20441246. Herrera Malatesta et al. Journal of Computer Applications in Archaeology DOI: 10.5334/jcaa.83 28 TO CITE THIS ARTICLE: Herrera Malatesta, E, Ulloa Hung, J and Hofman, CL. 2023. Looking at the Big Picture: Using Spatial Statistical Analyses to Study Indigenous Settlement Patterns in the North-Western Dominican Republic. Journal of Computer Applications in Archaeology, 6(1): 16–28. DOI: https://doi.org/10.5334/jcaa.83 Submitted: 13 August 2021 Accepted: 02 March 2023 Published: 23 March 2023 COPYRIGHT: © 2023 The Author(s). This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International License (CC-BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. See http://creativecommons.org/licenses/by/4.0/. Journal of Computer Applications in Archaeology is a peer-reviewed open access journal published by Ubiquity Press.