A user-friendly manual to apply the water availability predictive model from D'Agostini et al. 2023 and D'Agostini et al. 2024 to your phytolith dataset
Creators: name: D'Agostini, Francesca
name: Ruiz Perez, Javier
name: Lancelotti, Carla
Publisher: name: Zenodo
Year: 2025
Subjects: phytoliths
millets
water availability
sorghum
pearl millet
finger millet
irrigation
rainfed

description: List of contents:



What is this manual for

What you need to know before starting

You can apply the full model if…

Prepare your datasets         

Apply the model

Interpret the result

Make predictions on plant water availability using a subset of morphotypes

References

Do you need help? Contact us!


Required software and materials:



A spreadsheet editor capable of supporting .xlsx and .csv formats (e.g., LibreOffice https://www.libreoffice.org/).

R (https://cran.r-project.org/mirrors.html) or RStudio (https://dailies.rstudio.com/).

A dataset of modern phytoliths (training dataset) derived from plants grown under controlled watering conditions, where samples of both irrigated and non-irrigated plants are available. Access links are provided below or download the file “Modern.csv” directly from this upload.

A dataset of archaeological phytoliths (testing dataset).


1. What is this manual for 

Between 2021 and 2024, we developed a predictive model to reconstruct water availability for drought-tolerant crops  in archaeological contexts using phytoliths. The model is intended to build hypotheses on the presence or absence of irrigation systems. The model's output provides a measure of the likelihood that the archeological phytolith assemblage under observation had formed from irrigated plants. In this step-by-step guide, we show you how to apply the model to your dataset and the theoretical framework behind it. 

The predictive model has been previously published in two articles:



In D'Agostini et al. (2023, https://doi.org/10.1016/j.revpalbo.2022.104783), the production of phytoliths in three crop species was analyzed and modeled in relation to water availability, based on data from one experimental cultivation (conducted in 2019), and applied to an archaeological case study from the Indus Valley Civilization. 

D'Agostini et al. (2024, https://doi.org/10.1007/s00334-024-01012-9) combines the experimental results of 2019 with those obtained from an additional cultivation (replica) conducted in 2020. The 2019 and 2020 cultivations were conducted in the same location and under the same experimental conditions. Therefore,either dataset included in the 2023 and 2024 publications is valid as training data, but the 2024 dataset is more complete and based on a larger number of samples (361 entries).


2. What you need to know before starting

The model was developed using three C4 cereal species: finger millet (Eleusine coracana (L.) Gaertn.), pearl millet (Cenchrus americanus (L.) Morrone, syn. Pennisetum glaucum (L.) R.Br.), and sorghum (Sorghum bicolor (L.) Moench).

Ten landraces of each species from Ethiopia, Kenya, Pakistan, Sudan, and Tanzania were tested. For further information on the origin of the grains as well as on growth and productivity data, please refer to the supplementary database “FileS1-Experimental cultivation.xlsx” in D’Agostini et al. (2024) at the following link https://doi.org/10.5281/zenodo.10246367.

Each of the landraces analyzed was subjected to two distinct water treatments to simulate rainfed conditions in an arid context and irrigated conditions that also acted as a control. For data regarding the amount of water received by each replicate in total and weekly, harvest and flowering time, and additional experimental data refer to “FileS1-Experimental cultivation.xlsx” at the following link https://doi.org/10.5281/zenodo.10246367. To understand the basis on which we set the treatments, refer to the methodology explained in both D'Agostini et al. (2023) and D’Agostini et al. (2024).

Phytolith analysis was conducted on the leaves and chaff of 4 finger millet, 5 pearl millet, and 5 sorghum landraces. The morphotypes identified and counted include: ACUTE BULBOSUS, BILOBATE, BLOCKY, BULLIFORM FLABELLATE, CROSS, ELONGATE SINUATE clavate, ELONGATE SINUATE crenate, ELONGATE DENTATE, ELONGATE ENTIRE, ELONGATE SINUATE, POLYLOBATE, RONDEL, SADDLE, STOMA, TRACHEARY. For information on identification, classification, and photos of these morphotypes, please refer directly to the supplementary file “FileS2-Phytolith data.xlsx” (sheet Metadata - Phyto. classification) in D’Agostini et al. (2024), at the following link https://doi.org/10.5281/zenodo.10246367. You can consult the raw phytolith counts, the calculation of concentrations (expressed as grams of dry tissue), and the percentages of the total phytoliths counted in the same file in sheet Data. 

3. You can apply the full model if…

The water availability model can be applied to archeological assemblages of phytoliths using our phytolith dataset for model training (as done in D'Agostini et al., 2023 for sites of the Indus Valley Civilization) in two ways. Either utilizing the complete training dataset, which includes all morphotypes, or a subset of morphotypes. Check the conditions listed below to determine if your archaeological dataset meets the criteria for using the complete training dataset or just a subset of morphotypes. 

You can use the complete training dataset if…



Your archaeological phytolith assemblage includes ACUTE BULBOSUS, BILOBATE, BLOCKY, BULLIFORM FLABELLATE, CROSS, ELONGATE SINUATE clavate, ELONGATE SINUATE crenate, ELONGATE DENTATE, ELONGATE ENTIRE, ELONGATE SINUATE, POLYLOBATE, RONDEL, SADDLE, STOMA, TRACHEARY.

C₄ Panicoideae and Chloridoideae are considered predominant for the agro-environmental context/setting under study (e.g., in the Sahel regions of Africa where millets constitute the predominant part of agricultural production).

Finger millet, pearl millet, and sorghum were likely present in the past because they have been identified through phytolith analysis and/or alternative means (e.g., charred grains, isotopic signal of C₄). 


…otherwise you can apply the model using a subset of morphotypes if:



It is reasonable to assume that if most of the morphotypes listed here are part of your archaeological assemblage, then the presence of C₄ cultivated crops is likely. Therefore, by selecting only the phytoliths included in this list from your archaeological assemblage, the model can be applied. In addition, take into account that phytoliths from C₃ cereals are also effective in predicting water (Ermish and Boomgarden 2022; Jenkins et al., 2016, 2020; Madella et al., 2009). If C₃ cereals are identified in the archaeological context under study, morphotypes that are redundant in both C₃ and C₄ cultivated crops but that cannot be taxonomically distinguished can still be used (e.g., RONDEL).

You work in areas dominated by C₄ grasses but have a few of the listed morphotypes. In this case, the absence of certain morphotypes could be due to: (a) taphonomic factors that affect some morphotypes more than others (see for example Cabanes and Shahack-Gross 2015). In this instance, we suggest removing the missing morphotypes from the training dataset and apply the model using a subset of morphotypes; (b) the fact that not all three cultivated species used in the training model are/were present but the presence of at least one of these species can be confirmed through a combination of diagnostic morphotypes and/or other proxies.


If you are going to use the full model, proceed to sections 4-6. If you need to use a subset of morphotypes or only one morphotype (or several morphotypes but analyzed individually) for your prediction, then proceed directly to section 7. 

4. Prepare your datasets

Two datasets will be needed: the dataset of modern morphotypes to be used as training data and the dataset of morphotypes derived from an archaeological context to be used as test data.

The two datasets must meet the following criteria:



They should be two separate CSV files.

Both files must have the same column headers and number of columns. This means that the list of phytoliths identified in the archaeological record that are not present in the modern training dataset should be excluded (e.g., those clearly associated with C₃ like SPHEROID and CRENATE).

Both datasets must express the concentrations of morphotypes as a percentage of the total phytoliths counted (no absolute counts, no concentration). If some of the columns have been excluded from the analysis, it is essential to retain the percentages as they are (even if the total no longer adds up to 100%) so that the relative abundance of the phytoliths excluded from the analysis is still accounted for.


Prepare the training dataset

The training dataset can be downloaded directly from this guide (below files “Modern.csv”) or from the supporting data of D'Agostini et al. (2024) “FileS2-Phytolith data.xlsx" at https://doi.org/10.5281/zenodo.10246367. In the Data sheet, select and copy only the columns "Sample_Acronym", "Species", "Treatment", and "Origin" to obtain a unique identifier for each modern sample. Subsequently, select the columns corresponding to the percentages of each morphotype: TriPER, SaddlePER, RondelPER, BulliParallelPER, BulliFlabellatePER, BulliTOTPER, ClavatePER, SinuatePER, EntirePER, CrenatePER, StomaPER, CrossPER, BilobatePER, PolylobatePER. Create a new CSV with the column copied from the training dataset  and name it: Modern.csv. Be sure to use the filing name suggested, so that the R code will run without need to be modified. If you prefer to use a different file naming, make sure you change them in the R code as well, when appropriate.

Prepare the test dataset (archaeological)

Select the same morphotypes in your archaeological dataset and proceed with the creation of a CSV file, analogous to “Modern.csv”, with identical columns and name it: Archaeological.csv. Since the information regarding Species and Treatment is not available for the archaeological samples (and specifically Treatment is the variable we want to make predictions on), fill the Species column with a label of your choice that seems representative of the sample and useful for interpretation purposes. In D’Agostini et al. (2023), for example, we filled it with the type of analyzed context (ash, floor, fireplace, pit, street, and plaster), since we expected possible differences in the treatment of the plants depending on the context. However, the Treatment column can be filled with any label, which can be the same for all samples (e.g.,“archeo”), or different (e.g, context, site name, layer, etc). Do not use  “WW”, “WS”, “0” or “1”, which are the characters associated with the “irrigation” and “rainfed” treatments in the training dataset. To see, as an example, how the Archaeological.csv is organized in D’Agostini et al. (2023), please refer to “FileS1.2-Archaeological.csv” in the supplementary material available at https://doi.org/10.5281/zenodo.7120448.

5. Apply the model 

Our predictive model consists of two phases:



A binomial Generalised Linear Model with stepwise selection to identify the best explanatory variables (morphotypes) of the training dataset that can distinguish between watering treatments.

A logistic regression to predict the watering treatments, either WW or WS, using morphotype percentages selected by the stepwise selection.


Procedure to follow: 

1) Open R or if you prefer RStudio. For our analyses and for this tutorial we relied on the use of RStudio and all the following instructions refer to applying the model in RStudio, but you can run the exact same code directly in R as well.

2) Download the R code “FileS3-Rscript.R” from this guide or from the supplementary material of D'Agostini et al. (2023) at the following link  https://doi.org/10.5281/zenodo.7120448, and open it in RStudio or follow the R code explained in the following sections. The code has been written and tested using R version 4.2.3 and RStudio 2022.07 for MacOS 10.15.

3) Open the two files Modern.csv and Archaeological.csv in RStudio. To open the code you can either run the following code (copied from  “FileS3-Rscript.R”) or manually load the two files in Environment → Import Dataset → From text (base).

