order_of_code-doc.Rmd

---
title: "Adult Mortality in the Metropolis of London 1100–1850: a Bayesian View Based on Osteological Data"
subtitle: 'Supplement: Code structure, data source and processing'
author:
- "Nils Müller-Scheeßel^[Institute for Prehistoric and Protohistoric Archaeology - Kiel University nils.mueller-scheessel@ufg.uni-kiel.de]"
- "Katharina Fuchs^[Institute for Prehistoric and Protohistoric Archaeology - Kiel University, k.fuchs@ufg.uni-kiel.de]"
- "Christoph Rinne^[Institute for Prehistoric and Protohistoric Archaeology - Kiel University, crinne@ufg.uni-kiel.de]"
date: "`r format(Sys.time(), '%d. %B %Y')`"
output:
  pdf_document:
    toc: yes
    toc_depth: 3
  html_document:
    code_folding: hide
    toc: yes
    toc_depth: 3
    toc_float: yes
    number_sections: no
license: "CC-BY 4.0"
header-includes: null
papersize: a4
# knit rmdfig_caption: yes
urlcolor: blue
link-citations: yes
linkcolor: blue
bibliography: References.bib
csl: apa-6th-edition.csl
---

```{r knitr global options, include=FALSE}
# Global options for knitr
knitr::opts_chunk$set(echo=TRUE, fig.align="center", duplicate.label = "allow")
```

# Prerequisites {-}

The calculations were made in R using R-Studio. The structure of the code is essentially based on the structure of the text. The raw code is in the file ```order_of_code.R```. The file extended with Markdown is ```order_of_code-doc.RMD``` and the file ```order_of_code-doc.html``` is generated from it.

Note: The base path for rmd files is the folder in which they are located, not the r-project. Consequently, ```order_of_code.R``` and ```order_of_code-doc.RMD``` are both located in the root folder of the project.

> Depending on the hardware, the subsequent code can run for several hours or even a few days.

> Install "Just Another Gibbs Sampler" (JAGS) [@ref_103719] if you want to run the Bayesian analyses anew.
> Version 4.3 - as used here - can be downloaded in pre-compiled form for a number of OS here: https://sourceforge.net/projects/mcmc-jags/ . 
> The manual can be found here: https://people.stat.sc.edu/hansont/stat740/jags_user_manual.pdf

The code makes extensive use of the function ```source``` to call external code. Thus, the main part of the code remains slim, well structured and readable. 

Install required packages, set some options and link the sources for the helper functions. 

Remark: Depending on your R version the package osmplotr may be installed from github using ```devtools::install_github ("ropensci/osmplotr")```. 

```{r Packages and sources, echo=TRUE}
require(pacman) || install.packages("pacman")
pacman::p_load(coda, cowplot, demogR, dplyr, flexsurv, ggplot2, ggrepel, grid, 
               gridExtra, HMDHFDplus, kableExtra, Metrics, mortAAR, osmdata, 
               psych, readxl, reshape2, rjags, runjags, sf, tidyr, ggspatial)

options(scipen = 999)
options(dplyr.summarise.inform = FALSE)

source("./functions/gomp_MLE.R")
source("./functions/gomp_MLE_adapted.R")
source("./functions/gomp_MLE_interval.R")
source("./functions/gomp_anthr_age.R")
source("./functions/gomp_anthr_age_r.R")
source("./functions/gomp_known_age.R")
source("./functions/gomp_known_age_r.R")
source("./functions/helper_functions.R")
source("./functions/lt_MC.R")
source("./functions/lt_MC_Gomp.R")
RNGkind("L'Ecuyer-CMRG") # conservative random number generator to avoid periodicity
```
Important for saving time:  Decide to run extensive code anew (app. 6 h +). In addition, you can set the folder for preprocessed files.

```{r run code anew, echo=TRUE}
runCodeNew <- FALSE
#runCodeNew <- TRUE

# Ask for credentials of the Human Mortality Database if the code runs anew
if (runCodeNew){
  HMD_username <- readline(prompt = "Enter username: ")
  HMD_password <- readline(prompt="Enter password: ")
  credentials <- c(HMD_username, HMD_password)
}

# Specify filename prefix for saved files and create a folder if needed:
saveFileDir = "preprocessed_files"
if (saveFileDir %in% list.files(getwd())) {
  # Dir exists
}else{
  dir.create(file.path(".", saveFileDir), showWarnings = FALSE )
}
```

\newpage
# Chapter 01 Introduction

### Figure 1: Exemplary life table curves generated by Gompertz functions with different values for the $\beta$ parameter.

```{r Fig. 1: Gompertz distribution, echo=TRUE, warning=FALSE, include = TRUE, code = readLines("./chapter_01_introduction/gompertz_distribution.R")}
#source("./chapter_01_introduction/gompertz_distribution.R")
```

