Alexandre Hyafil

Nicolas Baumard

Evoked and Transmitted Culture models: Using bayesian methods to infer the evolution of cultural traits in history

A central question in behavioral and social sciences is understanding to what extent cultural traits are inherited from previous generations, transmitted from adjacent populations or produced in response to changes in socioeconomic and ecological conditions. As quantitative diachronic databases recording the evolution of cultural artifacts over many generations are becoming more common, there is a need for appropriate data-driven methods to approach this question. Here we present a new Bayesian method to infer the dynamics of cultural traits in a diachronic dataset. Our method called Evoked-Transmitted Cultural model (ETC) relies on fitting a latent-state model where a cultural trait is a latent variable which guides the production of the cultural artifacts observed in the database. The dynamics of this cultural trait may depend on the value of the cultural traits present in previous generations and in adjacent populations (transmitted culture) and/or on ecological factors (evoked culture). We show how ETC models can be fitted to quantitative diachronic or synchronic datasets, using the Expectation-Maximization algorithm, enabling estimating the relative contribution of vertical transmission, horizontal transmission and evoked component in shaping cultural traits. The method also allows to reconstruct the dynamics of cultural traits in different regions. We tested the performance of the method on synthetic data for two variants of the method (for binary or continuous traits). We found that both variants allow reliable estimates of parameters guiding cultural evolution, and that they outperform purely phylogenetic tools that ignore horizontal transmission and ecological factors. Overall, our method opens new possibilities to reconstruct how culture is shaped from quantitative data, with possible application in cultural history, cultural anthropology, archaeology, historical linguistics and behavioral ecology.

Introduction Evoked and transmitted culture

It is standard in behavioral sciences to distinguish between evoked culture and transmitted culture [START_REF] Gangestad | Toward an integrative understanding of evoked and transmitted culture: The importance of specialized psychological design[END_REF][START_REF] Nettle | Beyond nature versus culture: Cultural variation as an evolved characteristic[END_REF][START_REF] Tooby | The adapted mind: Evolutionary psychology and the generation of culture[END_REF][START_REF] Sng | The behavioral ecology of cultural psychological variation[END_REF]. Transmitted culture refers to behaviors and beliefs that are mostly due to cultural transmission, either vertically from previous generations (parents, teachers, elders) or horizontally from peers and neighbors. Language is a typical example of transmitted culture as the best predictor of an individual's language is the language of her parents. By contrast, evoked culture refers to behaviors and beliefs that are mostly due to the expression of an evolved behavioral program in response to some environmental cues. A typical example here is the number of children per woman: the March 4, 2021 1/28 environment in which this individual lives (in terms of income, social support, interpersonal violence) is a better indicator of an individual's reproductive preference than the inherited culture (religion, social norms). For example, in response to economic development, most middle-income countries (such as Iran or Korea) have experienced a rapid transition from 6 children per woman to 2 children per woman, regardless of their original religious tradition and socio-moral norms [START_REF] Roser | Fertility -Our World In Data[END_REF][START_REF] Rosling | Factfulness: Ten Reasons We're Wrong about the World-and why Things are Better Than You Think[END_REF].

Obviously, transmitted and evoked cultures are located on the same continuum. All evoked cultural behaviors are also partially influenced by cultural transmission. For instance, an evolved program may trigger an early onset of sexual debut, but the courtship involved in finding a mate may be heavily influenced by transmitted traditions (e.g. the 'date' tradition in the 20th c. US culture). In the same way, all transmitted behaviors are also partially influenced by evoked factors. For instance, phonology has recently been shown to be influenced by dietary and behavioral practices through changes in bite configuration [START_REF] Blasi | Human sound systems are shaped by post-Neolithic changes in bite configuration[END_REF]. Even the content of artistic works can be shaped by the environment. For instance, it has been shown that people living in more favorable environments tend to be more trustworthy and to invest more in cooperation. In line with this observation, individuals in European portraits tend to display more prosocial signals in more economically developed periods in history traditions [START_REF] Safra | Tracking historical changes in trustworthiness using machine learning analyses of facial cues in paintings[END_REF]. What brings all these behaviors together, what makes them considered cultural, is that they exist because of some sort of inheritance: either people are influenced by the information (language, beliefs, etc.) they inherited from the previous generation or they are influenced by the environment (resources, institutions, pathogens) they also inherited from the previous generations.

Computational tools for assessing models of cultural evolution

Disantangling evoked and transmitted factors is often difficult. Variation between any two cultures may result from separate traditions which lead to the transmission of two different transmitted cultures or from different local environments which trigger different behavioral programs leading to different evoked cultures. Conversely, uniformity between two cultures could result from similar traditions coming from the same source or from the triggering of similar evolved mechanisms in similar environments.

Typical quantitative cultural evolution studies use a phylogenetic method to reconstruct when a particular trait emerged in the past [START_REF] Blasi | Human sound systems are shaped by post-Neolithic changes in bite configuration[END_REF][START_REF] Jordan | Matrilocal residence is ancestral in Austronesian societies[END_REF][START_REF] Mace | A phylogenetic approach to cultural evolution[END_REF][START_REF] Watts | Ritual human sacrifice promoted and sustained the evolution of stratified societies[END_REF]. Phylogenetic methods are however limited by the differences between cultural and genetic evolution. These methods do not take into account horizontal tranmission of cultural traits, i.e. the geographical diffusion of culture by communication between populations. Moreover, genetic evolution does not include evoked component, i.e. the role of ecological factors such socio-economic development, latitude, or the size of the population. A possibility to take all these factors into account is to use instead a regression method [START_REF] Blasi | Human sound systems are shaped by post-Neolithic changes in bite configuration[END_REF][START_REF] Botero | The ecology of religious beliefs[END_REF] that incorporates the influence of ecology, shared cultural background and spatial proximity onto the observed cultural trait. However these methods may not be as accurate as methods that rely on an explicit model of cultural evolution through time, and are necessarily arbitrary in what they define as shared cultural background. Importantly, they do not permit to reconstruct the evolution of a cultural trait through time and across regions.

Cultural evolution with diachronic datasets

Another important characteristic of these studies is that cultural traits as well as ecological factors are observed only once (synchronic dataset), usually at the time of the study. However, the development of long-term series of socio-economic indicators [START_REF] Bosker | From baghdad to London: Unraveling urban development in Europe, the middle East, and North Africa[END_REF][START_REF] Buringh | Charting the "rise of the west": Manuscripts and printed books in Europe, a long-term perspective from the sixth through eighteenth centuries[END_REF][START_REF] Fouquet | Seven Centuries of European Economic Growth and Decline[END_REF][START_REF] Morris | The measure of civilization: How social development decides the fate of nations[END_REF], combined with quantitative record of cultural traits and cultural March 4, 2021 2/28 productions (diachronic dataset) [START_REF] Baumard | Increased Affluence Explains the Emergence of Ascetic Wisdoms and Moralizing Religions[END_REF][START_REF] Turchin | Seshat: The global history databank[END_REF], now makes it possible to directly assess the impact of ecological and transmitted factors in cultural evolution. To address these limitations, we describe a general statistical framework to infer the specific causes of cultural evolution in a diachronic or synchronic dataset. The framework relies on a generative model for the evolution of cultural traits that incorporates evoked component and spatial diffusion.

Materials and methods

Defining ETC models

We first distinguish between cultural traits and cultural artifacts. Cultural artifacts correspond to the material production of a given population: tools, songs, portraits, watches, novels, etc. (see Figure 1). By contrast, cultural traits are the cognitive and behavioral dispositions shared by the member of the population: social trust, belief in a moralizing god, conservativeness, etc.. While these traits are not directly observable to the social scientist, they strongly influence how cultural artifacts are produced. Cultural traits can themselves be influenced by several factors: they can be transmitted from generation to generation (vertical cultural transmission), transmitted from one society to its neighbor (horizontal cultural transmission) or evoked by ecological factors (socio-economic development, political organization, etc.).

A general question in cultural evolution is to understand the relative importance of each factor using the indirect observations of cultural artifacts. We formalize this problem under the generic name of Evoked and Transmitted Culture (ETC) models, using the general formalism of probabilistic generative models. ETC models are special forms of latent process models (also called state-space or Hidden-Markov models), where an unobserved (latent) process evolves through time and determines the value of some observed variables [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF]. Here in general, we take a set or K interconnected geographical regions and want to understand cultural evolution in each of these regions within a certain period of time. For practical reasons we use here a set of n discrete time points with time step dt. Cultural trait takes a certain value T t,r each for region r and each time point t, and its evolution is guided by its value at the preceding time step, as well as the value in neighbouring regions N(r) and the value of ecological factors E t,r :

