Why are population growth rate estimates of past and present hunter-gatherers so different?
- Data analyses

Table of Contents

1 Introduction

This document makes it possible to reproduce data manipulations and analyses in Tallavaara and Jørgensen (2020: 'Why are population growth rate estimates of past and present hunter-gatherers so different?' in Philosophical transactions of the Royal Society B) by running the following R-code (https://www.r-project.org/). It assumes that dataTJ2020.Rdata -file is placed in your working directory. This file contains all the required data that are not created within this code. The file can be downloaded from the same Zenodo repository, where this document was downloaded. This document contains citations to references that are listed in the original paper. The pure R-code is in the TJ2020.R-file available in the Zenodo repository.

2 Download R-packages

Install following packages, if they are not already installed to your computer:

install.packages(c("rcarbon", "truncnorm", "WaveletComp", "foreach", "doParallel"),

Next, load these packages into workspace:


3 Data

First, load the dataTJ2020.RData file into workspace. It is assumed that this file is in your working directory:


3.1 Simulated hunter-gatherer dynamics

Next, load Belovsky's simulation results. These data are digitised from the original publication (Belovsky 1988). Each individual dataset contains two columns: time (years) and density (#/100km2).

belo100 <- dataTJ2020[[1]]
belo200 <- dataTJ2020[[2]]
belo400 <- dataTJ2020[[3]]
belo800 <- dataTJ2020[[4]]

Then, different temporal patterns based on different combinations of Belovsky's simulation. These are used as templates for the simulated SPDs. The aim is to imitate long-term hunter-gatherer population dynamics in changing environment. Each individual data is vector of population densities.

tpattern1 <- dataTJ2020[[5]]
tpattern2 <- dataTJ2020[[6]]
tpattern3 <- dataTJ2020[[7]]
tpattern4 <- dataTJ2020[[8]]

3.2 Historical Sámi data

Data from N. Norway cover the n. of tax payers and derived population sizes for two Saami communities, Guovdageaidnu (Kautokeino) and Ávjovárri (Karasjok). The data are from Hood 2015. The data have four colums: year; kautok (n. of taxpayers in Kautokeino community); avjo (n. of taxpayers in Kautokeino community); kautoPop (estimated population size based on n. of taxpayers); and avjoPop (estimated population size based on n. of taxpayers). See Hood 2015 for more details.

saaminor <- dataTJ2020[[9]]

## Calculate total population size
saaminor$sum <- rowSums(saaminor[, 4:5], na.rm=TRUE)
saaminor$sum[saaminor$sum==0] <- NA 

Data from N. Finland cover the n. of taxpayers from the historical Kemi Lappi region, which includes nine historical Saami communities. The data are extracted from Vahtola 2003 and Enbuske 2008. The data have 11 columns: year; Inari (n. of taxpayers (n.o.t.) in Ánar or Inari community); Sodankyla (n.o.t. in Sodankylä community); Kittila (n.o.t in Kittilä); Sompio (n.o.t. in Sompio); Kemikyla (n.o.t. in Kemikylä); Kitkajarvi (n.o.t. in Kitkajärvi); Maanselka (n.o.t. in Maanselkä); Kuolajarvi (n.o.t in Kuolajärvi); Peltojarvi (n.o.t. in Peltojärvi); refs (reference).

## N. of taxpayers in historica Kemi Lappi region in N. Finland.
## Peltojärvi Sámi community is excluded.
saamifin <- dataTJ2020[[10]]
saamifin$sum <- rowSums(saamifin[,2:9], na.rm=TRUE)
saamifin$sum[saamifin$sum==0] <- NA 

3.3 Global human population growth

Global human population growth over times (for the figure 3 in Tallavaara and Jørgensen (2020)). These data are digitised from the Biraben 1979. The data have two columns: age and popsize (human population size).

wpop <- dataTJ2020[[11]]

## Let's make it a little shorter
wpop <- wpop[wpop$age < 30000,] 

3.4 Ethnographic population growth rate estimates

3.4.1 !Kung

Growth rate of Dobe !Kung is from Howell 2000 (1979) page 213-220. See esp. the table 11.1. It is the intrinsic rate of increase, which is calculated using net reproduction rate. Growth rate is 0.26%.

3.4.2 Agta

Agta growth rate is based on the change in population size between 1950-1965. During this period Agta were still foragers. Growth rate is based on population sizes in the table 6.1 in Early and Headland 1998 (page 83). Growth rate is 1.4%.

3.4.3 Asmat

Asmat growth rate is from Van Arsdale 1978 (page 457). Annual growth of 1.5% is based on the change in population size between 1956 and 1973. Crude rate of natural population increase would be 1.9%.

3.4.4 Hadza

Hadza growth rates is from Blurton Jones 2016. Different methods give slightly different estimates, but Blurton Jones considers annual growth of 1.6% to be the most reliable one. See e.g. pages 171-174, 181, and 188 in Blurton Jones 2016.

3.4.5 Ache

Growth rate of Ache is from Hill and Hurtado 1996 page 84. Growth rate of 2% is based on the change in population size between 1930 (241) and 1970 (544), i.e., during the forest phase.

3.5 Archaeological population growth rate estimates

3.5.1 Pleistocene and Holocene Australia

Growth rate in Williams 2013 is rather complex one. Williams fits a smoothing spline to archaeological SPD and then determines the growth rate at every time point in the curve. Then he averages these negative and positive rates and gets an averaged annual growth rate of 0.01%. However, we will use here maximum positive rate estimated from the figure 4 in Williams 2013 (taphonomically corrected). This growth rate is approximately 0.045% annually at 20,000 and 2000 years ago in Australia.

3.5.2 Holocene Australia

Growth rate in Johnson and Brook 2011 is based on the temporal distribution of radiocarbon dates from rock shelters for the last 5000 years. It does give an annual growth rate of 0.04%.

3.5.3 Wyoming and colorado

Growth rate in Zahid et al. 2016 is calculated by fitting an exponential model to the SPD covering the period between 13,000 and 6000 years ago in Wyoming and Colorado. This growth rate is 0.04%.

3.5.4 South America

Growth rate in Goldberg et al. 2016 is calculated by fitting a logistic model to SPD covering the hunter-gatherer period between 14,000 and 6000 years ago in South America. We use here the maximum growth rate of 0.132%.

3.5.5 Kuril Islands

Growth rate in Fitzhugh et al. 2016 is calculated from the SPD between 2500 and 2000 years ago, i.e.,over 500 years. Growth rate is 0.2%.

3.5.6 Big Horn Basin in Wyoming

Growth rate in Kelly et al. 2013 is calculated from smoothed SPD. Mean annual growth rates for five marjor phases of population growth range between 0.16% and 0.31% in Wyoming. We use the largest maximum growth rate.

3.6 Growth rates summarised

This presents ethnographic, archaeological and historical growth rates in tabular form.

Table 1: Ethnographic population growth rate estimates
Group Growth rate (%) Reference
Dobe !Kung 0.26 Howell, N., 2000
Agta 1.8 Early, J.D., Headland, T.N., 1998
Asmat 1.9 Van Arsdale, P.W., 1978.
Hadza 1.6 Blurton Jones, N.G., 2016.
Ache 2 Hill, K., Hurtado, A.M., 1996.
Table 2: Archaeological population growth rate estimates
Region Growth rate (%) Time range Reference
Australia 0.045 40,000 - 0 Williams, A.N., 2013
Australia 0.04 5000 - 0 Johnson, C.N., Brook, B.W., 2011
Wyoming and Colorado 0.04 13,000 - 6000 Zahid, H.J. et al., 2016
South America 0.132 14,000 - 6000 Goldberg, A., et al. 2016
Kuril Islands 0.2 2500 - 2000 Fitzhugh, B., et al. 2016
Big Horn Basin (WY) 0.31 Several periods Kelly, R.L., 2013

Growth rate (r) for historical Sámi populations in N. Norway (Kautokeino and Karasjok communities) and in N. Finland (Kemi Lappi region) are calculated as changes in population size through time during selected growth periods using the following equation:

\begin{equation} r = \ln(N_t/N_0)/t \end{equation}

where N0 is the size of the population at the start of growth and Nt is the size of the population at time t later. t is calculated as a difference between End date and Start date of the period over which the growth rate is calculated.

Table 3: Historical Sámi population growth rates
Region Start date End date t N0 Nt r (%)
N. Norway 1559 1602 43 137 302 1.84
N. Norway 1679 1752 73 106 523 2.19
N. Finland 1555 1605 50 72 111 0.87
N. Finland 1620 1645 25 85 188 3.18
N. Finland 1655 1692 37 86 229 2.65

4 Figure 1

Following code produces the figure 1 in Tallavaara and Jørgensen (2020). This figure compares population growth rates and dynamics between ethnographic, historical, and archaeological data.

4.1 Transparent colors

First, we build a transparent color -function by Mark Gardener (https://www.dataanalytics.org.uk/make-transparent-colors-in-r/)

  t_col <- function(color, percent = 50, name = NULL) {

#         color = color name
#       percent = % transparency
#          name = an optional name for the color

## Get RGB values for named color
  rgb.val <- col2rgb(color)

## Make new color using input color as base and alpha set by transparency
  t.col <- rgb(rgb.val[1], rgb.val[2], rgb.val[3],
               max = 255,
               alpha = (100-percent)*255/100,
               names = name)

## Save the color

## END

4.2 Plotting the figure

This code produces the actual multipanel figure.

## Population growth rates in ethnographic data
par(fig=c(0,0.35,0.60,1), mar=c(3,4,1,3))
plot(NA, xlim=c(0.5,2.5), ylim=c(0.01, 10), xlab="Index", ylab="",
     log="y", yaxt="n", xaxt="n")
text(-0.6, 0.3, "Population growth rate (%)", xpd=NA, srt=90)
at.y <- outer(1:9, 10^(-2:1))
lab.y <- ifelse(log10(at.y) %% 1 == 0, at.y, NA)
abline(h=at.y, col="lightgrey", lty=2)
axis(2, at=at.y, labels=lab.y, las=1, tcl=-0.3, mgp=c(3, 0.5,0))
text(-1.1, 17, "(a)", xpd=NA)

## Population growth rates in archaeological data
stripchart(ethnogr[,2], method="jitter", jitter=0.1, vertical=TRUE, pch=21,
           cex=1.5, bg=t_col("#FAAB18",40), col="black", add=TRUE, at=1)
text(1, 0.002, "ethno", xpd=NA, srt=90)
stripchart(archgr[,2], method="jitter", jitter=0.15, vertical=TRUE, pch=21,
           cex=1.5, bg=t_col("#588300", 40), col="black", add=TRUE, at=2)
text(2, 0.002, "arch", xpd=NA, srt=90)

## Simulated H-G population dynamics according to Belovsky (1988) 
par(fig=c(0.33,0.66,0.79,0.99), mar=c(1,2,0,1), new=TRUE)
plot(belo100, type="l", las=1, axes=FALSE, xlab="", ylab="", ylim=c(0,4))
axis(side=2, at=seq(0,4,1), las=2, tcl=-0.3, mgp=c(3, 0.5,0))
abline(h=mean(belo100$density, na.rm=TRUE))
text(250,2, "0.8%", cex=0.7)
text(150, 3.8, "Productivity 100 g/m2", cex=0.7, font=2)
text(-90, 0, "Population density (#/100km2)", xpd=NA, srt=90)
text(-136, 4, "(b)", xpd=NA)

par(fig=c(0.66,0.99,0.79,0.99), mar=c(1,2,0,1), new=TRUE)
plot(belo200, type="l", las=1, axes=FALSE, xlab="", ylab="", ylim=c(0,15))
axis(side=2, at=seq(0,15,5), las=2, tcl=-0.3, mgp=c(3, 0.5,0))
abline(h=mean(belo200$density[belo200$time > 70], na.rm=TRUE))
text(250,4, "1.2%", cex=0.7)
text(150, 14, "Productivity 200 g/m2", cex=0.7, font=2)

par(fig=c(0.33,0.66,0.6,0.8), mar=c(1,2,0,1), new=TRUE)
plot(belo400, type="l", las=1, axes=FALSE,  xlab="", ylab="", ylim=c(0,70))
axis(side=2, at=seq(0,60,20), las=2, tcl=-0.3, mgp=c(3, 0.5,0))
axis(side=1, at=seq(0,300, 50), tcl=-0.3, mgp=c(3, 0.5,0))
abline(h=mean(belo400$density[belo400$time > 70], na.rm=TRUE))
text(75,10, "3%", cex=0.7)
text(200,10,"0.4%", cex=0.7)
text(150, 66, "Productivity 400 g/m2", cex=0.7, font=2)
text(350, -30, "Time in years", xpd=NA)

par(fig=c(0.66,0.99,0.6,0.8), mar=c(1,2,0,1), new=TRUE)
plot(belo800, type="l", las=1, axes=FALSE, xlab="", ylab="", ylim=c(0,200))
axis(side=2, at=seq(0,160,40), las=2, tcl=-0.3, mgp=c(3, 0.5,0))
axis(side=1, at=seq(0,300, 50), tcl=-0.3, mgp=c(3, 0.5,0))
abline(h=mean(belo800$density[belo800$time > 90], na.rm=TRUE))
text(200,10,"3%", cex=0.7)
text(150, 190, "Productivity 800 g/m2", cex=0.7, font=2)

## Changes in Sámi population size in Northern Norway
par(fig=c(0,1,0.25,0.55),mar=c(3,5,1,1), new=TRUE)
plot(saaminor$year, saaminor$sum, type="n", lwd=1.7, las=1,xaxt="n", yaxt="n",
     xlim=c(1550,1750), xlab="", ylab="Population size")
abline(h=seq(50,500,50), v=seq(1550,1750,10), col="darkgrey", lty=3)
lines(saaminor$year, saaminor$sum, lwd=2, col="#990000")
axis(side=2, at=seq(50,500,50), las=2, tcl=-0.3)
text(1580, 300, paste(samigr[1, 7],"%"), cex=0.7)
text(1700, 350, paste(samigr[2, 7],"%"), cex=0.7)
text(1570,450, "N. Norway", cex=1)
text(1507, 520, "(c)", xpd=NA)

## Changes in the number of Sámi taxpayers in Northern Finland
par(fig=c(0,1,0.05,0.35), new=TRUE)
plot(saamifin$year, saamifin$sum, type="n", lwd=1.7, las=1, xlim=c(1550,1750),yaxt="n",
     xaxt="n", xlab="", ylab="N. of taxpayers")
mtext("Years (AD)", side=1, line=2, cex=1)
abline(h=seq(80,220,20), v=seq(1550,1750,10), col="darkgrey", lty=3)
lines(saamifin$year, saamifin$sum, type="l", cex=1, pch=19, lwd=2, col="#1380A1")
axis(side=1, at=seq(1550,1750,10), lab=rep("",21), tcl=-0.3)
axis(side=2, at=seq(80,220,20), las=2, tcl=-0.3)
axis(side=1, at=seq(1550,1750,50), tcl=-0.4)
axis(side=1, at=seq(1550,1750,10), lab=rep("",21), tcl=-0.3)
text(1580, 120, paste(samigr[3, 7],"%"), cex=0.7)
text(1620, 150, paste(samigr[4, 7],"%"), cex=0.7)
text(1670, 185, paste(samigr[5, 7],"%"), cex=0.7)
text(1570,200, "N. Finland", cex=1)
text(1507, 230, "(d)", xpd=NA)


Figure 1: Archaeological hunter-gatherer population growth rates as compared to ethnographic, simulated and historical estimates. (a) Ethnographic and archaeological estimates of annual change in population size (see table 1). (b) Belovsky’s simulation of hunter-gatherer population dynamics in environments with different productivity. Horizontal lines indicate mean population density. (c) Sámi population size in Guovdageaidnu and Ávjovárri communities in northern Norway. (d) Number of tax payers in the Kemi Lappi region in northern Finland. This is assumed to reflect the Sámi population size in the region. (b–d) Show annual population growth rates (%) during growth periods.

5 Figure 2

5.1 Simulating archaeological SPDs

Here, we simulate archaeological SPDs based on Belovsky's simulated hunter-gatherer dynamics (see tpatterns above). We sample archaeological radiocarbon dates assuming that their temporal distribution follows simulated hunter-gatherer population dynamics. Then, we create SPDs based on these dates. Following code uses functions from rcarbon packages (as is or slightly modified versions). Errors of the simulated dates are sampled from the truncated normal distributions (rtruncnorm).

The simulation uses foreach -loop from foreach package. It assumes that you have a multi-core setup. If your computer has a single core, you should switch dopar to do in the code or change the loop into traditional for -loop.

The results of the first simulation are used in fig 2b and supplementary figures 2 and 3.

## Parameter values for the simulation
startyr <- 10000
yearRange1 <- c(startyr,startyr - length(tpattern1)+1)
ndates1 <- length(yearRange1[1]:yearRange1[2])
dnorm <- TRUE
dnormspd <- TRUE
spdnorm <- TRUE
pweights1 <- tpattern1/sum(tpattern1) ## Temporal distribution follows tpattern1     
samplesize <- 5000 ## Number of calendar dates in each sample
error <- round(rtruncnorm(samplesize, a=20, b=80, mean = 50, sd = 15),0) 
nsim <- 9 ## Number of samples/simulation rounds
sim1 <- matrix(NA,nrow=ndates1,ncol=nsim)
n_cores <- detectCores()/2 ## parallel processing will use half of the available cores

The simulation itself starts here. Parallel execution using 4 cores on i7 machine would usually take less than 5min.

cl <- makeCluster(n_cores)

sim1 <- foreach(i = 1:nsim,
               export=c(yearRange1, samplesize, pweights1, error, dnorm, dnormspd, spdnorm),
               .packages="rcarbon", .combine=cbind) %dopar% {
                                              replace=TRUE, prob=pweights1))
                 calibdates1 <- calibrate(x=randomDates1.2, errors=randomSDs1,
                                          calCurves='intcal13', normalised=dnorm,
                                          ncores = n_cores, verbose = FALSE)
                 spd(calibdates1, timeRange=yearRange1, datenormalised = dnormspd,
                     spdnormalised = spdnorm, verbose = FALSE)$grid$PrDens


Make a data frame from the resulting matrix and add time variable date and the underlying "true" temporal pattern tpattern to the data.

res <- as.data.frame(sim1) # matrix as data frame
res$date <- yearRange1[1]:yearRange1[2] # add time variable
res$tpattern <- tpattern1 # add the underlying temporal pattern to the dataset

The results of the second simulation are used in fig 2c-d. Here, we produce only one round of simulations, so no loop here. This simulation uses also slightly different parameters.

This simulation and parameters are for the fig 2c

yearRange2 <- c(startyr,startyr - length(tpattern3)+1)
pweights2 <- tpattern3/sum(tpattern3)       

Simulation starts here:

calibdates2 <- calibrate(x=randomDates2.2, errors=randomSDs2, timeRange=yearRange2,
                         calCurves='intcal13', normalised=dnorm, ncore=n_cores, verbose=FALSE)
sim2 <- as.data.frame(spd(calibdates2, timeRange=yearRange2, datenormalised = dnormspd,
                          spdnormalised = spdnorm, verbose=FALSE)$grid)

This simulation and parameters are for the fig 2d

yearRange3 <- c(startyr,startyr - length(tpattern2)+1)
pweights3 <- tpattern2/sum(tpattern2)       

Simulation starts here:

                           replace=TRUE, prob=pweights3))
calibdates3 <- calibrate(x=randomDates3.2, errors=randomSDs3, timeRange=yearRange3,
                        calCurves='intcal13', normalised=dnorm, ncores=n_cores, verbose=FALSE)
sim3 <- as.data.frame(spd(calibdates3, timeRange=yearRange3, datenormalised = dnormspd,
                          spdnormalised = spdnorm, verbose=FALSE)$grid)

5.2 Long-term growth rates

Long-term growth rates are calculated by fitting exponential models to SPDs. The (exponential) growth rate is the coefficient for the time variable.

This is for the figure 2b

## Exponential fit
## Add new time variable "yr" to the simulated data just for calculating the exp fit.
## This is only to get positive coefficient
res$yr <- 1:4120 

## Fitting the exponential model to the firts sample
## of the simulated data
expfit1 <- lm(log(res[,1]) ~ yr, data=res)
grSIM1 <- summary(expfit1)$coefficients[[2]]

# Predicted 
exprob1 <- exp(predict(expfit1))

This is the rate of change in environmental productivity

## "Growth rate" for environmental productivity 
n0 <- 100 # productivity at the beginning
nT <- 800 # productivity in the end
t <- 3000 # time, 9500-6500 yrs ago
grEP <- log(nT/n0)/t

These are the growth rates for the figure 2c-d

## Exponential fits for the fig 2c-d

## Adding time variables just to get calculate positive growth rates
sim2$yr <- 1:nrow(sim2) 
sim3$yr <- 1:nrow(sim3)

fitSim2 <- lm(log(PrDens) ~ yr, data=sim2)
grSIM2 <-  summary(fitSim2)$coefficients[[2]]
exprobSim2 <- exp(predict(fitSim2))

fitSim3 <- lm(log(PrDens) ~ yr, data=sim3)
grSIM3 <- summary(fitSim3)$coefficients[[2]]
exprobSim3 <- exp(predict(fitSim3))

5.3 Plotting the figure 2

This code produces the actual 4 by 4 multipanel figure.

par(mar=c(3,5,1,5), mfcol=c(2,2), cex.lab=1)
##par(fig=c(0,1,0.45,0.95), mar=c(3,5,1,1))
plot(res$date*-1, res$tpattern, type="l", xlim=c(-10000,-5800),
     ylab="Population density (#/100km2)", las=1, xaxt="n", yaxt="n")
text(-10000, 130, "(a)")
axis(1, at=seq(-10000,-6000, by=1000), lab=seq(10000,6000, by=-1000), tcl=-0.3)
axis(2, at=seq(0,140, by=20), las=2, hadj=0.7, tcl=-0.3)

##par(fig=c(0,1,0.02,0.52),new=TRUE, mar=c(3,5,1,1))
plot(res$date*-1, res[, 1], type="l", xlim=c(-10000,-5800),
     xaxt="n", yaxt="n", xlab="", ylab="Summed probability")
text(-10000, 0.00058, "(b)")
mtext( "Years ago", side=1, line=2)
axis(1, at=seq(-10000,-6000, by=1000), lab=seq(10000,6000, by=-1000), tcl=-0.3)
axis(2, at=seq(0e-00,8e-04, by=2e-04), lab=c("0.0000","0.0002","0.0004","0.0006", "0.0008"),
     las=2, hadj=0.7,  tcl=-0.3)
lines(res$date*-1, exprob1, col="red")
text(-9000, 0.0005, paste("Growth rate=", round(grSIM1*100, 2),"%"))
text(-8600, 0.00045, paste("Rate of ghange in env. prod=", round(grEP*100, 2),"%"))

plot(sim2$calBP*-1, tpattern3, type="l", col="darkgrey",
     xlim=c(-10000,-7500), xaxt="n", yaxt="n", xlab="",
     ylab="Population density (#/100km2)", lwd=1.5)
text(-10000, 130, "(c)")
mtext("Summed probability", side=4, line=4, cex=1)
axis(side=2, at=seq(0,140,20), las=2, tcl=-0.3)
axis(side=1, at=seq(-10000,-7500, 500), lab=seq(10000,7500,-500),
     tcl=- 0.3)#, lab=c("10,000","9,500", "9,000", "8,500", "8,000", "7,500"))
plot(sim2$calBP*-1, sim2$PrDens, type="l",
     xlim=c(-10000,-7500), lwd=2, xaxt="n", yaxt="n", xlab="",
     ylab="", axes=FALSE)
axis(side=4, at=seq(0.0001,0.0004, 0.0001),
     lab=c("0.0001","0.0002","0.0003","0.0004"), las=2, tcl=-0.3)
text(-9000, 0.00015, paste("Growth rate=", round(grSIM2*100, 2), "%"))
lines(sim2$calBP*-1, exprobSim2, lwd=2, col="red")

plot(sim3$calBP*-1, tpattern2, type="l", col="darkgrey",
     xlim=c(-10000,-7500), yaxt="n", xaxt="n", xlab="",
     ylab="Population density (#/100km2)", lwd=1.5)
text(-10000, 130, "(d)")
mtext("Summed probability", side=4, line=4)
mtext("Years ago", side=1, line=2)
axis(side=2, at=seq(0,140,20), las=2, tcl=-0.3)
axis(side=1, at=seq(-10000,-7500, 500), lab=seq(10000,7500,-500), tcl=-0.3)
plot(sim3$calBP*-1, sim3$PrDens, type="l",
     xlim=c(-10000,-7500), lwd=2, xaxt="n",yaxt="n", xlab="", ylab="", axes=FALSE)
axis(side=4, at=seq(0,0.0012, 0.0002),las=2, tcl=-0.3)
text(-9200, 0.0003, paste("Growth rate=", round(grSIM3*100, 2),"%"))
lines(sim3$calBP*-1, exprobSim3, lwd=2, col="red")


Figure 2: Simulated archaeological proxies of hunter-gatherer population dynamics. (a) Simulated hunter-gatherer population trajectory based on Belovsky´s simulation. Different regimes of environmental productivity are shown below the curve. These regimes correspond to different long-term mean population densities. (b) Archaeological proxy (SPD) of the population pattern shown in (a). The proxy is based on dates sampled under the assumption that their distribution follows the underlying pattern shown in (a). (c—d) Underlying population patterns (grey curves) and their archaeological proxies (black curves) in environments with constant (c) and changing (d) productivity. Productivity regimes are shown above (c) and below (d) the curves. (b—d) Show exponential growth models and annual growth rates (%).

6 Figure 3

This figure illustrates differences between different scales of analysis/dynamics - global (partly "imaginary"), regional (archaeological proxies), and actual dynamics (not captured by archaeological proxies). Regional archaeological proxy in this figure is, again, based on simulated data. The simulation code is shown first.

## This is the regional long-term scale, i.e., archaeological proxy
## Requires some more simulation using yet another underlying pattern
## (tpattern4)

startyr2 <- 17000
yearRange4 <- c(startyr2,startyr2 - length(tpattern4)+1)
pweights4 <- tpattern4/sum(tpattern4)       

                           replace=TRUE, prob=pweights4))       