\newpage
# Chapter 02 Materials and methods

### Figure 3: Hazard (m) of the population of England and Wales, 1841. The blue rectangle in the upper panel corresponds to the extent of the lower panel. The turning point of 12 years is marked (from Human Mortality Database).

```{r Fig. 3: Hazard curve, echo=TRUE, warning=FALSE, include = TRUE, code = readLines("./chapter_02_materials_and_methods/hazard_curve.R")}
#source("./chapter_02_materials_and_methods/hazard_curve.R")
```

\newpage
# Chapter 03 Data

### Figure 4: Major cemeteries in Greater London 1100–1850 used in the present study. 

```{r Fig. 3: Map of London, echo=TRUE, warning=FALSE, include = TRUE, code = readLines("./chapter_03_data/London_places.R")}
#source("./chapter_03_data/London_places.R") 
```

### Table 2: Overview of major cemeteries in Greater London 1100–1850 used in the present study.

```{r Table 2: Overview of included cemeteries, echo=TRUE, warning=FALSE, results="asis"}
read.table("chapter_03_data/london_cemeteries.txt", header = T, sep = "\t") %>%
  knitr::kable(., col.names = c("Map no.",	"Name",	"Period",	"Excavation year",	
                                "Social character",	"Absolute n",	"n analysed",	
                                "n >= 12 years",	"References")) %>%
  kableExtra::kable_styling(latex_options = "HOLD_position") %>% unclass() %>% cat()
```

\newpage
### Figure 5: Population development of London, compiled from @ref_268519, 39 table 1; @ref_90775, 41; 179 table 5.7; @ref_115940, 655--657.

```{r Fig. 5: Population development, echo=TRUE, warning=FALSE}
source("./chapter_03_data/London_population.R")
grid::grid.newpage()
grid::grid.draw(rbind(london_pop1, london_pop2))
```

### Footnote 6: Re-calculation of population increase rates of London from @ref_278969.
Calculated in ./chapter_03_data/London_population.R

```{r Footnote 6: Razzell/Spence 2007, warning=FALSE}
knitr::kable(razz_df)%>%
  kableExtra::kable_styling(latex_options = "HOLD_position")
```

\newpage
# Chapter 04 Results

## Simulations

### Figure 6: Comparison of algorithms to estimate the original Gompertz $\beta$ from simulated data with known age-at-death (n = 1,000).

```{r Fig. 6. Simulations retrieving Gompertz parameters, fig.height=8.5, fig.width=6, warning=FALSE}
source("./chapter_supplement/simulations_run.R")
ggsave(filename = "fig06_lt_sim_plots.pdf",width = 8, height = 11.5,plot = gridExtra::grid.arrange (grobs = plot_list_shapes, ncol = 2), device = cairo_pdf,path = "documented"
)
```

\newpage

### Figure 7: Difference between observed and estimated Gompertz $\beta$-values by different algorithms with known age-at-death (n = 1,000).

```{r Fig 7. Simulations with difference between observed and estimated Gompertz parameters, fig.height=8.5, fig.width=6, warning=FALSE}
source("./chapter_supplement/simulations_run.R")
ggsave(filename = "fig07_lt_sim_plots.pdf", width = 8, height = 11.5, plot = gridExtra::grid.arrange (grobs = plot_list_diff, ncol = 2), device = cairo_pdf, path = "documented"
)
```

\newpage
### Table 3: Root mean square errors (RMSE) for different algorithms for fitting known age-at-death, in ascending order.

```{r T3 Simulation evaluation via RMSE}
# table of RMSEs
kable(rmse_result[order(rmse_result$RMSE) ,]) %>%
  kableExtra::kable_styling(latex_options = "HOLD_position") 
write.table(rmse_result[order(rmse_result$RMSE) ,], file = "./documented/table03_rmse_known_age.txt", 
            sep="\t", quote = FALSE)
```

\newpage
### Figure 8: Comparison of algorithms to estimate the original Gompertz $\beta$ from simulated data with osteological age categories (n = 1,000).

```{r Fig. 8. Simulations retrieving Gompertz parameters with age estimations, fig.height=8.5, fig.width=6, warning=FALSE}
source("./chapter_supplement/simulations_run.R")
ggsave(filename = "fig08_lt_sim_plots.pdf", width = 8, height = 8, plot = gridExtra::grid.arrange (grobs = plot_list_estim_shapes, ncol = 2), device = cairo_pdf,path = "documented"
)
```

\newpage
### Table 4: Root mean square errors in different algorithms to estimate the original Gompertz  $\beta$ from simulated data with osteological age categories, in ascending order.