p(T t,r |T t-1 , E t,r ) = f ((T t,r , {(T t,s , G r,s ), s ∈ N(r)}, E t,r ; θ T )) (1)
G r,s represents the connectivity between regions r and s in the set of regions N(r),

i.e. G r,s is the level of cultural communication between the two regions. θ T represents a set of parameters of cultural evolution. The overall communication between regions is summarized in the cultural connectivity matrix G, which here is taken to be constant in time. It is straightforward to extend the framework to use instead time-varying connectivity G(t), for example to account for the opening and closing of channels of communication between regions.One important implication of Equation 1 is that cultural trait is a Markov process, i.e. it depends on the past only through its value at the previous time step.

Defining the temporal limits and the borders of a specific 'culture' is tricky. What count as a culture? What counts as a new culture? For pragmatic reasons, we will consider two cases for creation of a new culture. Either it is created de novo (from populations for which we do not have previous records in the database), or it is the offspring of another region (i.e. when a population moves en masse to the new area or split into different populations whose cultural dynamic will now diverge). In case of a region created de novo, we will assume that the value of the cultural trait for the new region at its time of creation t 0 is sampled from some initial distribution ). In the case of an offspring region, we will assume that the value of the new region r is directly transferred from its parent region s, i.e. T r,t = T s,t .

The production of artifacts is guided by the value of the cultural traits, and possibly as well by some other (known) factors F r,t that can relate to the material conditions of artifact productions in region r at time t (e.g. the number of universities in the region if artifacts are the number of scientific publications in a certain field) and to the data collection process (e.g. the number of registered archaeological sites when assessing the fabrication of a particular tool in a region)

p(A t,r ) = g(T t,r , F t,r ; θ A ) (2)
where θ A are the parameters of the artifact production process.

While in general we want to infer the trajectory of unobserved cultural traits from the observations of cultural artifacts, it is also possible that the values of cultural traits are accessible at some moment (typically at contemporary times, i.e. at the final time step), and we can include such information to further constraint the ETC model.

In the following sections we discuss two distinct classes of ETC models: when the cultural traits take the form of a binary value (T t,r ∈ {0, 1}), and when they take the form of a continuous value (T t,r ∈ R).

Binary cultural traits

In many circumstances it is convenient to model a cultural trait as a binary value (T t,r ∈ {0, 1}), such as the presence or absence of the trait in the population [START_REF] Baumard | Increased Affluence Explains the Emergence of Ascetic Wisdoms and Moralizing Religions[END_REF]. In this case the probability of trait transitions between the two consecutive time steps t and t + dt is defined by time-dependent rates of transitions R up (t, r) and R down (t, r):

p(T t+1,r = 1|T t,r = 0) = R up (t, r)dt p(T t+1,r = 0|T t,r = 0) = 1 -R up (t, r)dt p(T t+1,r = 0|T t,r = 1) = R down (t, r)dt p(T t+1,r = 1|T t,r = 1) = 1 -R down (t, r)dt (3)
Here the approximation is valid only if both rates are much smaller at each time than the inverse of the time step (R up (t, r), R down (t, r) 1/dt ∀(t, r)).

Initial values are parametrized by p(T 1,r = 1) = x 0 (and also for regions created de novo). The rates of transitions depend on the environmentally shaped susceptibility of innovation (transition from trait absent to trait present) and loss of trait (transition from trait present to trait absent), as well as diffusion from neighbouring areas. We propose here linear effects:

R up (t, r) = b up + λE t,r + ξ( s∈N(r)|Tt,s=1 G r,s - s∈N(r)|Tt,s=0 G r,s ) R down (t, r) = b down -λE t,r + ξ( s∈N(r) |Tt,s=0 G r,s - s∈N(r)|Tt,s=1 G r,s ) (4)
where b up is the intrinsic rate of cultural innovation, and b down is the intrinsic rate for loss of the cultural trait, λ is the sensitivity of the cultural transition to the ecological factor and ξ determines the rate of cultural diffusion. These five parameters θ T = (x 0 , b up , b down , λ, ξ) determine the dynamics of the cultural traits. Equation 4 can be rewritten in terms of the multivariate rate of transitions for all regions R up (t) and R down (t) :

R up (t)dt = b up 1 + λE t + ξG(2T t-1 -1) R down (t)dt = b down 1 -λE t -ξG(2T t-1 -1) (5)
March 4, 2021 5/28

Combining equations 3 and 5 allows to compute the time-dependent transition rate p(T t+1 |T t ; θ T ). This system is equivalent to an input-output HMM [START_REF] Bengio | An Input Output HMM Architecture[END_REF] where the latent variable is defined as the combination of the cultural trait for all regions, thus taking 2 K possible values: T t ∈ {0, 1} K .

Continuous cultural traits

Here we present the case when trait takes a continuous value instead of binary value.

Again, we can define the dynamics of the cultural trait as following the linear influence of the ecological factor and value of the trait in neighbouring regions, as captured in the following differential equation (in continuous time)

dT t,r dt = -ρT t,r + γ + λE t,r + ξ s∈N(r) G r,s (T ts -T t,r ) + η (6)
ρ is the cultural lability (or cultural leak, i.e. the inverse of the intrinsic time scale of cultural evolution), γ is a bias towards positive or negative trait values and η is some gaussian white noise capturing the stochastic component in cultural evolution. As for the binary case, λ represents susceptibility to the external factor, and ξ the rate of cultural diffusion. Equation 6 is transformed in discrete time as a linear-Gaussian evolution equation, taking the form:

T t,r = (1 -ρdt)T t-1,r + dt γ + λE t,r + ξ s∈N(r) G r,s (T t-1,s -T t-1,r ) + η √ dt (7)
Cultural trait is initiated independently at each region from a gaussian distribution centered on x 0 and of variance σ 2 0 . The same applies to regions created de novo, while regions created from a parent simply inherit the value of the trait for the parent region.

Equation 7 can be reframed in terms of the K-dimensional vector T t :

T t = (1 -ρdt)T t-1 + dt(γ1 + λE t + ξ GT t-1 )) + η √ dt, η ∼ N (0; I) ( 8
)
where G is defined by Grr = G rrs G rs and Grs = G rs for r = s.

Thus T is a multivariate auto-regressive model of order one, its dynamics is given by:

   T t = M T t-1 + J t + η √ dt, with M = (1 -ρdt)I + ξdt G and J t = (λE t + γ1)dt (9) Initial state is T 0 ∼ N (x 0 1, σ 2 0 I).
The set of parameters for the cultural trait evolution is θ T = (γ, λ, ξ, ρ, σ 2 , x 0 , σ 2 0 ).

The first five parameters relate to the cultural trait evolution (θ T e ), while the last two refer to the initialization (θ T i ). Depending on the particular ETC problem, there may also be parameters associated with the artifact production model, that we denote generically θ A .

General fitting procedure

We are generally interested in estimating from a dataset (including cultural artifacts A, ecological factors E as well as the cultural connectivity G) the set of parameters θ = (θ T , θ A ) from the ETC model. Here we define A and T as the entire collection of cultural artifacts and traits, i.e. A = (A 1 ..A n ) and T = (T 1 ..T n ). Best fitting parameters report the degree of ecological determination, vertical and horizontal transmission of the cultural traits, and allow estimating the dynamics of the cultural March 4, 2021 6/28 traits across time. They also allow comparing different models, for example to compare how each ecological factor can account for the variation of cultural traits.

There is in general no analytical form to find a solution to this problem, except for some cases such as HMM (see Appendix A) or Kalman filters (i.e. continuous trait variable with continuous artifacts generated as gaussian observations). For other types of model, we will have to recur to approximate solutions, using a variety of techniques such as analytical approximations, variational approaches or sampling methods [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF].