Code to open/load the datasets in RStudio:

setwd("/Users/Francesca/desktop") 

### In our case, the files are saved on the computer's desktop (/Users/Francesca/desktop), but modify this section by indicating where your datasets are saved in your computer (e.g., /Users/jrp/Documents/Model)

data<-read.table("Modern.csv",sep=",",header=TRUE)

dataset<-read.csv(file="Archaeological.csv", sep=",")

### load the following libraries by running:

library(ggplot2)

library(ggpubr)

library(ggrepel)

library(tidyverse)

library(svglite)

library(MASS)

4) Run the code from ###LOGISTIC REGRESSION### onwards copying and pasting the following section into the console and pressing Enter, or alternatively selecting the code in the script and clicking Ctrl/Cmd + Enter. 

#LOGISTIC REGRESSION###

#MODERN

data<-read.table("Modern.csv",sep=",",header=TRUE)

data[, c(3)] <- sapply(data[, c(3)], as.numeric) 

### all variables (columns) are transformed into numeric 

data$Treatment [data$Treatment == "1"] <- 0

data$Treatment [data$Treatment == "2"] <- 1

data$Treatment<- as.integer(data$Treatment) 

### the values of “Treatment” are modified so that all “WW” correspond to the value 1 and “WS” correspond to 0

full.model <- glm(Treatment ~ TriPER + SaddlePER + RondelPER + BulliParalPER + BulliFlabPER+ StomaPER + CrossPER + BiloPER + PolyPER,data = data, family = binomial(link="logit"))