```{r T4 Simulation evaluation via RMSE: anthropological categories}
# table of RMSEs
kable(rmse_estim_result[order(rmse_estim_result$RMSE) ,]) %>%
  kableExtra::kable_styling(latex_options = "HOLD_position") 
write.table(rmse_estim_result[order(rmse_estim_result$RMSE) ,], file = "./documented/table04_rmse_osteo_age.txt", 
            sep="\t", quote = FALSE)
```

### Figure 9: Simulation of population increase with known age-at-death and Maximum Likelihood Estimation (MLE) (top four) and osteological estimates, Bayesian model and including rate of increase (bottom four).

```{r Fig. 9. Simulation of population increase, fig.height=6, fig.width=8, warning=FALSE, include = TRUE, code = readLines("./chapter_04_results/simulations_pop_incr_run.R")}
#source("./chapter_04_results/simulations_pop_incr_run.R")
```


## Written sources

Preprocessing of data used in figure 9: Estimated modal ages.

### Basic statistics

The data is referenced and aggregated in "./chapter_04_results/historical_lifetables.R". In this file, all records from individual preprocessing files located in "./liftables_preprocessed/" are ```sourced```.  The corresponding data files are stored in "./data/".

London_1728_1840.R, Mortality_bills_1728_1840.txt, Source: @ref_301236, 304 Table 6.5; > 100 years and < 1 year collapsed

```{r Basic statistics Mortality bills 1728-1840 age modes, warning=FALSE}
source("./chapter_04_results/historical_lifetables.R")
kable(london_1728_1840_ranges, 
      caption = "London Mortality bills 1728-1840.") %>%
  kableExtra::kable_styling(latex_options = "HOLD_position")
```

```{r Basic statistics Mortality bills 1728-1840 rate, warning=FALSE}
kable(london_1728_1840_ranges_r, 
      caption = "London Mortality bills 1728-1840, corrected for population growth.") %>%
  kableExtra::kable_styling(latex_options = "HOLD_position")
```

London_1841_raw_all.R, London_1841_raw.txt, Source: @ref_94043, 19 table q.

```{r Basic statistics census data for London from 1841, warning=FALSE}
kable(London_1841_ranges, 
      caption = "Census data for London from 1841.") %>%
  kableExtra::kable_styling(latex_options = "HOLD_position")
```

English_Mortality.R, wrigley_et_al_1997_england_1640-1809.txt, Source: @ref_73169, 290 table 6.19

```{r Basic statistics English mortality, warning=FALSE}
kable(eng_mort_ranges, 
      caption = "English mortality data.") %>%
  kableExtra::kable_styling(latex_options = "HOLD_position")
```

HMD_UK_ranges.R

