esm.Rmd

---
title: "How cultural transmission through objects impacts inferences about cultural evolution: supplementary material"
author: "Crema,E.R., Bortolini, E., Lake,M.W."
date: "28 October 2022"
output: html_document
---

```{r, setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
options(scipen=9999)
```

```{r, settings_libraries_functions, include=FALSE}
# Load Required Packages
library(foreach)
library(doParallel)
library(RColorBrewer)
library(here)

#Source Functions
source(here('src.R'))
```


# Models

Consider a population of $N$ individuals each associated with a mental representation of a cultural variant. At each generation step, each individual produce $n$ copies of a physical realisation of the mental representation, with $n$ drawn from a Poisson distribution with intensity $lambda$. After this _production event_ all individuals update their mental representation by randomly selecting an object and copying the associated cultural variant. New variants are introduced into the population in two distinct ways:

## Decoding Error

Under this scenario, with probability $\mu_d$, the focal individual fails to correctly reconstruct the mental representation associated to the object it copies from, and as consequence introduces a new cultural variant in the population of the mental representations. The new variants are introduced following the infinite allele model (i.e. there are 0 probability of convergent error).  

## Encoding Error

With probability $\mu_e$ an object produced by the focal individual is no longer directly associated to its parent mental representation. As for the decoding error scenario described above new variants are introduced following the infinite allele model. However, the mental representation of the focal individual remains unchanged (unless it subsequently select the mutant object in the transmission process) and an individual could produce more than one mutant object, each associated with a different mental representation that does not exist in the population of the producers yet.  

## Model Scheduling and Data Collection

Both models are initiated with $N$ individuals associated with different mental representation represented by an integer value. At each time-step the following processes occur:

1. *Production*. For each agent $i$, $n_i$ identical objects with the same numerical value as the mental representation of the agent are created, with $n_i \sim Poisson(\lambda)$.

2. *Encoding Error* With probability $\mu_e$, each object value is associated with a new integer value that has not been introduced in the simulation run. 

3. *Cultural Transmission*. Each agent randomly selects an object and updates its mental representation number with the integer value associated with the object.

4. *Decoding Error* With probability $\mu_d$, each agent updates its mental representation value with a new number that has not been previously used in the simulation run. This portrays an instance where the individual fails to correctly reconstruct the mental representation of an object.

Models are executed first for $T_{burnin}+T_{collection}$ time-steps. The first $T_{burnin} steps ensure that the simulation reach equilibrium conditions. During $T_{collection}$, the frequency of different cultural variants associated to the objects are retained for analysis. 

Notice that both $mu_e$ and $mu_d$ can be positive, and two errors can co-exist. However, in this case will explore each error type separately, referring to instances of $\mu_e>0,\mu_d=0$ as _encoding error model_ and $\mu_e=0,\mu_d>0$ as _decoding error model_.

### Wright Fisher Model

For comparative purposes we also simulated the outputs of a Wright-Fisher model where individuals directly copy the mental representation and with a probability $\mu$ in the copying process.

# Experiment Design

We compared the output of our three scenarios (wright-fisher, unbiased object-mediated transmission with encoding error, and unbiased object-mediated transmission with decoding error) using different descriptive statistics and visualisation techniques (diversity, richness, and progeny distribution) employed in cultural evolutionary studies.

## Experiment 1a

In experiment 1a we ran our simulations using different combinations of population size $N$ and production rate $\lambda$ for each of the two error type (i.e., encoding and decoding), and recorded the richness $k$ (the number of unique variants) and homogeneity $H$ (equivalent to $\sum_{i}^{k}p_{i}^2$, with $p_i$ being the proportion of the $i$-th variant) for the population of objects and mental representations. We used the following parameter settings, recording the two summary statistics after 5,000 time-steps.

```{r,parameter_setting_exp1}
nsim = 1000 #number of repetitions for each parameter combination
N  = c(100,500,1000) #Population size parameter sweep
L  = c(0.5,1,5,10) #Productivity rate parameter sweep
mu.e =mu.d = 0.01 #error rates
# crate parameter space
pspace.obj = expand.grid(nsim=1:nsim,lambda=L,pop=N) 
timesteps <- 5000 #number of timesteps
```

```{r,register_clusters_exp1,include=FALSE}
#Setup parallel processing using n-1 threads, with n being maximum possible for the specific machine
cl <-  makeCluster(detectCores()-1)  
registerDoParallel(cl)
```