summary(full.model)

### the logistic regression model is tested without excluding any predictors = logistic regression (family = binomial(link="logit")) is run, using all morphotypes of the dataset as predictors of “Treatment”

#Stepwise selection

step.model<- stepAIC(full.model, direction = "both", trace = FALSE)

summary(step.model)

### stepwise selection is applied to determine if the predictive variability can be improved by using a subset of predictors = a stepwise selection function (stepAIC) is applied to the full.model to select the best predictors (i.e., morphotypes), making the model more robust (i.e., showing lower AIC and better performance metrics)

#ARCHAEOLOGICAL

dataset<-read.csv(file="Archaeological.csv", sep=",") 

### archaeological dataset loading

dataset[, c(1,3,4)] <- sapply(dataset[, c(1,3,4)], as.numeric) 

### all variables are transformed into numeric factors as done with the modern training dataset

X<-dataset[,c(5:14)] 

### an additional dataset (X) is created, which contains only the variables corresponding to the morphotypes count values in percentage (i.e., excluding the descriptive columns Sample_Acronym, Species, Treatment, and Origin)

p <- predict(step.model, X,  type="response") 

### the predictive model is applied using the predict function, which allows the user to use the results of the step.model to make predictions about the samples in the X dataset

output <- cbind(dataset, p)