Most of these methods rely on variations of the Expectation-Maximization (EM) algorithm (although see [START_REF] Paninski | A new look at state-space models for neural data[END_REF] and [START_REF] Buesing | Spectral learning of linear dynamics from generalised-linear observations with application to neural population data[END_REF] for alternative approaches). The EM algorithm runs recursively between two stages. In the Expectation stage, we compute the posterior probability for the value of the (unobserved) cultural traits, given the (observed) cultural artifacts and the current value of the parameters, i.e. q(T ; θ) = p(T |A, θ). In case the value of cultural traits are observed at some points, inference will also be based on these observed traits T obs : q(T t ; θ) = p(T t |A, T obs , θ). While exact inference can be performed when the cultural trait is binary (see Appendix A), it cannot in general when the cultural trait is continuous, so we have to resort to approximate methods to evaluate the posterior q(T t ; θ). We have investigated three types of approximations: forward-backward method, Laplace approximation and Expectation-Propagation (EP) (see Appendix B for details). In the Maximization stage, we look for the parameters θ that best account for the observed data, given the estimation of the value of cultural artifacts from the E-step, i.e we maximize Q(θ; θ old ) = q(T ; θ old ) log p(A, T |θ)dT w.r.t θ:

Q(θ; θ old ) = q(T ; θ old ) log p(A, T |θ)dT = q(T ; θ old ) log p(T |θ T )dT + q(T ; θ old ) log p(A|T , θ A )dT (10)
We define Q T (θ T ) and Q A (θ A ) as resp. the first and second integrals in the right hand side. We can see that trait parameters θ T and production parameters θ A can be estimated independently by maximizing Q T and Q A separately. Since trait evolution is a Markov process, p(T ) = p(T 0 )Π t p(T t |T t-1 ), Q T can further be dissociated into initialization and evolution terms:

Q T (θ T ; θ old ) = q(T 0 ) log p(T 0 |θ T i )dT + n t=1 q(T t , T t-1 ) log p(T |T t-1 , θ T e )dT = Q T i (θ; θ old ) + Q T e (θ T e ; θ old ) (11) Q A (θ A ) can
also be decomposed into a sum of terms depending on each artifact independently:

Q A (θ A ; θ old ) = q(T ) log p(A|T , θ A )dT = t,r q(T t,r ) log p(A t,r |T t,r , θ A )dT t,r (12)
The algorithm stops when the change in log-likelihood log L(θ) = log p(A|T , θ) of the parameters between two successive iterations is below a certain threshold. Since the (log-) likelihood function is in general non-convex, the EM procedure can converge to a local (non global) minimum. To mitigate this problem, multiple EM loops are launched, each with a different set of initival values for parameters θ. The set of final parameters with overall higher likelihood is retained. Once we have identified best-fitting parameters θ, we can infer the trajectory of cultural traits across time in the different region, which is given directly by the result of the E-step q(T ; θ). We also define measures of influence, which compare how much variance in the cultural trait value can be explained by the value of the cultural trait in the same region one generation before (influence of vertical transmission), by the value of the cultural trait in neighbouring regions (influence of horizontal transmission) and the value of the ecological factor in the same region (influence of environment, see Appendix B.5 for details).

March 4, 2021 7/28

Confidence interval and model comparison

Once we have obtained the best-fitting values θ for the ETC parameters, we are generally interested in deriving some confidence intervals for these values. We will use two different ways to compute them. Since in a latent model observations are not independently drawn, classical bootstrapping methods cannot be used; we can use parametric bootstrapping instead [START_REF] Visser | Confidence intervals for hidden Markov model parameters[END_REF]. In parametric bootstrapping, we evaluate the bias in parameter estimation that we should expect if the true model corresponded indeed to the ETC model with best-fitted values. We will generate n B different samples from the ETC parametrized with best fitting value, i.e. we sample from A| θ. Each sample has the exact same number of observations at each data point as in the real dataset. Then the fitting procedure is applied to each of the synthetic datasets. The distribution θ bootstraps of the parameters reflects the intrinsic uncertainty of the estimation process for our dataset. It can be used directly to extract confidence interval for each parameter.

An alternative strategy is to use the Laplace approximation to approximate the posterior distribution over parameters p(θ|A) by a Gaussian distribution centered on the MAP parameters N ( θ, Σ) [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF]. Posterior distribution for each parameter is the approximated by N (θ i , Σ ii ), and significance for the parameter vs. null hypothesis θ i = 0 can be tested using the Wald test. Σ is taken to be the inverse of the Hessian of -log L(θ) taken at the MAP parameters, i.e. it is defined by the local curvature of the log-likelihood function.

Following [START_REF] Jamshidian | Standard errors for EM estimation[END_REF], the Hessian H of the log-likelihood function evaluated through an EM algorithm can be computed following