```{r,execution_exp_1}
res.decode = foreach (i = 1:nrow(pspace.obj), .combine='rbind') %dopar% {
unlist(objTr(N=pspace.obj$pop[i],lambda=pspace.obj$lambda[i],mu.e=0,mu.d=mu.d,timesteps = 5000,output="sumstat"))}
pspace.obj.decode = cbind.data.frame(pspace.obj,res.decode)

res.encode = foreach (i = 1:nrow(pspace.obj), .combine='rbind') %dopar% {
unlist(objTr(N=pspace.obj$pop[i],lambda=pspace.obj$lambda[i],mu.e=mu.e,mu.d=0,timesteps = 5000,output="sumstat"))}
pspace.obj.encode = cbind.data.frame(pspace.obj,res.encode)
```

```{r, close_clusters_exp1, include=FALSE}
#stop cluster
stopCluster(cl)
```
For a given population size, richness show an increase whilst homogeneity show a decrease with higher production rates, although the patterns is detectable only with larger differences in $\lambda$. In the case of homogeneity, no differences has been identified between statistics calculated on the frequency of mental representations or objects nor on whether the mutation occurred during encoding or decoding, effectively yielding similar values for a given combination of $N$ and $\lambda$. This was not the case for richness, particular with larger settings of $N$. In particular, when the error occurred during the encoding process, the number of unique variants was significantly higher than the number of mental representations. The difference between the two summary statistics is most likely caused by the fact that, under encoding error regimes, the number of new variants introduced at a given time-steps is on average $N$ times $\lambda$ times $\mu_e$. Most of these variants are, however, not adopted by the social learners and hence go extinct after one generation (see also experiment 3 below), leading to a marked difference between the richness of objects and mental representations. In contrast, under decoding error regimes, the number of new variants is simply given by $N$ times $\mu_d$, with $\lambda$ this time defining the frequencies of these new variants. These settings lead to lower values of richness, as the number of novel variants is $\lambda$ times smaller compared to encoding error regime.   

```{r,exp1a_figure1,results='hide',fig.width= 10, fig.height= 4}
cols = brewer.pal(4, 'Set2')
par(mfrow=c(1,3))
for (i in 1:length(N))
{
	tmp.decode = subset(pspace.obj.decode,pop==N[i])	
	tmp.encode = subset(pspace.obj.encode,pop==N[i])	
	plot(NULL,xlim=c(0.5,4.5),ylim=c(0,1),xlab=expression(lambda),ylab='Homogeneity',axes=FALSE,main=paste0('N=',N[i]))
	for (j in 1:length(L))
	{
		boxplot(at=j-0.32,x=subset(tmp.decode,lambda==L[j])$hom.obj,col=cols[1],add=T,axes=F,boxwex=0.3,outline=FALSE,whisklty=1)
		boxplot(at=j-0.12,x=subset(tmp.decode,lambda==L[j])$hom.mental,col=cols[2],add=T,axes=F,boxwex=0.3,outline=FALSE,whisklty=1)
		boxplot(at=j+0.12,x=subset(tmp.encode,lambda==L[j])$hom.obj,col=cols[3],add=T,axes=F,boxwex=0.3,outline=FALSE,whisklty=1)
		boxplot(at=j+0.32,x=subset(tmp.encode,lambda==L[j])$hom.mental,col=cols[4],add=T,axes=F,boxwex=0.3,outline=FALSE,whisklty=1)
	}
	axis(1,at=c(1,2,3,4),labels=L)
	axis(2)
	box()
}

legend('topleft',legend=c('Decoding Error - Objects','Decoding Error - Mental Representations','Encoding Error - Objects','Encoding Error - Mental Representations'),fill=cols,bty='n')
```