### create a dataframe (output) that combines the archaeological test dataset with these probability predictions

#Plotting results

Regre<-ggplot(output, aes(x = output$Origin_Spe, y = output$p, label=output$Sample)) + 

  geom_point(aes(color=factor(output$Origin)), size=3) +

### placing the sample origins on the x-axis while the prediction result p is on the y-axis

  facet_wrap(~Species, scales="free_x") + 

 ### dividing samples by archaeological site of origin

  geom_text_repel(size=2) + 

### each sample, representing an archaeological phytolith assemblage, is depicted by a point

labs(title="",x="",y="probability of WW") 

jpeg("Prob.jpg", width=2500, height=2000, res=300)

### save high-definition jpeg figures of the plotted results (it is also possible to save the graphs in TIFF or SVG by replacing “jpeg” in the code with “tiff” or “svg”)

Regre +theme_bw()+ scale_color_brewer(palette="Set2")+ theme(legend.title=element_blank(), legend.position="bottom", aspect.ratio=1,text=element_text(size=15),axis.title=element_text(size=13),axis.text.x=element_blank(),axis.ticks.x=element_blank())

dev.off()

### aesthetics of the graph (e.g., the chosen color palette, a white background with grids, and the font sizes)

6. Interpret the result

The graph represents the output obtained after running the code using phytolith data from multiple archeological sites from the Indus Valley Civilization (Kanmer, Shikarpur, Alamgipur, and Harappa) (see D’Agostini et al., 2023).

Each point in the graph represents an archaeological sample. The y-axis corresponds to the probability that the assemblage tested is composed of phytoliths from plants grown under well-watered conditions, i.e., likely irrigated in areas with scarce rainfall.  In our case, WS corresponds to an environmental setting where, during the growing season, an average of 150-155 mm of water was received. The y-axis ranges from 0, indicating no probability of WW, to 1, representing 100% probability that an assemblage derives from plants grown under well-watered conditions (over 600 mm of water received). The inflection point is at 0.5, where there is an equal probability (50%) that the plants were grown either under well-watered conditions or water scarcity. If the sample distributes between 0-0.5 on the y-axis, it possibly originates from an assemblage of plants grown under scarce water conditions, stimulating a drought response, as is the case with the majority of our samples depicted in the graph. Whereas, if the sample falls within the graph area between 0.5 and 1, it is more likely that the assemblage derives from plants grown under good water conditions.

7. Making predictions on plant water availability using a subset of morphotypes

If both the conditions of your context and your archaeological dataset indicate the need to use a reduced training dataset, there are two possible options: 