H = ∇∇ log L(θ) = ∇∇Q( θ) (I -∇ arg max θ Q(θ; θ old )| θ old = θ (13)
The Jacobian can be evaluated by methods of finite differences, while the Hessian of Q, because of the decomposition of Q in equation 11, decomposes also into Hessians related to Q T and Q A :

∇∇Q( θ) = ∇∇Q T ( θT ) 0 0 ∇∇Q A ( θA ) (14)
Both Hessians can be computed analytically in a variety of cases (see sections A.2 and B.2).

Model comparison permits to compare the full version of the model to reduced versions, where either diffusion between regions and/or evoked component is removed.

To perform model comparison, we can use the log model evidence (marginalized over parameter vlaues), which can be approximated using the Laplace approximation (and assuming a broad prior for parameters θ) by:

log p(A) ≈ log p(A| θ) + 1/2 log |H| (15)
Other methods that penalize the log-likelihood based on the numbers of parameters such as the Akaike Information Criterion (AIC) or Bayesian Information Criterion (BIC). Alternatively, model comparison can be performed using cross-validation, which is more expensive computationally but will be more accurate when the Laplace approximation is not appropriate, notably when there are multiple local maxima in the LLH. Here we used k-fold cross-validation, where model performance was assessed using Cross-Validated Log-Likelihood.

March 4, 2021

Results

Testing with binary ETC model

We tested our estimation method on synthetic data generated from an ETC with binary traits and binary artifacts. The model was run over 6 interconnected regions over 1000 time steps. We compared the true value with the estimated values for cultural traits and parameters of the model (Figure 2a). The fitting algorithm captured the overall dynamics of the up and down states of the cultural traits in all regions. It also provided reliable estimates of the parameters of the ETC model (Figure 2b). Finally we performed Bayesian model comparison of the full ETC model with reduced model where we either removed horizontal transmission (evoked-only model, where we enforced ξ = 0) or evoked component (diffusion-only model, where we enforced λ = 0). BIC was larger in either reduced model than in the full model (∆BIC ≈ 15), showing that the Bayesian procedure correctly identified that both horizontal transmission and evoked component contributed to the evolution of the cultural trait (Figure 2c).

Testing continuous ETC model

Subsequently, we tested the estimation method for ETC models with a continuous trait, and artifacts with count values (e.g. the number of paintings of a specific kind produced in a given region and year). We simulated R = 5 regions over 2000 years with time step of dt = 50 years. We fitted ETC models from a batch of 100 simulated datasets using the same parameter set, and for each dataset we compared three different methods for approximating the posterior latent distribution during the E-step of the EM algorithm: forward-backward method, Laplace approximation and Expectation-Propagation (EP) (see Appendix B).

We found that model fitting based all approximations provided very good estimates of true parameters and of the trajectory of cultural traits (Figures 3 and4f). All three approximation methods also provided good estimates of true parameters, with little to no bias (Figure 4a estimates that the forward-backward method, with less bias and less variance in parameter estimates. The Laplace and EP approximations provided indeed very similar results, with very high correlation between parameter estimates between the two methods: r > 0.99 for all 7 parameters expect noise parameter σ (r = 0.94). By contrast, correlation estimates were much lower between estimates given by the forward-backward method and either the EP or Laplace approximation (r in the range 0.2-0.8, except for σ 2 0 and x 0 where r values were in the range 0.94-0.97). We also tested the capacity of the method to detect from the data whether diffusion of the cultural trait between regions and/or evoked component are present. For such hypothesis testing, we compared which of four alternative models was better supported by the synthetic data: the full ETC model, that included diffusion and evoked component; the no diffusion model, that included evoked component but no diffusion (i.e. parameter ξ is set to 0); the no evoked model, that included diffusion but no evoked component (i.e. parameter λ set to 0); and the control model, where diffusion and evoked component are both absent. We compared different metrics to perform model selection: metrics based on penalized log-likelihood (Akaike Information Criterion, AIC; Bayesian Information Criterion, BIC), on cross-validation (Cross-Validated Log-Likelihood, CVLL), or on permutation tests on individual parameters on the full model (boostraps). We report in Figure 5 the proportion of datasets where each of the model is selected, for each metric and each approximation algorithm. Results were similar using either Laplace or EP approximation. For both approximations, the full model was correctly identified in a majority of cases for all metrics. Cross-validation (CVLL) provided the most sensitive measure, as it correctly identified the full model in almost all datasets. Using the forward-backward estimation instead, the full model was correctly identified in a smaller proportion of datasets irrespective of the metrics used. This confirms that the EP and Laplace approximations outperformed the forward-backward method for parameter estimation and model selection. ETC generates better estimates than two-stage analysis

Next, we compared our approach to a simpler two-stage approach where cultural traits are first estimated using cultural artifacts produced in the same region and at the same time, and parameters are later estimated directly based on these estimated latents. This second stage resembles the fitting procedure used in [START_REF] Whitehouse | Complex societies precede moralizing gods throughout world history[END_REF], applied directed to cultural traits rated by different experts for each region and time. The estimated cultural traits using this method differ from the EM algorithm in two important aspects: they are only constrained by cultural artifacts (whereas in the EM algorithm they are constrained by both artifacts and estimated parameters); they do no include any level of uncertainty.

Based on both of these aspects, we can expect poorer estimates: ignoring reliability of the estimates leads to weighing equally data points with many cultural artifacts and data points with a scarcity of artifacts. We compared the two methods on synthetic data generated from an ETC model with 5 regions over 500 time steps, producing continuous-valued artifacts. We separated the time frame into a longer period of scarcity of artifacts (first 400 times steps; mean artifact per time point: 0.1), and a shorter period of abundance of artifacts (last 100 time steps; mean number of artifacts: 1).

Cultural traits estimated from both methods are compared in Figure 6a. We see that the one-stage method provides good estimates of the true cultural traits. It also displays more uncertainty for data points with scarce cultural production (left to the black vertical bar) compared to data points with an abundance of cultural artifacts (right to the black vertical bar). By contrast, the two-stage method provides less smooth and less accurate estimates of cultural traits. We also compare the estimated parameters from both methods in Figure 6b for 30 different simulations. The one-stage method provides very reliable and unbiased estimates, especially for the parameters defining the dynamics of cultural trait evolution (γ, λ, η, ρ, σ 2 ). Moreover, as predicted, the two-stage method provides very biased and unreliable estimates. We also found that the mean-squared error term over the cultural trait was on average much lower in the one-stage than in the two-stage procedure (Figure 6c). We attribute this to the fact that the two-stage method estimates the dynamic regime of the cultural trait without taking into account the uncertainty about cultural trait.

Cultural evolution can also be inferred from synchronous data

Finally, we studied whether ETC model parameters could also be reliably inferred from synchronous instead of diachronous data. We simulated the evolution of a cultural trait in regions forming a phylogenetic tree. The evolution of the cultural trait was influenced by the value of a binary ecological factor (Figures 7a and7b). We then fitted ETC parameters to the cultural traits observed in all regions at the final time of the simulation. We compared two versions of the model: either the full ETC model, or the model without spatial diffusion, which is very similar to the classical phylogenetic method previously used on cultural datasets [START_REF] Jordan | Matrilocal residence is ancestral in Austronesian societies[END_REF][START_REF] Mace | A phylogenetic approach to cultural evolution[END_REF][START_REF] Watts | Ritual human sacrifice promoted and sustained the evolution of stratified societies[END_REF]. Despite using only the final points of the cultural evolution, the method reliably reconstructed the trajectory of the cultural traits along the tree (Figure 7c). The reconstruction error increased when cultural diffusion was not included into the dynamics of cultural evolution (Figures 7d and7f). Parameters were also estimated relatively well using the full ETC model (Figure 7e), with little to no bias. By contrast, ignoring cultural diffusion lead to systematic biases in parameter estimation, especially for the cultural lability parameter ρ, which controls vertical transmission. We interpret this as an incorrect assignment of cultural transmission: because horizontal transmission is not taken into account, similarities between the cultural trait in different regions are uniquely attributed to vertical transmission from a common parent.

March 4, 2021 12/28

Discussion

There has been a recent interest in using quantitative approaches, and in particular model fitting, to test hypotheses about the evolution of culture and cultural traits at a microscopic [START_REF] Beppu | Iterated Learning and the Cultural Ratchet[END_REF][START_REF] Kandler | Inferring individual-level processes from population-level patterns in cultural evolution[END_REF] or macroscopic scale [START_REF] Jordan | Matrilocal residence is ancestral in Austronesian societies[END_REF][START_REF] Mace | A phylogenetic approach to cultural evolution[END_REF][START_REF] Botero | The ecology of religious beliefs[END_REF]. While phylogenetic approaches have been applied in certain cases [START_REF] Jordan | Matrilocal residence is ancestral in Austronesian societies[END_REF][START_REF] Mace | A phylogenetic approach to cultural evolution[END_REF][START_REF] Watts | Ritual human sacrifice promoted and sustained the evolution of stratified societies[END_REF][START_REF] Botero | The ecology of religious beliefs[END_REF], there are two major differences between the dynamics of cultural and genetic evolution. First, unlike genetic evolution, horizontal transmission is an important factor of cultural evolution, competing with vertical transmission (i.e. transmission from previous generation). Second, unlike Darwinian evolution, cultural evolution may be influenced directly by the environment via phenotypic plasticity [START_REF] Baumard | Psychological origins of the industrial revolution[END_REF]. Thus, quantitative models of cultural evolution should capture these three types of driving forces: vertical transmission, horizontal transmission and influence of the environment. We have presented here a new quantitative model termed ETC which captures how these three driving forces can shape the evolution of a cultural trait. In ETC, the importance of these factors are captured by specific parameters. Importantly, we have shown that our model fitting procedure allows to estimate reliably model parameters from a diachronic database, where the evolution of cultural traits are not directly observed but inferred from the production of cultural artificats in different regions and different times. The ETC framework allows to capture different types of cultural traits (as binary or continuous variable) and cultural artifacts (binary, counts, continuous variable, ...). At the current time it models the cultural trait as a single metric over the entire population in one region at each time step, but possible extensions could include richer representations of the cultural traits such as stratified cultural traits by generation [START_REF] Fogarty | The life history of learning: Demographic structure changes cultural outcomes[END_REF] or the encoding of various cultural traits [START_REF] Watts | Ritual human sacrifice promoted and sustained the evolution of stratified societies[END_REF].

The fitting procedure relies on an Expectation-Maximization algorithm, where in the E-step two approximate inference algorithms have shown to yield the best results (for continuous traits): the Laplace approximation and Expectation-Propagation. These approximations work better than the forward-backwardmethod [START_REF] Brown | A statistical paradigm for neural spike train decoding applied to position prediction from ensemble firing patterns of rat hippocampal place cells[END_REF] probably because they are global approximations, i.e. they approximate directly the distribution of the collection of trait values at all time steps T . By contrast, the forward-backward method relies on local approximations, i.e. approximate the trait value at each time step T t sequentially, allowing for errors to accumulate from time step to time step. The forward-backward method has a lower computational cost however, and especially it scales better with the number of time steps. It thus provides a good alternative for the other methods for large datasets.

We also found that the fitting procedure with EM algorithm is much more reliable than a two-stage procedure where cultural traits are first estimated from cultural artifacts, and ETC parameters are then fitted on these cultural traits (Figure 6). The major reason is that this second method ignores altogether the uncertainty about the cultural trait, which strongly biases the estimates of the dynamic rules that dictate the evolution of these traits. This suggests that fitting and comparing models of cultural evolution based on datasets of cultural traits that ignores the necessary error-in-measure (whether such error comes from limited samples of cultural artifacts or any imprecision in manual ratings of cultural traits) can lead to incorrect conclusions [START_REF] Whitehouse | Complex societies precede moralizing gods throughout world history[END_REF]. As a side note, the direct procedure has also three important benefits. First, the uncertainty about the cultural traits naturally reflects the abundance of cultural artifacts in the region (or neighbouring regions) at the time. Second, the temporal and spatial scale over which cultural artifacts are integrated to produce latent trait estimates are not defined a priori but directly adjusted to the ETC parameters (as implicitly defined from the E step of the EM algorithm). Finally, the dynamics of cultural traits is jointly constrained by cultural artifacts, and the ETC estimated parameters through dynamic rules of the ETC model. This constrasts with the two-step procedure where cultural traits are only inferred from cultural artifacts. March 4, 2021 15/28

While it is naturally designed to assess cultural evolution from diachronic datasets, the method also works to reconstruct the evolution of a cultural trait from a synchronous dataset, i.e. from observing this trait in different regions at a single time (Figure 7). In this sense the method bears many similarities to existing phylogenetic methods that have been previously employed for cultural datasets. However we found that it is better adapted than phylogenetic methods, as it includes the specificities of cultural evolution with respect to genetic evolution: notably ETC models can incorporate the influence of evoked component and cultural diffusion, unlike standard phylogenetic methods. We found that ignoring horizontal transmission (i.e. cultural diffusion) leads to misattributing shared common traits to vertical transmission (Figure 7). This induces a biased estimation of the cultural lability parameter that controls vertical transmission, as well as unreliable reconstructed evolution of the cultural trait. This suggests a cautious interpretation of previous studies that have used the phylogenetic method to infer co-evolution of cultural traits. For instance, while it is true that Austronesian societies have remained relatively isolated from each other during the last millennium [START_REF] Watts | Ritual human sacrifice promoted and sustained the evolution of stratified societies[END_REF], cultural transmission probably occurred many times before the arrival of the Europeans [31]. We believe that taking this horizontal transmission into account could significantly change the way we think cultural traits have evolved in these regions.

The natural application of the method is to study how cultures evolve in different populations or regions, but ETC models could be equally well deployed for quantitative analyses in many other fields. For example, in archaelogy, it could also be used for reconstructing the evolution of cultural traits in a network of connected settlements, where data is accessible to us through fragmented records with usually limited coverage in time and space [START_REF] Prignano | Wiring the Past: A Network Science Perspective on the Challenge of Archeological Similarity Networks[END_REF]. Comparison between ETC models with different underlying diffusion networks could be used in this context to infer the pathways of cultural diffusion. In behavioral ecology, ETC models could be applied to understand how cultural traits can diffuse between social groups and whether they can be shaped by ecological factors [START_REF] Connor | Two levels of alliance formation among male bottlenose dolphins (Tursiops sp.)[END_REF][START_REF] Mundinger | Animal cultures and a general theory of cultural evolution[END_REF][START_REF] Jones | Cultural Transmission[END_REF]. In studies of human cultures, the different nodes of the network do not need to correspond to geographical regions but could correspond to any kind of social group. For example, ETC models could be used to infer how and why specific cultural traits (e.g. linguistic traits) spread between different social classes, based on a language corpus. The nodes of the network could even be interconnected individuals. For example, ETC models could be used to understand how socioeconomic conditions shape the diffusion and the evolution of artistic practices within a group of artist [START_REF] Fraiberger | Quantifying reputation and success in art[END_REF].

Conclusion

We have proposed a class of model, called Evoked and Transmitted Culture (ETC) models, to capture the dynamics of evolution of cultural traits in a set of interconnected networks. Bayesian methods were shown to provide accurate tool to reconstruct such evolution from a (possibly fragmented) diachronic dataset of cultural artifacts, and also perform well on synchronic datasets. These methods will be available as a Matlab toolbox. We hope this new quantitative tool will help unveil the determinants of cultural evolution among societies, social groups and individuals.

March 4, 2021 16/28

It should be noted that the inference problem is only tractable when the number of regions K remains low. Indeed the transition matrix has size 2 2K , which is already around 10 6 for K = 10 regions. If dt is sufficiently small, we can make the assumption that at most one transition can occur at each time step. Thus each multivariate state can only map onto K + 1 states (the state itself plus all the states corresponding to transition in one region), which means that the transition matrix has no more than 2 K (K + 1) non-zero values. We note C t the transition from t to t + 1, with C t = r if the value of the trait in region r (and only region r) has changed between t -1 and t, and C t = 0 if no region has changed its state. The transition matrix then rewrites to:

p(C t |T t ) =    R up (t, r)dt if C t > 0 and T t,r = 0 R down (t, r)dt if C t > 0 and T t,r = 1 1 - K r=1 p(C t =r|T t ) if C t = 0 (17)
Then we note ζ t (ε, r) the posterior probability that T t = ε and that cultural trait in r has changed between t and t + dt for r ≥ 1, and ζ t (ε, 0) the posterior probability that T t+1 = T t = ε (i.e. no trait change between t and t + dt).

For large number of regions (more than 10-20), we must resort to approximate solutions using variational inference or sampling methods [START_REF] Ghahramani | Factorial hidden Markov models[END_REF].

A.2 M-step for binary traits

In the M-step we maximize Q T i (θ; θ old ), Q T e (θ T e ; θ old ) and Q A (θ A ; θ old ) with respect to new parameters θ:

             Q T i (θ; θ old ) = K r=1 (γ r,1 (0) log(1 -x 0 ) + γ r,1 (1) log(x 0 )) Q T e (θ T e ; θ old ) = N -1 t=1 ε∈{0,1} K ς∈{0,1} K ζ t (ε, ς) log p(T t+1 = ς|T t = ε, θ T ) Q A (θ A ; θ old ) = T p(T |θ old , A) log p(A|T , θ A ) = N t=1 T ∈{0,1} γ t,r (T ) log p(A t,r |T t,r = a, θ A ) ( 18
)
March 4, 2021 17/28

Maximizing Q T i (θ; θ old ) yields x 0 = r γ r,1 (1)/R.

Maximizing Q T e (θ T e ; θ old ) with respect to trait evolution parameters (b up , b down , λ, ξ) can be made use gradient methods. When dt is small enough that multiple trait changes in one time step are negligible, Q T e (θ T e ; θ old ) rewrites as:

Q T e (θ T e ; θ old ) = N -1 t=1 ε∈{0,1} K K Ct=0 ζ t (ε, C t ) log p(C t |T t = ε, θ T ) (19)
The gradient and Hessian (necessary for the Laplace approximation) are obtained analytically by extracting first-and second-order derivatives from equations 19, 17 and 5.

Maximizing Q A in equation 18 will depend on the particular form of the artifact production model. If each artifact can take a finite set of values A t,r ∈ {1, ..N A }, then we can define parameters p A (A, T ) as the probability that a given artifact takes value A if cultural trait is in state T . By noting a(t, r, A) the counts of cultural artifacts of value A observed in region r at time step t, Q A can be rewritten as:

Q A (θ A ; θ old ) = N A A=1 T ∈{0,1} a(t, r, A) log p A (A, T ) r,t γ t,r (T ) (20)
Maximizing Q A with the constraints A p A (A, T ) = 1 yields:

p A (A, T ) = t,r a(t, r, A)γ t,r (T ) t,r,B a(t, r, B)γ t,r (T ) (21)
A more general case will be that trait T t,r parametrizes as an exponential family distribution for artifact T t,r . For example it will set the probability of observing one value if T t,r is binary (or can take a finite set of values), or it will determine the expected value of T t,r if T t,r can take continuous values or represents counts. Such framework allows to enter other covariates F t,r into the production model. In such cases maximizing Q A to find the new values for θ A is equivalent to a GLM with observations {A t,r } and observations weights γ t,r .

Appendix B Fitting continous ETC model B.1 E-step for continuous traits

If traits are continuous and the production rule is itself linear-Gaussian, then the classical Kalman filter tools can be used to infer parameters and the trajectory of cultural traits through time and space. Otherwise the inference problem is intractable, and a variety of approximation methods exist. Here we will restrict to Gaussian approximations since they allow analytical treatment, i.e. we will look for posterior distribution for cultural traits in the form:

p(T t |A (1,.,n) , θ) ≈ q(T ) = N (T ; m, V ).
In particular we will be interested in the marginalized distribution q(T t ) = N (T t , m t , V t ) and the covariance between T t and T t-1 noted V t,t-1 .

We describe three different methods that look for gaussian approximations to the posterior: a forward-backward method adapted from [START_REF] Smith | Estimating a State-Space Model from Point Process Observations[END_REF], and two variants derived from the Gaussian Processes (GP) formalism, the Laplace approximation and Expectation Propagation.

The forward backward method is largely adapted from [START_REF] Brown | A statistical paradigm for neural spike train decoding applied to position prediction from ensemble firing patterns of rat hippocampal place cells[END_REF][START_REF] Smith | Estimating a State-Space Model from Point Process Observations[END_REF]. The E-step consists of three successive algorithms: a nonlinear recursive filter to compute p(T t |A (1,.,t) , θ) , which is where the Gaussian approximation is introduced; a Fixed Interval Smoothing We can also define the E-step as a GP inference problem, where we again look for a gaussian approximation to the cultural trait posterior. While the Forward-Backward algorithm models only works on the distribution of traits on successive time steps, GP inference approximates the whole history of cultural traits as a multivariate normal distribution. We define a gaussian prior for T , p(T |θ T ) = N (µ, K), which is provided by the generative model (see below), and inference is run by looking at an approximate posterior p(T |θ T , A) ≈ q(T ) = N (m, V ) [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF]. One possible concern with the forward-backward method is that errors may accumulate across time steps while we run gaussian approximations separately at each time step. By looking at an gaussian approximation directly over the entire dataset, we may significantly reduce this error. This is however at the expense of higher computations demands: while the forward-backward algorithm runs linearly in time (i.e. O(n)), the different GP algorithm will in general run cubically in the number of observations, and thus of time (i.e. O(n 3 ). This constraint can be alleviated in two directions: first, instead of using all of the time points in the dataset, we can restrict to time point with corresponding observations (artifacts). Second we may use sparse approximations to GP inference, which provide much more scalable solutions [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF]. We have not however explored this second possibility in this work.

From the auto-regressive model for cultural trait represented by equation 9, it follows that the prior over T will be normally distributed. The prior mean and covariances can be evaluated iteratively:

µ t = M µ t-1 + J t (22)
with µ 0 = x 0 1.

The prior covariance for T t follows:

K t,t = σ 2 dt t-1 v=0 (M T M ) v + σ 2 0 (M T M ) t (23)
The covariance between traits at two different time steps t and t + u (0 ≤ u ≤ n -t) is given by:

K t,t+u = M u-1 K t,t = M u-1 σ 2 dt t-1 v=0 (M T M ) v + σ 2 0 I(M T M ) t (24)
If the connectivity matrix G is symmetric, then M will be as well and

(M T M ) t = M 2t .
Once the mean and covariance for T are built, we select the sub-vector m obs and sub-matrix K obs corresponding to data points T obs with at least one observation (artifact). Then the GP inference corresponds to finding an approximate posterior N (T obs ; m obs , V obs ) for T obs , given a prior p(T obs ) = N (T obs ; µ obs , K obs ) and a series of observations {A tr } that each depends only on the value of the trait at one data point. The GP framework is extremely flexible as to the type of observations: they can be binary (GP classification), continuous (GP regression), or as for GLMs be parametrized by any exponential-family distribution, such as Poisson distribution for count data. The GPML toolbox for Matlab [START_REF] Rasmussen | Gaussian Processes for Machine Learning (GPML) Toolbox Carl[END_REF] allows to run inference for all these cases using a variety of different algorithms that each represent one case of gaussian approximation to the posterior. Here we use two different algorithms, Laplace approximation and Expectation Propagation (EP). Laplace approximation consists of identifying the mode m of the posterior distribution, i.e. the maximum a posteriori (MAP) value for the traits arg max T (p(T |θ T )p(A|T )). The March 4, 2021 19/28 posterior covariance matrix is given by the Hessian of the log-posterior function evaluated at m. EP is another approximation which usually provides more accurate approximation than Laplace for binary observations [START_REF] Nickisch | Approximations for Binary Gaussian Process Classification[END_REF]. See [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF] for a detailed presentation of these algorithms.

The GP inference algorithm outputs the posterior for traits corresponding to observed data points q(T obs ) = N (m obs , V obs ). To recover the full posterior for all points in the data set, we used the prediction equations:

m t = µ t + K t,obs K -1 obs (m obs -µ obs ) (25)
V t = K t,t -K T t,obs V -1 obs K t,obs (26)
V t-1,t = K t-1,t -K T t-1,obs V -1 obs K t,obs (27)
where K t,obs is the submatrix from K representing the prior covariance between traits at time point t and traits at all observed data points.

Once the gaussian approximation have been computed, using any of the methods presented above, we will interested in computing second moments of T t , and of T t and T t-1 jointly. We have directly:

W t = <T t T T t > q(Tt) = V t + m t m T t ( 28
) W t-1,t = <T t T T t-1 > q(Tt,Tt-1) = V t-1,t + m t-1 m T t (29)
B.2 M-step for continuous traits

In the M-step we maximize Q T i (θ; θ old ), Q T e (θ T e ; θ old ) and Q A (θ A ; θ old ) with respect to parameters θ.

Q T i (θ; θ old ) = q(T 0 ) log N (T 0 ; x 0 , σ 2 0 I) = r < -(x -x 0 ) 2 /(2σ 2 0) -log σ 2 0 /2> x∼q0r + const (30)
(note that actually the sum is not restricted to time step 0 but also applies to all regions created de novo and their time of creation)

<(x -x 0 ) 2 > x∼q0r = <x 2 > x∼q0r -2x 0 <x 2 > x∼q0r + x 2 0 = W 0 (r, r) -2x 0 µ 0 (r) + x 2 0 (31)
Plugging equation 31 into equation 30 and maximizing over x 0 and σ 2 0 , we obtain:

x 0 = r µ 0 (r)/R σ 2 0 = tr(W 0 )/K -x 2 0 (32)
tr is the trace operator: tr(M ) = r M (r, r). The corresponding maximum value is

Q T i (θ; θ old ) = R(1 -log σ 2 0 )/2.
We now seek to maximize Q T e (θ T e ; θ old ). From equation 9, we have:

log p(T t |T t -1, θ T e ) = log N (T t -M T t-1 -J t , σ 2 dtI) = -1 2σ 2 dt ||T t -M T t-1 -J t || 2 -K 2 log(2πσ 2 dt) (33)
Plugging this into equation 11, we obtain:

   Q T e (θ T e ; θ old ) = -Kn 2 log(2πσ 2 dt) -1 2σ 2 dt Z Z = n t=1 < T T t T t + T T t-1 M T M T t-1 + J T t J t -2T T t M T t-1 -2T t J T t + 2J T t M T t-1 > q (34) March 4, 2021 20/28
The sum in Z is defined in terms of first and second moments of the distributions q(T t ) and q(T t , T t-1 )

Z = t tr(W t ) + s( t W t-1 • M T M ) t J T t J t -2s(M • t W t-1,t ) -2 t µ t J t + 2s(M T • t µ t-1 J t ) (35)
• is the Hadamard product, i.e. (M • N ) ij = M ij N ij , and s refers to the sum of all elements in a matrix (s(M ) = ij M ij ). The range of subscript t in this equation and the followings is from 1 to n.

By substiting the definitions of M and J t (equation 9) into equation 35, we find that Z is a bilinear function of the subset of parameters ν = (γ, λ, ξ, ρ), that is:

                           Z = dt 2 ν T Cν -2dt d T ν + const, with C =     Kn 1 T E t 1 T G µ t-1 -1 T µ t-1 1 T E t E T t E t s( G • µ t-1 E T t ) -µ T t-1 E t 1 T G µ t-1 s( G • µ t-1 E T t ) β 2 -β 1 -1 T µ t-1 -µ T t-1 E t -β 1 β 0     d =     1 T (µ n -µ 0 ) (µ t -µ t-1 ) T E t s( G • W t-1,t ) -β 1 β 0 -tr( W t-1,t )     (36) We have defined χ = t W t-1 , β 0 = tr(χ), β 1 = s(χ • G) and β 2 = s(χ • G2
). The subindex t has been omitted from the sums. When some of the parameters in ν are fixed (i.e. ξ if we look at models without cultural diffusion), we simply remove rows and columns corresponding to the fixed parameters in C and d.

Maximizing Q T e (θ T e ; θ old ) w.r.t. ν is equivalent to minimizing Z, which is achieved at

ν = 1 dt C -1 d (37)
In order to find the new value of σ 2 , we look for the maximum of Q T e (θ T e ; θ old ) w.r.t. σ 2 using equation 34, yielding:

σ 2 = Z Kndt ( 38
)
where Z is computed using the values of M and J t with the new values of ρ, ξ, λ and γ. In summary we can use equations 32,37, 38 to update the values of parameters θ T . In case the artifact production model also contains some parameters θ A , we will seek to maximize:

Q A (θ A ) = t,r N (T tr ; m t (r), V t (r, r)) log p(A tr |T tr , θ A )dT tr (39)
In the case of count observations, where the observations are generated from Poisson distribution where the rate is defined as λ(r, t) = exp(λ 0 + λ 1 T t,r ), we can use Newton's method to find best fitting values for parameters (λ 0 , λ 1 ) as described in [START_REF] Smith | Estimating a State-Space Model from Point Process Observations[END_REF]. For the general case the integrals and their derivatives w.r.t parameters θ A have no analytical forms, but various approximations schemes can be applied, i.e. using Gaussian-Hermite quadrature [START_REF] Liu | A note on Gauss-Hermite quadrature[END_REF]. Such approximations can be used to find new set of parameters θ A that maximizes Q A . March 4, 2021 21/28

B.3 Hessian

Confidence intervals and marginalized evidence both rely on computing the Hessian of the log-likelihood, which is related to the Hessian of Q(θ) through equation 14. Using equations 30, 31,34 and 36, we find that the Hessian of Q T has the following form:

H Q T = ∇∇Q T (θ T ) =     -dtC/σ 2 0 0 0 0 -Kn 2σ 4 0 0 0 0 -K/σ 2 0 0 0 0 0 -K 2σ 4 0     (40)

B.4 Forward backward method

These equations are directly derived when adapting the formalism of [START_REF] Brown | A statistical paradigm for neural spike train decoding applied to position prediction from ensemble firing patterns of rat hippocampal place cells[END_REF][START_REF] Smith | Estimating a State-Space Model from Point Process Observations[END_REF] to the generative model of continouous ETC.

B.4.1 Nonlinear recursive filter

We define the approximations p(T t |A) = N (x t|t , W t|t ) and

p(T t+1 |A (1,.,t) ) = N (x t+1|t , W t+1|t )
, where the distributions of traits are conditioned on all artifacts up to time step t. Mean and covariance matrices are initiated at the first time step by x 0|0 = x 0 and W 0|0 = σ 2 0 I. From equation 9, we have the following prediction rule:

x t|t-1 = M x t-1|t-1 + J t (41)
W t|t-1 = M W t-1|t-1 M + σ 2 dtI (42)
In the special case where a new region r is initiated at t with parent region s (i.e.

split or a political state in two), the value of T t is inherited from the parent region, i.e.

x t|t (r) = x t|t (s) and W t|t (r, r) = W t|t (r, s) = W t|t (s, r) = W t|t (s, s). For creations de novo, we have x t|t (r) = x 0 and W t|t (r, r) = σ 2 0 .

In the case of continuous artifacts generated from gaussian noise with variance 2 , i.e. gaussian observations A tr ∼ N (T tr , 2 )), the posterior distribution after integrating new observation A t is itself gaussian, and we can calculate the new distribution p(T t |A 1..t ) = N (T t ; x t|t , W t|t ) by following the classical Kalman filter:

x t|t = x t|t-1 + K t (A t -x t|t-1 ) (43)
W t|t = (I -K t )W t|t-1 (44)
where we have used the Kalman gain matrix

K t = W t|t-1 (W t|t-1 + 2 I) -1 .
If observations are not gaussian, the step is not exact and we look for Gaussian approximations p(T t |A 1..t ) = N (T t ; x t|t , W t|t ). If artifacts take binary values, the distribution set by: log

p(T t |H t ) ∝ (T t -x t|t-1 ) T W -1 t-1|t-1 (T t -x t|t-1 )+ r (A + tr log σ(T tr )+A - tr log σ(-T tr )) (45
) We recognize the mode x t|t as being the solution of a logistic regression problem with observations A t , identity design matrix I and gaussian prior over weights p(x) = N (x; x t|t-1 , W t|t-1 ) (this can also be viewed as a Gaussian Process classification problem with prior mean x t|t-1 and prior covariance W t|t-1 ). The mode can simply be found solving this logistic regression problem, e.g using the IRLS March 4, 2021 22/28 algorithm. W t|t is the inverse of the Hessian of log p(T t |H t ) evaluated at the mode, i.e. it is the inverse of the observed Fisher information matrix.

In general, if observation process (artifacts) are derived from an generalized family distribution, i.e. Poisson counts, we can find the mode x t|t and W t|t by solving a GLM with prior mean gaussian prior over weights p(x) = N (x; x t|t-1 , W t|t-1 )). The solution is equivalent to that found by the recursive method proposed by [START_REF] Smith | Estimating a State-Space Model from Point Process Observations[END_REF].

B.4.2 Fixed Interval Smoothing

Defining B t = W t|t M W -1 t+1|t [START_REF] Söderkvist | An algorithm for Kalman filtering and smoothing with diagonal input covariance matrices[END_REF], we have the backward recursive equations:

m t = x t|t + B t (m t+1 -x t+1|t ) (46)
V t = W t|t + B t (V t+1 -W t+1|t )B T t (47) V t,t+1 = B t V t+1 (48)
(this last equation from [START_REF] Jong | Covariances for Smoothed Estimates in State-space Models[END_REF])

B.5 Relative influence of vertical transmission, horizontal transmission and ecological factors

We show here how the relative influence of vertical transmission I v , horizontal transmission I h and ecology I e onto the evolution of the cultural trait can be estimated from a dataset. Once maximum likelihood ETC parameters have been identified, we estimate the proportion of the variance of T that can be explained by each factor:

   I v = V v /SV I h = V h /SV I e = V e /SV (49)
where SV = V v + V h + V e + V n is the sum of variance, V v , V h , V e and V n are the variance terms of T due to vertical transmission, horizontal transmission, ecology and noise respectively. The different terms can be computed analytically by looking how the cultural trait in one region T (r, t), in different regions T (s, t) and the ecological factor E(r, t) determine cultural trait in the same region the cultural trait after one generation T (r, t + T g ), where T g is the duration of a generation (typically ∼ 25 years). For the case of continuous traits, we decompose the dynamics of T from equation 9 into

T = T v + T h + T e + T n
, where T h = Th -T v , and the dynamics of the different components follow the different equations:

       dTv dt = -ρT v + γ1 d Th dt = -M Th + γ1 dTe dt = -ρT v + E(t) dTn dt = -M T n + η(t) (50)
with initial conditions T v (t 0 ) = Th (t 0 ) = T (t 0 ) and T e (t 0 ) = T n (t 0 ) = 0. The solution to first (deterministic) differential equations are:

   T v (t) = T (t 0 ) + e -ρ(t-t0) (γ1 -T (t 0 )) Th (t) = T (t 0 ) + e -M (t-t0) (γ1 -T (t 0 )) T e (t) = (I -e -M (t-t0) )E(t 0 ) (51)
Here we assumed that the ecological factor is nearly constant within the time range of one generation. We use these equations to compute the value of T v , T h and T e after March 4, 2021 23/28

one generation (at t = t 0 + T g ) for all starting points t 0 in the dataset. Variance terms V v , V h and V e are obtained by computing the variance of the respective terms over regions and time points. Finally, we recognized that T n is a multivariate Orstein-Uhlenbeck process, whose variance is Var(T n (t)) = σ 2 M -1 (I -e -M (t-t0) ).

Appendix C Simulations C.1 Binary cultural trait

We simulated an ETC model with R = 6 regions over a period of 1000 times steps. The duration of each time step was drawn from a uniform distribution between 0.1 and 0.2 arbitrary units of time. Region 1 was created from region 3 at time step 501. Region 3 was created de novo at time step 301. All other regions were created de novo at time step 1. Region 3 was destroyed at time step 600, while all other regions were maintained for the rest of the simulation. The graphical connectivity was generated as a Erdös-Rényi random graph, i.e. the connection between each pair of regions was drawn with a fixed probabilty p = 0.3 independently of other connections. The ecological factor E was generated from independent stochastic process in each region, obtained by filtering a gaussian white noise process η rt : τ E E r,t + E r,t-1 = 5.65

τ E u=1 η r,t-u .
Parameters of the ETC were b up = b down = .1 (transition probabilities), λ = .03 (evoked parameter), ξ = 0.03 (diffusion parameter), x 0 = 0 (initial state of T ).

We used binary artifacts with value 0 or 1. The number of artifacts per data point was drawn from a Poisson distribution with rate 0.4. For each artifact, the probability that the artifact takes value 1 was determined by the value of the trait at the same time and region: p(A t,r = 1|T t,r = 1) = C up = .6, and p(A t,r = 1|T t,r = 0) = C down = 0.

After times series for ecological factors, cultural traits and cultural artifacts were synthesized, we used the algorithm described in section ?? to estimate ETC parameters (x 0 , b up , b down , λ, ξ, C up , C down ) from the observed data (E, A).

C.2 Continuous cultural trait

We simulated an ETC model with R = 5 regions over a period of 400 times steps, with fixed time step duration dt = 0.05 (for example representing a total span of 2000 years with time steps of 50 years). Region 2 was created from region 3 at time step 31.

Region 1 was created de novo at time step 51. All other regions were created de novo at time step 1. Region 3 was destroyed at time step 300, while all other regions were maintained for the rest of the simulation. The graphical connectivity was generated as a Erdös-Rényi random graph, i.e. the connection between each pair of regions was drawn with a fixed probabilty p = 0.91 independently of other connections. The ecological factor E was generated from an independent stochastic process in each region, similarly to binary cultural trait simulations. Parameters of the ETC were ρ = .08 (cultural lability), λ = 1 (evoked parameter), ξ = 0.1 (diffusion parameter), γ = .3 (trait default value), σ 2 = .2 (noise variance), x 0 = .6 and σ 2 0 = .5 (mean and variance of initial state of T ).

We used count variable for artifacts. The number of artifacts per data point was drawn from a Poisson distribution with rate 0.5. The value associated with each each artifact was sampled from a Poisson distribution with rate exp(T t,r ).

After times series for ecological factors, cultural traits and cultural artifacts were synthesized, we used the algorithm described in section B to estimate ETC parameters (x 0 , σ 2 0 , ρ, λ, ξ, γ, σ 2 ) from the observed data (E, A), using the three different approximation techniques for the E-step: forward-backward method, Laplace March 4, 2021 24/28 approximation and Expectation-Propagation. The analysis was repeated for 100 different simulations of the model.

C.3 Two stage analysis

For the analysis comparing one-stage vs two-stage approaches, we simulated an ETC model with R = 5 regions over a period of 500 times steps, with fixed time step duration dt = 0.05. Region 2 was created from region 3 at time step 31. All other regions were created de novo at time step 1, and all regions were maintained until the end of the simulation. The graphical connectivity was generated as a Erdös-Rényi random graph with a fixed probability of connection p = 0.91. The ecological factor E was generated as in previous analyses, and we used the same parameters as for the continuous analysis above.

We used a continuous variable for artifacts. The number of artifacts per data point was drawn from a Poisson distribution whose rate changed over time: the rate was set to .1 for the first 400 time steps (scarcity of artifacts period) and jumped to 1 for the remaining 100 time steps (abundance of artifacts period). The value associated with each each artifact was sampled from a normal distribution with mean T t,r and variance 2 = 1, so that the average of artifacts at any time point provides a direct unbiased estimate of the associated cultural trait.

After times series for ecological factors, cultural traits and cultural artifacts were synthesized, we compared the estimated provided by standard ETC method (with exact E step, since the observations are gaussian) with a two-stage procedure. In the two-step procedure, cultural traits are first estimated directly from cultural artifacts: estimated cultural trait T (r, t) is taken as the average of the value of all artifacts produced in region r in a window of 10 time steps centered on t. Then we estimate the ETC parameters by maximizing the likelihood over latent evolution p( T |θ T ) and the likelihood over cultural production p(A| T ; θ A ). This is indeed equivalent to maximizing the log of the joint-likelihood log p(A; T |θ) = log p(A| T , θ A ) + log p( T |θ T ), which is equivalent to the M step of the E algorithm where the variance of posterior cultural traits is null. The analysis was repeated for 100 different simulations of the model.

C.4 Synchronous data analysis

We simulated data from an ETC model with continuous traits over 100 regions over 400 time steps (dt = 0.02). Regions were created according to a random phylogenetic tree, which was generated as follow. We started the tree from a single region at time step 1 (created de novo), at coordinates (0, 0) on a map. At each iteration, we randomly selected one region out of the alive regions. This region had a random number between 2 and 4 offspring regions, after a random duration following its creation. This region was extinguished at the moment it generated offspring regions. Ecological variables in each region were drawn from independent binary markov processes with a probability of 0.5dt of changing value at each time step. The location for each new region was taken at a distance Inferring individual-level processes from population-level patterns in cultural evolution. from the normal distribution both in the x and y direction.

Fig 2 .

 2 Fig 2. Fitting a binary ETC model on synthetic data. A. True vs. inferred cultural traits across time in all regions (bottom to top, from region 1 to region 6). Full lines represent the true binary cultural trait T t,r , while dotted line represent the posterior probability of the cultural trait given by estimated parameters p(T t,r | θ), where parameters are estimated using the EM procedure from observed data (E, A). B. Value of parameters estimated with the EM procedure (gray bars) against true value of the parameters (red ticks). C Model comparison between full and reduced models. Value of ∆BIC = BIC(reduced) -BIC(full), for two different reduced models: evoked-only model and diffusion-only model.

28 Fig 3 .

 283 Fig 3. Trait variables estimated from continuous ETC model, using forward-backward, Laplace or Expectation-Propagation method in E-step. True trait value is represented by thick line, inferred trait posterior is represented as posterior mean (dotted line) and standard deviation (shaded area). Each colour represent values in a different region.

Fig 4 .

 4 Fig 4. A-G Estimated ETC parameters for each type of approximate method in E-step (Laplace approximation, EP or forward-backward method). Dotted line represent the true value. H: mean squared error on cultural trait. Bars: average value over datasets; error bars: standard deviation across datasets.

Fig 5 .

 5 Fig 5. Model selection for simulated datasets using different metrics (AIC, BIC, CVLL or boostrapping p-values) and approximate inference methods for E-step (forward-backward: blue bars; Expectation-Propagation: red bars; Laplace method: green bars). Each bar represents the proportion of simulated datasets where the corresponding model (full, no evoked model, no diffusion model or control model) was selected. CVLL was computed using 5-fold cross-validation.

Fig 6 . 28 Fig 7 .

 6287 Fig 6. Comparison of one-stage and two-stage estimation procedures. A. Left: Estimates of cultural trait using standard ETC model. Legend as in Figure 3. Vertical lines marks the transition from low to high cultural production period. Right: Direct estimates of cultural trait by averaging cultural artifacts in same region and within a certain window of time. B. Parameter estimates using one-stage (blue bars) vs two-stage (red-bars) method. Vertical dashed line indicate true value of the parameters. Error bars: standard deviation of estimates across datasets. C. Average mean-squared-error on cultural traits using one-stage vs. two-stage method (log-scale).

  a State-Space Covariance algorithm. The details of these algorithms are provided in Appendix B.4.

  Non-diagonal elements in the connectivity matrix G scaled exponentially with the euclidian distance between each regions d rs , i.e. G rs = e -drs . Parameters of the model were the following: ρ = 0.5, λ = 2, ξ = 1, σ 2 = 0.002, γ = -1, x 0 = -0. After simulating the data, we fitted the ETC model using synchronous data, i.e. the value of the cultural trait observed at the last time steps for all alive regions. We compared two variants of the ETC model: either the full model, or the model without horizontal transmission. March 4, 2021 25/28

Acknowledgments

Supporting information Appendix A Fitting binary ETC model A.1 E-step for binary traits

In the case of binary traits, exact inference can be performed using the standard forward-backward algorithm, irrespective of the model of artifact production. The algorithm yields the log-likelihood for current parameters θ. In the E-step, we compute the probability mass functions for the trait value and transitions given current values of the parameters [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF][START_REF] Rabiner | Dune Patterns, Namib Desert[END_REF], i.e. γ t (ε) = p(T t =ε|A; θ) = p(T t,1 =ε 1 , .., T t,K =ε K |A; θ) and ζ t (ε, ς) = p(T t =ε, T t+1 =ς|A; θ), defined over trait vectors ε and ς in {0, 1} K . The marginal posterior probabilities for the cultural trait γ r,t = p(T t,r |A; θ) and marginal posterior probabilities of transition ζ r,t = p(T t,r , T r,t+1 |A; θ) are then computed simply by marginalizing:

This research was supported by the Spanish Ministry of Economy and Competitiveness (RYC-2017-23231 to A.H.) and the Agence Nationale pour la Recherche (ANR-17-EURE-0017 Frontcog, ANR-10-IDEX-0001-02 PSL* to N.B.).