```{r,exp1a_figure2,results='hide',fig.width= 10, fig.height= 4}
cols = brewer.pal(4, 'Set2')
par(mfrow=c(1,3))
maxK = max(c(pspace.obj.decode$k.obj,pspace.obj.decode$k.mental,pspace.obj.encode$k.obj,pspace.obj.encode$k.mental))
for (i in 1:length(N))
{
	tmp.decode = subset(pspace.obj.decode,pop==N[i])	
	tmp.encode = subset(pspace.obj.encode,pop==N[i])	
	plot(NULL,xlim=c(0.5,4.5),ylim=c(0,maxK),xlab=expression(lambda),ylab='Richness',axes=FALSE,main=paste0('N=',N[i]))
	for (j in 1:length(L))
	{
		boxplot(at=j-0.32,x=subset(tmp.decode,lambda==L[j])$k.obj,col=cols[1],add=T,axes=F,boxwex=0.3,outline=FALSE,whisklty=1)
		boxplot(at=j-0.12,x=subset(tmp.decode,lambda==L[j])$k.mental,col=cols[2],add=T,axes=F,boxwex=0.3,outline=FALSE,whisklty=1)
		boxplot(at=j+0.12,x=subset(tmp.encode,lambda==L[j])$k.obj,col=cols[3],add=T,axes=F,boxwex=0.3,outline=FALSE,whisklty=1)
		boxplot(at=j+0.32,x=subset(tmp.encode,lambda==L[j])$k.mental,col=cols[4],add=T,axes=F,boxwex=0.3,outline=FALSE,whisklty=1)
	}
	axis(1,at=c(1,2,3,4),labels=L)
	axis(2)
	box()
	if (i==1)
	{
		legend('topleft',legend=c('Decoding Error - Objects','Decoding Error - Mental Representations','Encoding Error - Objects','Encoding Error - Mental Representations'),fill=cols,bty='n')
	}
}
```


## Experiment 1b

In most real-world situations we do not have estimates of $N$ or $\lambda$, and we simply infer modes of transmission from the frequency distribution of the observed objects. In experiment 1b we emulate this condition by exploring four different combinations of $\lambda$ and $N$, ensuring that the number of objects produced at each time step was, on average, equal to 1,000. We compare the output of these settings against a conventional Wright-Fisher Infinite Allele model with a population size randomly drawn from a Poisson distribution with the same intensity 1,000 to determine the consequence of implicitly assuming the population size to be comparable, on average, to the number of objects. Detailed parameter settings are as follows:

```{r, parameters_exp2}
nsim = 1000 #number of repetitions for each parameter combination
n.objects = 1000 #target number of objects produced at each time-step 
N <- c(50,500,1000,2000) #number of individuals (producers)
lambda <- n.objects/N # production rate
# create parameter spaces
pspace = data.frame(N=rep(N,each=nsim), lambda=rep(lambda,each=nsim)) 
mu.e = mu.d = mu = 0.01 #error rates
timesteps <- 5000 #number of time-steps
```

```{r,register_clusters_exp2,include=FALSE}
#Setup parallel processing using n-1 threads, with n being maximum possible for the specific machine
cl <-  makeCluster(detectCores()-1)  
registerDoParallel(cl)
```

```{r execution_exp2}
res.wf = foreach (i = 1:nsim, .combine='cbind') %dopar% {
  unlist(wf(N=rpois(1,n.objects),mu=mu,timesteps = 5000,output="sumstat"))
}
res.wf.hom = res.wf[1,]
res.wf.k = res.wf[2,]


res.decode = foreach (i = 1:nrow(pspace), .combine='cbind') %dopar% {
  unlist(objTr(N=pspace$N[i],lambda=pspace$lambda[i],mu.e=0,mu.d=mu.d,timesteps = 5000,output="sumstat"))}
res.decode.hom=matrix(res.decode[1,],ncol=length(N),nrow=nsim)
res.decode.k=matrix(res.decode[2,],ncol=length(N),nrow=nsim)

res.encode = foreach (i = 1:nrow(pspace), .combine='cbind') %dopar% {
  unlist(objTr(N=pspace$N[i],lambda=pspace$lambda[i],mu.e=mu.e,mu.d=0,timesteps = 5000,output="sumstat"))}
res.encode.hom=matrix(res.encode[1,],ncol=length(N),nrow=nsim)
res.encode.k=matrix(res.encode[2,],ncol=length(N),nrow=nsim)

```

```{r, close_clusters_exp2, include=FALSE}
#stop cluster
stopCluster(cl)
```