Reduce the number of morphotypes using a subset assemblage of morphotypes in the training dataset, while replicating the same methodology explained above (sections 4, 5, and 6). The output will be a percentage corresponding to the probability that the archaeological assemblage of phytoliths originates from plants grown in a well-irrigated context. If this is the case go to section 7.1.

Base the prediction on one single key morphotype. The logistic regression formula can be applied to a single morphotype predictor to obtain the specific probability percentage number (quantitative result), considering that, although the prediction result is based on the concentration of a single morphotype, the training dataset still accounts for the entire assemblage.If this is the case go to sections 7.2 and then 7.2.2. 


7.1) Reduce the number of morphotypes both in training and testing datasets

From the modern dataset, remove all morphotypes that are not present in  the archaeological assemblage under study (see section 3 for discussion) from Modern.csv dataset. Subsequently, follow the model application procedures explained in sections 5 and 6 but, in RStudio, modify the code line “Full.model <- glm(Treatment ~ TriPER + SaddlePER + RondelPER + BulliParalPER + BulliFlabPER + StomaPER + CrossPER + BiloPER + PolyPER, data = data, family = binomial(link="logit"))” by deleting from the independent variables (in this case, the percentages of morphotypes highlighted in yellow) those that have been removed from the training dataset. The resulting output will be the same as shown in sections 5 and 6, but the model will be based on a more restricted number of predictors.

Keep in mind that even if specific columns of morphotypes are removed from the training data, the relative abundances represented there do not change, as the total number of phytoliths counted for the training remains the same. Therefore, the analysis will take into account that the three species also produce morphotypes that, due to pre- or post-depositional taphonomic reasons, have not appeared in the archaeological record.

7.2) Using individual morphotypes (i.e., just one or one at a time) to make predictions

From now on, we will use RONDEL as an example, but any of those morphotypes included in the full model (i.e., ACUTE BULBOSUS, BILOBATE, BLOCKY, BULLIFORM FLABELLATE, CROSS, ELONGATE SINUATE clavate, ELONGATE SINUATE crenate, ELONGATE DENTATE, ELONGATE ENTIRE, ELONGATE SINUATE, POLYLOBATE, RONDEL, SADDLE, STOMA, TRACHEARY) can be used. RONDEL has been chosen as an example  because it has shown to be crucial in predicting WW and WS conditions (see D'Agostini et al., 2024).

Option 7.2.1) Estimate the probability by graph comparison 

To compare the concentration of RONDEL (or any other morphotype of interest) in your dataset, follow these steps:

1) Open R or if you prefer RStudio. The following instructions refer to applying the model in RStudio, but you can run the exact same code directly in R as well.

2) Download the R code “FileS3.R” file from this upload or from the supplementary information of D'Agostini et al. (2024) at  https://doi.org/10.5281/zenodo.10246367, and open it in RStudio. The code has been written and tested using R version 4.2.3 and RStudio 2022.07 for MacOS 10.15.

3) Download “Phitoliths.csv” form this upload or the “FileS2-Phytolith data.xlsx” data from the supplementary information of D'Agostini et al. (2024) at  https://doi.org/10.5281/zenodo.10246367.

4) In the Data sheet of “FileS2-Phytolith data.xlsx”, select only the columns "Sample_Acronym", "Species", "Treatment", and "Origin" to obtain a unique identifier for each modern sample. Subsequently, select the columns corresponding to the percentages of each morphotype (all of them, the morphotype of interest will be selected at a later stage): TriPER, SaddlePER, RondelPER, BulliParallelPER, BulliFlabellatePER, BulliTOTPER, ClavatePER, SinuatePER, EntirePER, CrenatePER, StomaPER, CrossPER, BilobatePER, PolylobatePER. Create a new CSV file and name it: Phytoliths.csv. If you are interested in only one of the three species under examination, copy only the data related to the species of interest. For more details on the preparation of the dataset, see section 4.1.

5) Open Phytoliths.csv in RStudio. To open the code you can either run the following code (copied from “FileS3.R”) or manually load the file in Environment → Import Dataset → From text (base).