calibdates4 <- calibrate(x=randomDates4.2, errors=randomSDs4, timeRange=yearRange4,
                         calCurves='intcal13', normalised=dnorm,
                         ncores=n_cores, verbose=FALSE)
simpal <- as.data.frame(spd(calibdates4, timeRange=yearRange4,
                            datenormalised = dnormspd, spdnormalised = spdnorm,

Following code produces the figure.

par(mfrow=c(3,1), mar=c(3,5,1,1))
## Global scale population growth based on Biraben 1979
plot(wpop$age*-1, wpop$popsize, type="l", xlim=c(-30000,0), axes=F,
     xlab="", ylab="Global population size (billions)")
axis(side=1, at=seq(-30000,0, 10000), lab=seq(30000,0,-10000),
     pos=-0.3e+09, tcl=-0.3, mgp=c(3,0.5,0))
axis(side=1, at=c(seq(-25000,-5000, 10000)), lab=rep("",3),
     tck=-0.03, pos=-0.3e+09, tcl=-0.2)
axis(side=2, at=seq(0,8e+09, 2e+09), lab=seq(0,8,2), las=2, tcl=-0.3)
axis(side=2, at=seq(1e+09,7e+09, 2e+09), lab=rep("",4), las=2, tcl=-0.2)
text(-30000,  7.5e+09, "(a)")

## Plotting the regional long-term scale
plot(simpal$calBP*-1, simpal$PrDens, type="l", axes=FALSE,
     ylab="Summed probability", xlab="", xlim=c(-17000,-11000))
## mtext("Summed probability", side=2, line=3, cex=1)
axis(1, at=seq(-17000,-11000, by=1000), lab=seq(17000,11000, by=-1000),
     tcl=-0.3, mgp=c(3,0.5,0))
axis(1, at=seq(-16500,-11500, by=1000), lab=rep("",6), tcl=-0.2, mgp=c(3,0.5,0))
axis(2, at=seq(0e-00,4e-04, by=1e-04), las=2, hadj=0.7, tcl=-0.3)
text(-17000,  3.5e-04, "(b)")

## And finally te "true" population dynamics
plot(simpal$calBP*-1, tpattern4, type="l", axes=F, ylab="Population density (#/100km2)",
     xlab="", xlim=c(-14700,-13200), ylim=c(0,150))
mtext("Years ago", side=1, line=2)
## mtext("Population density (#/100km2)", side=2, line=3, cex=1)
axis(1, at=seq(-14700,-13200, by=300), lab=seq(14700,13200, by=-300),
     tcl=-0.3, mgp=c(3,0.5,0))
axis(1, at=seq(-14550,-13350, by=300), lab=rep("",5), tcl=-0.2, mgp=c(3,0.5,0))
axis(2, at=seq(0,140, by=10), las=2, hadj=0.7, tcl=-0.3)
text(-14700,  125, "(c)") 


Figure 3: Figure 3. Schematic representation of different scales of analyses and population dynamics of hunter-gatherers. (a) Global scale trajectory of human population size [47] implies that hunter-gatherer populations have been stationary. Prehistoric sections of such reconstructions are usually assumed rather than inferred from data. (b) Continental-to-regional scale, which is trackable with archaeological proxies, suggest that long-term population dynamics are characterised by periods of growth and decline. (c) The scale of actual population size is maybe the most dynamic, but it is usually beyond the resolution of archaeological methods.

7 Supplementary figure 1

This figure shows the change in the number of taxpayers through time in the historical Kemi Lappi region in northern Finland. See Figure 1 for the combined curve.

## Plotting individual Saami communities (excluding Peltojärvi)

## Vector of individual communities in 'correct order'
communities <- c("Inari", "Kittila", "Sodankyla", "Sompio", "Kemikyla", "Kuolajarvi",
                 "Kitkajarvi", "Maanselka")
commNames <- c("Inari", "Kittiä", "Sodankylä", "Sompio", "Kemikylä", "Kuolajärvi",
               "Kitkajärvi", "Maanselkä")

## For loop to create a multipanel plot
par(mfrow=c(4,2), mar=c(4,4,1,1))
for(i in 1:length(communities)){
  plot(saamifin$year, saamifin[,communities[i]], type="l", xlab="", ylab="", col="black",
       lwd=1, las=1, xaxt="n", tcl=-0.3, main=commNames[i])
  axis(side=1, at=seq(1550, 1700,50), tcl=-0.3)
  axis(side=1, at=seq(1550, 1700, 10), lab=rep("", length(seq(1550, 1700,10))), tcl=-0.2)
text(1530, -10, "Year AD", xpd=NA)
text(1330, 75, "N. of taxpayers", xpd=NA, srt=90)


Figure 4: (Supplementary figure 1) Number of taxpayers in individual Sámi communities in the historical Kemi Lappi region in northern Finland.

8 Supplementary figure 2

This figure shows all the nine samples of simulated SPDs created in the section 5. In addition, it shows the underlying population dynamics (on the top panel of the figure).

## First, a layout with the top panel having only the
## underlying simulated population dynamics (tpattern)
                8,9,10), 4,3, byrow = TRUE))