```{r,plot_res_exp2, results='hide',fig.width= 9, fig.height= 9}

par(mfrow=c(2,1))
xs=jitter(rep(c(1,4,7,10,2,5,8,11,14),each=nsim),factor=1.2)

# Diversity
ys.hom=(c(res.decode.hom,res.encode.hom,res.wf.hom))

lo.hom=quantile(res.wf.hom,0.025)
hi.hom=quantile(res.wf.hom,0.975)
col.hom=vector(length=length(ys.hom))

col.hom[ys.hom>=lo.hom&ys.hom<=hi.hom]=rgb(0,0,0,0.05)
col.hom[ys.hom<=lo.hom]=add.alpha("royalblue",0.05)
col.hom[ys.hom>=hi.hom]=add.alpha("indianred",0.05)

plot(xs,ys.hom,col=col.hom,pch=20,axes=F,xlab="",ylab="Homogeneity")
abline(h=c(lo.hom,hi.hom),col=c("royalblue","indianred"),lty=4)
axis(side=1,at=c(1.5,4.5,7.5,10.5,14),labels=c(N,n.objects))
mtext(side=1,line=2.5,"N")
axis(side=2)
axis(side=3,at=c(1.5,4.5,7.5,10.5,14),labels=c(n.objects/N,"NA\n (Wright-Fisher)"))
mtext(side=3,line=2.5,expression(lambda))
box()

text(4,0.8,"Decoding Error",srt=90,cex=0.8)
text(5,0.8,"Encoding Error",srt=90,cex=0.8)

# Richness
ys.k=(c(res.decode.k,res.encode.k,res.wf.k))

lo.k=quantile(res.wf.k,0.025)
hi.k=quantile(res.wf.k,0.975)
col.k=vector(length=length(ys.k))

col.k[ys.k>=lo.k&ys.k<=hi.k]=rgb(0,0,0,0.05)
col.k[ys.k<=lo.k]=add.alpha("royalblue",0.05)
col.k[ys.k>=hi.k]=add.alpha("indianred",0.05)

plot(xs,ys.k,col=col.k,pch=20,axes=F,xlab="",ylab="Richness (k)")
abline(h=c(lo.k,hi.k),col=c("royalblue","indianred"),lty=4)
axis(side=1,at=c(1.5,4.5,7.5,10.5,14),labels=c(N,n.objects))
mtext(side=1,line=2.5,"N")
axis(side=2)
axis(side=3,at=c(1.5,4.5,7.5,10.5,14),labels=c(n.objects/N,"NA\n (Wright-Fisher)"))
mtext(side=3,line=2.5,expression(lambda))
box()

text(1,44,"Decoding Error",srt=90,cex=0.8)
text(2,44,"Encoding Error",srt=90,cex=0.8)
```

Richness showed consistently lower values, whilst homogeneity showed higher values compared to the Wright-Fisher model for both versions of object-mediated transmission, even in situations where $N$ was larger than the value set for the Wright-Fisher model. 
With lower values of $N$ (and consequently in this case higher values of $\lambda$) the number of mental representations is limited and hence there is constraint in the amount of variability in the system. 
Under the decoding error scenario with $N=$ `r N[1]` the highest richness values possible (assuming all individuals have a different mental representation) is `r N[1]` ($k=N$), whilst the lowest homogeneity value is equal to `r N[1]*(1/N[1])^2` ($\left(\frac{\lambda}{\lambda N}\right)^{2}N$, or more simply $\left(\frac{1}{N}\right)^{2}N$). 

The situation is slightly different under encoding error regimes as one has to account, on average, on $\lambda N \mu_e$ newly introduced unique variants. 
Given $\lambda=$`r lambda[1]` and $\mu_e=$`r mu.e`, richness would be capped at approximately to `r round(N[1] + lambda[1]*N[1]*mu.e)`, equivalent to $N+\lambda N\mu_{e}$. 
The lowest limit for homogeneity (in this case equal to  `r N[1]*((1-mu.e)/N[1])^2 + mu.e / (lambda[1]*N[1])`), can be estimated with the equation $\left( \frac{\lambda-\lambda\mu_e}{\lambda N} \right)^2N+\left(\frac{1}{\lambda N}\right)^2\lambda N \mu_e$ (or more simply $\left( \frac{1-\mu_e}{N} \right)^2N+\frac{\mu_e}{\lambda N}$), where the first term is the contribution of objects reflecting the mental representations, and the second terms are the new variants introduced in the pool. 
In contrast to decoding error regimes, where limits in richness and homogeneity depend exclusively on $N$, both $\mu_e$ and $\lambda$ also condition cultural diversity, with lower homogeneity associated with higher error and production rates. 
These differences (or lack thereof) can be observed when we compare the two forms of error for each parameter combination of objected-mediated transmission. 