Code to open/load the datasets in RStudio:

setwd("/Users/Francesca/desktop")

### In our case, the files are saved on the computer's desktop (/Users/Francesca/desktop), but modify this section by indicating where your datasets are saved in your computer (e.g., /Users/jrp/Documents/Model)

data<-read.table("Phytoliths.csv",sep=",",header=TRUE)

### load the following libraries by running:

library(ggplot2)

library(ggpubr)

library(ggrepel)

library(tidyverse)

library(svglite)

library(MASS)

6) To obtain the concentration values of RONDEL only (or the morphotype of your choice) corresponding to WS and WW probabilities, modify the code #4 Plotting results according to the morphotypes you want to abalyse. For more details on the function of the code lines to modify, see the file PDF of this upload "User friendly manual.pdf".

Run the modified code from ###BEST FITTED MODEL & LOGISTIC REGRESSION### onwards copying and pasting the following section into the console and pressing Enter, or alternatively selecting the code in the script and clicking Ctrl/Cmd + Enter. 

###BEST FITTED MODEL & LOGISTIC REGRESSION###

data<-read.table("Phytoliths.csv",sep=",",header=TRUE)

data[, c(11,59:75)] <- sapply(data[, c(11,59:75)], as.numeric)

### all variables (columns) selected are transformed into numeric

X<-data[,c(1,8,10,11,24,59:75)] 

X$Treatment [X$Treatment == "1"] <- 0

X$Treatment [X$Treatment == "2"] <- 1

data$Treatment <- as.numeric(data$Treatment)

### the values of “Treatment” are modified so that all “WW” correspond to the value 1 and “WS” correspond to 0

#1 Full model

full.model <- glm(Treatment ~ RatioS/F_percentage + TriPER + SaddlePER + RondelPER  + BulliParallelPER + BulliFlabellatePER + BulliTOTPER + ClavatePER + SinuatePER + EntirePER + DentatePER + CrenatePER + ElongTOTPER + StomaPER + BilobatePER + PolylobatePER, data = X, family = binomial(link="logit"))

summary(full.model)

### the logistic regression model is tested without excluding any predictors = logistic regression (family = binomial(link="logit")) is run, using all morphotypes of the dataset as predictors of “Treatment”

#2 Stepforward model

step.model<- stepAIC(full.model, direction = "both", trace = FALSE)

summary(step.model)

### stepwise selection is applied to determine if the predictive variability can be improved by using a subset of predictors = a stepwise selection function (stepAIC) is applied to the full.model to select the best predictors (i.e., morphotypes), making the model more robust (i.e., showing lower AIC and better performance metrics)

#3 Predicting

p <- predict(step.model, X,  type="response")

### the predictive model is applied using the predict function, which allows the user to use the results of the step.model to make predictions about the samples in the X dataset

#4 Plotting results 

Regre<-ggplot(X, aes(x = X$RondelPER , y = X$Treatment)) + 

### placing the RONDEL concentration percentages on the x-axis while the Treatment is represented on the y-axis

 geom_point(size=3)+

 geom_smooth(color="black", method = "glm", method.args = list(family = "binomial"), se = FALSE)+

###  the trend of the probability of a given concentration of RONDEL belonging to WW or WS is represented by the curve coded with geom_smooth

  labs(title="",

       x="Rondels",

       y="probability of WW") 

jpeg("RondelPER.jpg", width=2500, height=2000, res=300)

### save high-definition jpeg figures of the plotted results (it is also possible to save the graphs in TIFF or SVG by replacing “jpeg” in the code with “tiff” or “svg”)