The data from the Human Mortality Database (https://mortality.org/) were retrieved with a personal account using the R package HMDHFDplus. Therefore, we only provide the processed data here.  

```{r Basic statistics HMD UK, warning=FALSE}
kable(HMD_UK_ranges, caption = "Human Mortality Database UK.") %>%
  kableExtra::kable_styling(latex_options = "HOLD_position")
```

English_Peers.R, russell.txt, Source: @ref_45151, table 2

```{r Basic statistics English Peers, warning=FALSE}
kable(peers_ranges, caption = "English Peers") %>%
  kableExtra::kable_styling(latex_options = "HOLD_position")
```

Medieval_England.R, Christ_church_monks.txt, Source: @ref_54763, 28 table 2

```{r Basic statistics Christ Church monks, warning=FALSE}
kable(monks_ranges, caption = "Christ Church monks") %>%
  kableExtra::kable_styling(latex_options = "HOLD_position")
```

### Extended statistics

```{r Extended statistics written sources, warning=FALSE}
kable(london_1728_1840_result, 
      caption = "London Mortality bills 1728-1840.") %>%
  kableExtra::kable_styling(latex_options = c("HOLD_position","scale_down"))
kable(london_1728_1840_result_r, 
      caption = "London Mortality bills 1728-1840, corrected for population growth.") %>%
  kableExtra::kable_styling(latex_options = c("HOLD_position","scale_down"))
kable(London_1841_result, 
      caption = "Census data for London from 1841.") %>%
  kableExtra::kable_styling(latex_options = c("HOLD_position","scale_down"))
kable(eng_mort_result, caption = "English mortality data.") %>%
  kableExtra::kable_styling(latex_options = c("HOLD_position","scale_down"))
kable(HMD_UK_result, caption = "Human Mortality Database UK.") %>%
  kableExtra::kable_styling(latex_options = c("HOLD_position","scale_down"))
kable(peers_result, caption = "English Peers.") %>%
  kableExtra::kable_styling(latex_options = c("HOLD_position","scale_down"))
kable(monks_result, caption = "Christ Church monks.") %>%
  kableExtra::kable_styling(latex_options = c("HOLD_position","scale_down"))
```

## London cemeteries

The data is mainly hard coded in the file ./chapter_04_results/Wellcome_DB.R. 

Only St. Bride's crypt is excluded but available from the Museum of London upon request. For general information: [https://www.museumoflondon.org.uk](https://www.museumoflondon.org.uk) go for: Collections > Archaeology at the Museum of London > Wellcome Osteological Research Database > St. Bride's Church Fleet Street. 
If runCodeNew == TRUE  the file ./lifetables_processing/stbrides_crypt.R will ask for the location of the retrieved dataset (Excel sheet) and process the data. In any other case pre-processed data will be loaded. 

```{r cemeteries data pre-processed, warning=FALSE}
source("./lifetables_processing/stbrides_crypt.R")
source("./chapter_04_results/Wellcome_DB.R")
```

```{r wellcome results, warning=FALSE}
kable(wellcome_result) %>%
  kableExtra::kable_styling(latex_options = c("HOLD_position","scale_down"))
kable(wellcome_result_r, caption = "London cemeteries data, corrected for population growth.") %>%
  kableExtra::kable_styling(latex_options = c("HOLD_position","scale_down"))
```

### Figure 10: Estimated modal ages from written sources and osteological data compared, upper panel: without population growth correction, lower panel: with population growth correction. Horizontal bars indicate the time span the data point covers. Vertical bars indicate 95% HDI for credible ranges and are only displayed for small n, i.e. English Peers, Christ Church monks and osteological data.

```{r figure 10 modal ages plot, fig.height=8.5, fig.width=11, echo=TRUE, warning=FALSE, , include = TRUE, code = readLines("./chapter_04_results/english_wellcome.R")}
#source("./chapter_04_results/english_wellcome.R")
```

The following data overview is build during pre-processing in ./chapter_04_results/Wellcome_DB.R and saved to a textfile (sep = \\t).

### Table 5: Major cemeteries of London, without and with (r) compensation of population growth. beta – Gompertz beta parameter; M – modal age; ex20 – life expectancy at age 20; ex25 – life expectancy at age 25. Ranges computed with credible HDIs of 95%.

```{r wellcome_overview_all, warning=FALSE}
kable(wellcome_overview_all) %>%
  kableExtra::kable_styling(latex_options = c("HOLD_position","scale_down"))
write.table(wellcome_overview_all, file = "./documented/table05_osteological_estimates.txt", 
            sep="\t", quote = FALSE)
```


Data are hard coded in the code. Sources: @ref_221997, 97--103 table 32 (St Marylebone); @ref_242052, 81 (St Marylebone north of Paddington street)

```{r Basic statistics Marylebone, warning=FALSE}
source("./lifetables_processing/Marylebone.R")
kable(Marylebone_ranges, 
      caption = "St Marylebone, corrected with population growth rate of 2.75 per-cent.") %>%
  kableExtra::kable_styling(latex_options = "HOLD_position")
```

The following plot is build in ./lifetables_processing/stbrides_crypt.R within the if-statement on ```runCodeNew``` (s. data limitations above).

### Figure 11: St. Bride's Crypt. Density of actual ages and Bayesian model of Gompertz distribution of actual ages and osteological estimates (without correction for population growth).

```{r Fig. 11. St. Bride Crypt, fig.height=3, fig.width=6, warning=FALSE}
plot(stbrides_crypt_plot)
```

\newpage
# Supporting information

The chapter 'Supporting information' provides details about the London cemeteries included in the study, the Gompertz parameters of the Coale & Demeny life tables, and the simulations and their results.

## The Coale & Demeny life tables

Calculation of the lowest $\beta$-value for any of the Coale & Demeny life tables (@ref_22426) which is 0.0391 (the female table "West", level 1).

```{r Coale and Demeny life tables beta-value, warning=FALSE, include = TRUE}
source("./chapter_supplement/coale_demeny_life_tables_gompertz.R")
min(gompertz_df$Gompertz_shape)
```

## Simulations

### S.Fig 1: Bayesian model of simulated data with osteological age categories. Difference of estimated to original Gompertz $\beta$ in relation to original $\beta$ (left) and sample size (right).

```{r S.Fig. 1. plot Bayesian model, fig.height=2.5, fig.width=5,warning=FALSE}
#source("./chapter_supplement/simulations_run.R")
gridExtra::grid.arrange(grobs = plot_list_bayes_diff, ncol = 2)
```

\newpage

### S.T.1. Bayesian model with simulated data-set to compare the impact of thinning and additional steps. n = 500, Gompertz $\beta$ = 0.05.

```{r S.T.1. Bayesian modeling stable means}
source("./chapter_supplement/bayes_complete.R") # can take a few minutes
kable (bayes_complete) %>%
  kableExtra::kable_styling(latex_options = c("HOLD_position","scale_down"))
```

\newpage
# References