The higher homogeneity and lower richness of the object-mediated transmission compared to the Wright-Fisher model is, however, not just explained by the combination of low $N$ and high $\lambda$. In two of the parameter combinations, we considered instances where $N$ was equal or higher for the object mediated transmission, yet homogentiy yielded higher and richness yielded lower figures. This is explained by the additional drift caused by stochasticity in the production event. For example with $\lambda=1$, approximately `r round(dpois(0,1)*100,2)`% of the individuals will not produce an object, limiting the amount of richness and diversity that could be potentially realised.


## Experiment 2

Bentley et al [1](https://doi.org/10.1098/rspb.2004.2746) first noted that given a temporal interval $T$, $log(P_k$), the log frequencies of variants appearing $k$ times, follow a power-law distribution under neutrality. It follows that the number of variants occurring once (i.e. $k=1$) within $T$ is far greater than those appearing twice ($k=2$), three times ($k=3$),etc., and that this rate of decline with increasing $k$ can be predicted. O'Dwyler and Kandler[2](https://doi.org/10.1098/rstb.2016.0426) recently examined the shape of such *progeny distribution* and noticed that:

* there is indeed a power law distribution, but with a constant exponent of $-3/2$, which *does not* depend on $N$ and $\mu$;
* the power law is actually followed by an exponential cut-off;
* the cut-off value depends on the innovation rate $\mu$;

This can be demonstrated by the following simulation:


```{r,execution_progeny_wf}
# changing population size
wf.prog.n300.m01 = wf(N=300,mu=0.01,warmup=10000,timesteps=20000,output="progeny")
wf.prog.n1000.m01 = wf(N=1000,mu=0.01,warmup=10000,timesteps=20000,output="progeny")
wf.prog.n3000.m01 = wf(N=3000,mu=0.01,warmup=10000,timesteps=20000,output="progeny")

# changing innovation rate
wf.prog.n1000.m005 = wf(N=300,mu=0.005,warmup=10000,timesteps=20000,output="progeny")
wf.prog.n3000.m001 = wf(N=300,mu=0.001,warmup=10000,timesteps=20000,output="progeny")

```


```{r, plot_progeny_wf, results='hide',fig.width= 9, fig.height= 5}
par(mfrow=c(1,2))
plot(wf.prog.n300.m01$d2,pch=20,log="xy",ylab="Probability P(k) of number variants >= k",xlab="k",col="darkgrey",type="b")
lines(wf.prog.n1000.m01$d2,pch=2,col="indianred",type="b")
lines(wf.prog.n3000.m01$d2,pch=3,col="royalblue",type="b")
legend("bottomleft",legend=c(expression(paste("N=300; ",mu,"=0.01")),
                             expression(paste("N=1000; ",mu,"=0.01")),
                             expression(paste("N=3000; ",mu,"=0.01")))
       ,pch=c(20,2,3),col=c("darkgrey","indianred","royalblue"),bty="n")


plot(wf.prog.n300.m01$d2,pch=20,log="xy",ylab="Probability P(k) of number variants >= k",xlab="k",col="grey",type="b")
lines(wf.prog.n1000.m005$d2,pch=2,col="indianred",type="b")
lines(wf.prog.n3000.m001$d2,pch=3,col="royalblue",type="b")
legend("bottomleft",legend=c(expression(paste("N=300; ",mu,"=0.01")),
                             expression(paste("N=300; ",mu,"=0.005")),
                             expression(paste("N=300; ",mu,"=0.001")))
       ,pch=c(20,2,3),col=c("darkgrey","indianred","royalblue"),bty="n")
```

### Progeny Distribution under object-mediated transmission

We explore whether and how the shape of the progeny distribution is affected by an unbiased object-mediated transmission with encoding and decoding errors. With a sufficiently large value of $T$ stochastic differences between runs become negligible, hence below we explore parameter combinations running each model one time over 10,000 generations.   


```{r parameters_progeny_exp1}
mu <- 0.01
N <- 300
warmup <- 10000
timesteps <- warmup + 10000
lambda <- c(0.1,1,5,10)
```

We first keep $N$ fixed to `r N[1]` and $\mu_{d}=\mu=$`r mu`, exploring different settings of $\lambda$:

```{r, execution_progeny_object_exp1}
decode.prog.varlambda <- vector("list",length=length(lambda))
encode.prog.varlambda <- vector("list",length=length(lambda))
legItems.progeny1=vector(length=length(lambda))

for (i in 1:length(lambda))
{
  decode.prog.varlambda[[i]] <-objTr(N=N,mu.e=0,mu.d=mu,
                                     lambda=lambda[i],warmup=10000,
                                     timesteps=20000,output="progeny")
  encode.prog.varlambda[[i]] <-objTr(N=N,mu.e=mu,mu.d=0,
                                     lambda=lambda[i],warmup=10000,
                                     timesteps=20000,output="progeny")  
  legItems.progeny1[i] = as.expression(bquote(paste(lambda,"=",.(lambda[i]),"; N=",.(N),sep="")))
}

legItems.progeny1=c(as.expression("Wright-Fisher"),legItems.progeny1)
```



```{r, parameters_progeny_exp2}
n.objects = 1000
N <- c(50,500,1000,2000)
lambda <- n.objects/N
```

We then look at instances where the number of objects produced is hold constant to approximately `r n.objects`, that is holding constant the product $\lambda \cdot N$:


```{r, execution_progeny_object_exp2}
encode.prog.fixobjects <- vector("list",length=length(N))
decode.prog.fixobjects <- vector("list",length=length(N))
legItems.progeny2=vector(length=length(N))


for (i in 1:length(N))
{
  decode.prog.fixobjects[[i]] <- objTr(N=N[i],mu.e=0,mu.d=mu,
                                       lambda=lambda[i],warmup=10000,
                                       timesteps=20000,output="progeny")
  encode.prog.fixobjects[[i]] <- objTr(N=N[i],mu.e=mu,mu.d=0,
                                       lambda=lambda[i],warmup=10000,
                                       timesteps=20000,output="progeny")  
  legItems.progeny2[i] = as.expression(bquote(paste(lambda,"=",.(lambda[i]),"; N=",.(N[i]),sep="")))

}
legItems.progeny2=c(as.expression("Wright-Fisher"),legItems.progeny2)
```

#### Results Decoding Error

```{r,plot_progeny_decode,fig.width=9,fig.height=5}
cc <- brewer.pal(length(lambda),'Set1')
par(mfrow=c(1,2))
plot(wf.prog.n300.m01$d2,type="n",lty=4,log="xy",ylab="Probability P(k) of number progeny >= k",xlab="k",col="darkgrey")
for (i in 1:length(decode.prog.varlambda))
{
  lines(decode.prog.varlambda[[i]]$d2.objects, type="b",pch=20,col=cc[i])
}
lines(wf.prog.n300.m01$d2,lty=4,lwd=2)
legend("bottomleft",legend=legItems.progeny1,col=c(1,cc),pch=c(NA,rep(20,length(lambda))),lty=c(4,rep(NA,length(lambda))),lwd=c(2,rep(NA,length(lambda))),bty="n")

plot(wf.prog.n1000.m01$d2,type="n",lty=4,log="xy",ylab="Probability P(k) of number progeny >= k",xlab="k",col="darkgrey")
for (i in 1:length(decode.prog.fixobjects))
{
  lines(decode.prog.fixobjects[[i]]$d2.objects, type="b",pch=20,col=cc[i])
}
lines(wf.prog.n1000.m01$d2,lty=4,lwd=2)
legend("bottomleft",legend=legItems.progeny2,col=c(1,cc),pch=c(NA,rep(20,length(N))),lty=c(4,rep(NA,length(N))),lwd=c(2,rep(NA,length(N))),bty="n")
```

The results on the left panel indicate that with low values of $\lambda$ the progeny distribution is very similar to the one expected from the Wright-Fisher model, with a power law section for lower values of $k$ followed by an exponential cut-off. When $\lambda\geq 5$, the exponential relationship starts at larger values of $k$, albeit with a slope similar to the one expected by the Wright-Fisher model. This is simply caused by the fact that when $\lambda>1$, obtaining $k=1$ instances of a variant becomes more difficult, as even if a mental representation is associated to a single individual for one generation, the number of objects with such mental representation is more likely to be larger than 1 (e.g. with $\lambda=5$ the probability of producing more than one copy of a variant is `r 1-ppois(1,lambda=5)`).

The results of right panel considers different permutations of $\lambda$ and $N$ producing approximately the same number of objects for each generation. The progeny distributions confirm the effect of large values of $\lambda$ --- in this particular case the smallest observed case of $k$ is just `r decode.prog.fixobjects[[1]][[1]][1,1]` which is observed only `r decode.prog.fixobjects[[1]][[1]][1,2]` `r if(decode.prog.fixobjects[[1]][[1]][1,2]==1){"time"}` `r if(decode.prog.fixobjects[[1]][[1]][1,2]>1){"times"}`. The other parameter settings have all comparatively low values of $\lambda$, and hence the resulting distribution have a close match to the expectation under Wright-Fisher. 

#### Results Encoding Error

```{r,plot_progeny_encode,fig.width=9,fig.height=5}
cc <- brewer.pal(length(lambda),'Set1')
par(mfrow=c(1,2))
plot(wf.prog.n300.m01$d2,type="n",lty=4,log="xy",ylab="Probability P(k) of number progeny >= k",xlab="k",col="darkgrey")
for (i in 1:length(encode.prog.varlambda))
{
  lines(encode.prog.varlambda[[i]]$d2.objects, type="b",pch=20,col=cc[i])
}
lines(wf.prog.n300.m01$d2,lty=4,lwd=2)
legend("bottomleft",legend=legItems.progeny1,col=c(1,cc),pch=c(NA,rep(20,length(lambda))),lty=c(4,rep(NA,length(lambda))),lwd=c(2,rep(NA,length(lambda))),bty="n")

plot(wf.prog.n1000.m01$d2,type="n",lty=4,log="xy",ylab="Probability P(k) of number progeny >= k",xlab="k",col="darkgrey")
for (i in 1:length(encode.prog.fixobjects))
{
  lines(encode.prog.fixobjects[[i]]$d2.objects, type="b",pch=20,col=cc[i])
}
lines(wf.prog.n1000.m01$d2,lty=4,lwd=2)
legend("bottomleft",legend=legItems.progeny2,col=c(1,cc),pch=c(NA,rep(20,length(N))),lty=c(4,rep(NA,length(N))),lwd=c(2,rep(NA,length(N))),bty="n")
```

As for the decoding error scenario we first explore the effect of different $\lambda$ values on the shape of the progeny distribution, then combinations of $\lambda$ and $N$ yielding the same number of objects per generation (i.e. holding $\lambda \cdot N$ constant). The progeny distribution of the objected mediated transmission again deviates from the Wright Fisher model with increasing values of $\lambda$ (left panel). While the power-law and exponential cut-offs sections are retained in all cases, the former is disrupted by two processes. Firstly even with high values of $\lambda$, mutations occurring at the moment of production allows for a large number of variants appearing only once ($k=1$). However, at $k=2$ we observe a sudden drop in the number of variants. This is because an object can appear more than once if it has also a mental representation in the population (i.e. there was a successful case of transmission), but with larger $\lambda$, lower values of $k$ are hard to achieve for the same reason as in the decoding error scenario. For example with $\lambda=5$, an object can appear in the assemblage twice ($k=2$) either because an individual produced two instances of its mental representation in a given generation (a probability of `r dpois(2,5)`) but possessed the mental representation only for that particular generation, or produced two instances in two generations (`r dpois(1,5)*dpois(1,5)`) before updating its mental representation. Either option is comparatively rare, and hence the number of variants appearing only few times in the assemblage is rare with larger values of $\lambda$. Indeed, the power-law section of the progeny distribution seems to approximately start from $k=\lambda$. 

The results of the experiment with a constant values for $\lambda \cdot N$ (right panel) shows again how low values of $\lambda$ yield a progeny distribution similar to Wright Fisher model, although with a higher drop in the proportion of variants after $k>1$. In the extreme case of $\lambda=20$ we do not observe any variants between $k=1$ and $k= `r encode.prog.fixobjects[[1]][[1]][2,1]`$, and after that the progeny distribution does follows a different trajectory compared to the Wright-Fisher model.

# References
1. Bentley, R.A., Hahn, M.W., Shennan, S.J., 2004. [Random drift and culture change](https://doi.org/10.1098/rspb.2004.2746). Procedings of The Royal Society B 271, 1443–1450.
2. O’Dwyer, J.P., Kandler, A., 2017. [Inferring processes of cultural transmission: the critical role of rare variants in distinguishing neutrality from novelty biases](https://doi.org/10.1098/rstb.2016.0426). Phil. Trans. R. Soc. B 372, 20160426.