Regre + theme_bw()+ theme(legend.title=element_blank(legend.position="bottom", aspect.ratio=1,text=element_text(size=15), axis.title=element_text(size=13),)

dev.off()

aesthetics of the graph (e.g., the chosen color palette, a white background with grids, and the font sizes)

The graph represents the output obtained after running the code. On the x-axis is the relative abundance as a percentage of RONDEL for each sample, while on the y-axis is the corresponding water treatment. The curve represents the probability that a specific concentration of RONDEL comes from plants grown under well-watered conditions, i.e., likely irrigated in areas with scarce rainfall. In our case, WS corresponds to an environmental setting where, during the growing season, an average of 150-155 mm of water was received. The y-axis ranges from 0, indicating no probability of WW, to 1, representing 100% probability that an assemblage derives from plants grown under well-watered conditions (over 600 mm of water received). The inflection point is at 0.5, where there is an equal probability (50%) that the plants were grown either under well-watered conditions or water scarcity. The curve represents the trend of probability, and helps us visualize in which range of RONDEL concentrations there is a probability that they derive from WW or WS plants, depending on its positive or negative trend, and which value corresponds to the inflection point of 0.5 on the y-axis. In our example shown below, if RONDEL concentration is between 0-1%, it possibly originates from plants grown irrigated. Whereas, if  RONDEL concentration exceeds 1% of the total assemblage, it is more likely that the assemblage derives from plants grown under WS conditions. This means also that you can compare the relative concentrations of RONDEL in your archaeological assemblage with the outputs of this logistic regression derived from our modern dataset and hypothesize whether they originate from a rainfed context or not: if your RONDEL exceed 1%, they most likely come from a rainfed context.

In D'Agostini et al. (2024), the results included graphs representing the outcome of the logistic regressions (the trend of the prediction curve) for all the morphotypes examined and considering all the three species simultaneously and separately (Figure 6 in the paper and Figures S5 to S15 in the supplementary information available at https://doi.org/10.5281/zenodo.10246367).

Option 7.2.2) Applying the logistic regression formula to calculate the probability

1) Open R or if you prefer RStudio. For our analyses and for this tutorial we relied on the use of RStudio and all the following instructions refer to applying the model in RStudio, but you can run the exact same code directly in R as well.

2) Download the R code “FileS3.R” file from this upload or from the supplementary information of D'Agostini et al. (2024) at https://doi.org/10.5281/zenodo.10246367, and open it in RStudio. The code has been written and tested using R version 4.2.3 and RStudio 2022.07 for MacOS 10.15.

3) Download “Phytoliths.csv” from this upload or the “FileS2-Phytolith data.xlsx” data from the supplementary information of D'Agostini et al. (2024) at https://doi.org/10.5281/zenodo.10246367.

4) In the Data sheet of “FileS2-Phytolith data.xlsx”, select only the columns "Sample_Acronym", "Species", "Treatment", and "Origin" to obtain a unique identifier for each modern sample. Subsequently, select the columns corresponding to the percentages of each morphotype: TriPER, SaddlePER, RondelPER, BulliParallelPER, BulliFlabellatePER, BulliTOTPER, ClavatePER, SinuatePER, EntirePER, CrenatePER, StomaPER, CrossPER, BilobatePER, PolylobatePER. Create a new CSV file and name it: Phytoliths.csv. If you are interested in only one of the three species under examination, copy only the data related to the species of interest. For more details on the preparation of the dataset, see section 4.1.

5) Open Phytoliths.csv in RStudio. To open the code you can either run the following code (copied from “FileS3.R”) or manually load the file in Environment → Import Dataset → From text (base).

Code to run to open/load the datasets in RStudio:

setwd("/Users/Francesca/desktop")

### In our case, the files are saved on the computer's desktop (/Users/Francesca/desktop), but modify this section by indicating where your datasets are saved in your computer (e.g., /Users/jrp/Documents/Model)

data<-read.table("Phytoliths.csv",sep=",",header=TRUE)

### load the following libraries by running:

library(ggplot2)

library(ggpubr)

library(ggrepel)

library(tidyverse)

library(svglite)

library(MASS)

6) Run the code from ###BEST FITTED MODEL & LOGISTIC REGRESSION### onwards copying and pasting the following section into the console and pressing Enter, or alternatively selecting the code in the script and clicking Ctrl/Cmd + Enter.  

###BEST FITTED MODEL & LOGISTIC REGRESSION###

data<-read.table("Phytoliths.csv",sep=",",header=TRUE)

data[, c(11,59:75)] <- sapply(data[, c(11,59:75)], as.numeric)