## Plotting the underlying dynamics
plot(res$date*-1, res$tpattern, type="l", xlim=c(-10000,-5800),
xaxt="n", xlab="", ylab="Population density (#/100km2)", las=2)
axis(1, at=seq(-10000,-6000, by=1000), lab=seq(10000,6000, by=-1000), tck=-0.04, mgp=c(3,0.2,0))

## Plotting the rest using for loop
for(i in 1:9){
  plot(res$date*-1, res[,i], type="l", ylim=c(0,0.0008), xlim=c(-10000,-5800),
       xaxt="n", yaxt="n", xlab="", ylab="")
  axis(1, at=seq(-10000,-6000, by=1000), lab=seq(10000,6000, by=-1000),
       tck=-0.04, mgp=c(3,0.2,0))
  axis(2, at=seq(0e-00,8e-04, by=2e-04), las=2, hadj=0.7, tcl=-0.2)

## Adding axis labels
text(-13800, -2.2e-04, "Years ago", xpd=NA)
text(-23000, 2.2e-03, "Summed probability", xpd=NA, srt=90)


Figure 5: (Supplementary figure 2) SPDs based on nine different random samples (N=5000) of calendar dates from the population pattern shown on the top row of the figure.

9 Supplementary figure 3

This figure shows the wavelet coherence between the underlying population dynamics and archaeological SPDs. It uses the simulated data created in the section 5. First part of the code runs a foreach -loop that compares each archaeological SPD to the underlying dynamics. Parallel execution of the loop using 4 cores on i7 machine would take c. 10min.

cl <- makeCluster(n_cores)

## strt<-Sys.time()
    wc_list <- foreach(i = 1:9, .packages="WaveletComp") %dopar% {
        analyze.coherency(res, my.pair=c(i, 11),
                          loess.span=1, n.sim=100, window.type.t=3, window.type.s=3,
                          window.size.t=1, window.size.s=1, verbose=FALSE)


## print(Sys.time()-strt)

The second part of the code produces the multipanel figure.

# Plotting the wavelet coherence

par(mfrow=c(3,3), mar=c(4,5,1,5))
for(i in 1:9){
    wc.image(wc_list[[i]], which.image="wc", color.key="interval",
         spec.time.axis=list(at=c(1,1000,2000,3000,4000), labels=c(10000,9000,8000,7000,6000)),
         periodlab="", plot.ridge=FALSE, graphics.reset=FALSE,
         legend.params = list(width=1.2, shrink = 0.9, mar = 12.1,
                              n.ticks = 6, label.digits = 1,
                              label.format = "f",lab = NULL,
                              lab.line = 1.5))

text(-17200, 19, "Period", xpd=NA, srt=90, cex=1.5)
text(5900, 19, "Wavelet coherence level", xpd=NA, srt=-90, cex=1.5)
text(-4300, -1, "Years ago", xpd=NA, cex=1.5)


Figure 6: (Supplementary figure 3) Wavelet coherence plots showing the similarity between sample SPDs and the underlying actual population pattern. Significant similarities in the periodicity are indicated by the white contour and red colour. Horizontal arrows pointing to the right indicate that the two series (actual pattern and SPD) are in phase at the respective period. Arrows pointing to the left indicate that the two series are in anti-phase.

Author: Miikka Tallavaara

Created: 2020-03-31 ti 09:47