### all variables (columns) selected are transformed into numeric

X<-data[,c(1,8,10,11,24,59:75)] 

X$Treatment [X$Treatment == "1"] <- 0

X$Treatment [X$Treatment == "2"] <- 1

data$Treatment <- as.numeric(data$Treatment)

### the values of “Treatment” are modified so that all “WW” correspond to the value 1 and “WS” correspond to 0

#1 Full model

full.model <- glm(Treatment ~ RatioS/F_percentage + TriPER + SaddlePER + RondelPER  + BulliParallelPER + BulliFlabellatePER + BulliTOTPER + ClavatePER + SinuatePER + EntirePER + DentatePER + CrenatePER + ElongTOTPER + StomaPER + BilobatePER + PolylobatePER, data = X, family = binomial(link="logit"))

summary(full.model)

### the logistic regression model is tested without excluding any predictors = logistic regression (family = binomial(link="logit")) is run, using all morphotypes of the dataset as predictors of “Treatment”

#2 Stepforward model

step.model<- stepAIC(full.model, direction = "both", trace = FALSE)

summary(step.model)

### stepwise selection is applied to determine if the predictive variability can be improved by using a subset of predictors = a stepwise selection function (stepAIC) is applied to the full.model to select the best predictors (i.e., morphotypes), making the model more robust (i.e., showing lower AIC and better performance metrics)

The result is a table (see on files below "User friendly manual.pdf"). The values needed in the formula to calculate the probability are β0 and β1. Replace the coefficients obtained from running the model in the logistic regression formula below to obtain the probability that your samples were grown under WW conditions:

p(x) = 11 + e^-(β0 +β1x)

p(x) represents the probability of a specific morphotype concentration that the plants producing it grew under WW conditions.

β0 corresponds to the intercept estimate value generated by the model (see the figure above of the console output). 

β1 corresponds to the estimate of the morphotype of interest (see the figure above of the console output). In this example where we used RONDEL, it corresponds to the RondelPER estimate.

The resulting value represents the probability that a specific morphotype concentration has been produced by a plant that grew under WW conditions. It reflects the percentage concentration obtained from the archaeological dataset, which serves as the input, specifically, the RONDEL percentage in the testing dataset. This value is derived from the phytolith analysis following slide-based counting.

References

Cabanes D, Shahack-Gross R. 2015. Understanding Fossil Phytolith Preservation: The Role of Partial Dissolution in Paleoecology and Archaeology. PLOS ONE 10(5): e0125532.

D’Agostini F, Ruiz-Pérez J, Madella M, Vadez V, Kholova J, Lancelotti C. 2023. Phytoliths as indicators of plant water availability: The case of millets cultivation in the Indus Valley civilization. Review of Palaeobotany and Palynology 309.

D’Agostini F, Ruiz Pérez J, Madella M, Vadez V, Lancelotti C. 2024. Predicting plant water availability from phytolith assemblages: an experimental approach for archaeological reconstructions in drylands. Vegetation History and Archaeobotany.

Ermish BJ, Boomgarden SA. 2022. Identifying water availability with maize phytoliths in Range Creek Canyon, Utah. Journal of Archaeological Science: Reports 41.

Jenkins E, Jamjoum K, Nuimat S, Stafford R, Nortcliff S, Mithen S. 2016. Identifying ancient water availability through phytolith analysis: An experimental approach. Journal of Archaeological Science 73: 82–93.

Jenkins EL, Predanich L, Al Nuimat SAMY, Jamjoum KI, Stafford R. 2020. Assessing past water availability using phytoliths from the C4 plant Sorghum bicolor: An experimental approach. Journal of Archaeological Science: Reports 33.

Madella M, Jones MK, Echlin P, Powers-Jones A, Moore M. 2009. Plant water availability and analytical microscopy of phytoliths: Implications for ancient irrigation in arid zones. Quaternary International 193: 32–40.

 

Do you need help? Contact us!

Francesca D’Agostini 

francesca.dagostini@upf.edu 

 

Javier Ruiz-Pérez

javier.ruizperez.academic@gmail.com