Mathsy Coding

Using Bayesian Regression to Quantify Sexual Dimorphism in Dinosaurs

The full text of a paper I presented a poster on at the Alberta Palaeontological Society's Paleo 24 conference, examining the use of Bayesian regression to quantify sexual dimorphism in dinosaurs.

Introduction

Sexual dimorphism is the difference in size, colour, or presence of a feature between different sexes in a species. For instance, generally human males are larger than females, male mallard ducks are colourfully ornamented while females are drab, and male deer have antlers while females do not. Sexual dimorphism is very common in extant species, but detection of sexual dimorphism in extinct taxa, especially vertebrates, has been more challenging. Essentially, there are two difficulties in determining sexual dimorphism in extinct taxa: the ‘sexual’ part and the ‘dimorphism’ part. Determining the sex of the individual fossil, in the absence of rare features (such as medullary bone or bacula) is extremely challenging (Saitta et al. 2020). Similarly, determining the mass or size of a feature from a fossil (such as those of dinosaurs) can be difficult (Brusatte 2012 p. 126). When trying to determine the level of sexual dimorphism, these problems compound each other.

There are several analytical methods which are commonly used. However, each of these has their weaknesses. The use of tt tests (for differences in the mean) or Hartigans’ dip test (for unimodality) (Hartigan and Hartigan, 1985) can suffer from low statistical power (Mallon 2017) and in the case of tt tests a further difficulty is the prolonged growth period of most dinosaurs (Hone and Mallon 2017). As a result, positive results for dimorphism in dinosaurs are rare.

However, a proposed method in (Saitta et al. 2020) attempts to sidestep these difficulties by explicitly constructing growth curves for the different sexes. Herein is proposed a modification to that method in which Bayesian regression is used to recover a distribution of possible growth curves for each sex and subsequently applied to calculate the level of dimorphism.

Methods

Bayesian regression is quite different from classical regression. To estimate a parameter, there needs to be a prior representing the information that you have before seeing the data. You can then calculate the likelihood of the data that you observed for the different parameter values. The optimal combination of the prior and the likelihood is given by Bayes’ Theorem and is termed the posterior. In general, the result of a Bayesian analysis will be a probability distribution, not a single number (McElreath 2020). An illustration of this method is shown in Figure 1.

To estimate human height, we start with what we know before seeing any data --- the prior. Then a sample of heights is collected and the likelihood of that sample for each of the parameter values is calculated. The combination of those, our final belief about the value of the parameter, is the posterior distribution. Figure 1

To estimate human height, we start with what we know before seeing any data - the prior. Then a sample of heights is collected and the likelihood of that sample for each of the parameter values is calculated. The combination of those, our final belief about the value of the parameter, is the posterior distribution.

The modification to the method in (Saitta et al. 2020) is presented below.

Saitta et al. 2020Modification
Gather data: age (or a correlate) and character of interest (e.g. length)
Fit growth curve to population as a whole
Predict sex by assigning all samples above the population curve to one sex and all below to the other.
Fit separate growth curves to each predicted sexCreate a posterior distribution for sex-specific parameters (Bayesian regression)
Calculate maximum likelihood level of dimorphismCalculate posterior distribution for dimorphism

In both methods, the amount of dimorphism is calculated as the difference in the asymptotes of the growth curves.

Method Validation

Validating the method began with a simulated population of American alligators (Alligator mississippiensis). The original von Bertalanffy growth curves were determined in Wilkinson and Rhodes (1997) and a simplified model using a constant standard deviation (rather than sex-specific standard deviation of the form αlog⁡(mean mass at age)+β\alpha \log(\text{mean mass at age}) + \beta) was used to generate different samples. The final model was

Lengthm/f∼N(μm/f,σm/f)μm=3.79∗(1−0.94e−0.0695∗age)μf=2.78∗(1−0.91e−0.0926∗age)σm/f=0.25\begin{aligned} \text{Length}_{ m / f } &\sim \mathcal{N}\left(\mu_{ m / f }, \sigma_{m/f}\right) \\ \mu_m &= 3.79 * \left(1 - 0.94e^{-0.0695 \ast \text{age}}\right) \\ \mu_f &= 2.78 * \left(1 - 0.91e^{-0.0926 \ast \text{age}}\right) \\ \sigma_{m/f} &= 0.25 \end{aligned}

For the von Bertalanffy equation

Length=L(1−Ae−K∗age)\text{Length} = L(1 - Ae^{-K \ast \text{age}})

the parameter LL is the asymptotic length, and so the actual level of sexual dimorphism is 3.79 m - 2.78 m = 1.01 m.

Simulated Alligator Populations & Method Figure 1: Simulated Alligator Populations & Method

A Bayesian model with wide priors centred around the population-level curve values was used to recover the values (R Core Team 2023; Stan Development Team 2023).

Lengthm/f,i∼N(μm/f,i,σm/f)μm/f,i=von Bertalanffy(agei,Lm/f,Am/f,Km/f)Lm/f∼N(population L,0.5)Am/f∼N(population A,0.025)Km/f∼N(population K,0.05)σm/f∼N(0.25,0.05)\begin{aligned} \text{Length}_{m/f,i} &\sim \mathcal{N}\left(\mu_{m/f, i}, \sigma_{m / f}\right) \\ \mu_{m/f, i} &= \text{von Bertalanffy}(\text{age}_i, L_{ m/f }, A_{ m/f }, K_{ m/f }) \\ L_{ m/f } &\sim \mathcal{N}\left(\text{population L}, 0.5\right) \\ A_{ m/f } &\sim \mathcal{N}\left(\text{population A}, 0.025\right) \\ K_{ m/f } &\sim \mathcal{N}\left(\text{population K}, 0.05\right) \\ \sigma_{m/f} &\sim \mathcal{N}\left(0.25, 0.05\right) \\ \end{aligned}

Since the priors used for the male and female parameters are identical (and initially have means equal to the population-level parameters), the results are slightly biased toward a finding of minimal or no sexual dimorphism.

Both when using only single sex (female) or when including both, the model was able to recover all of the population parameters with a high degree of accuracy (in all cases within the 95% credible interval for the posterior of each parameter).

The sample size also affects the ability of the method to recover the true values. The model was run with several different sample sizes and the posterior distribution for the level of dimorphism Lm−LfL_m - L_f is shown below.

How different sample sizes affect the estimation of the level of dimorphism Figure 2: How different sample sizes affect the estimation of the level of dimorphism

This method reliably generates mean levels which are close to the true value and the 95% credible intervals overlap it as well. The general pattern of approaching the value from below is due to the influence of the prior (centred at zero); as the sample size becomes larger they drag the posterior toward the correct value.

We can also ask about how the magnitude of dimorphism affects our ability to detect it. For this, the female parameters were kept at their original value and simulated male populations were created with an effect size EE representing the scaled difference between the actual male and female parameters. An effect size of E=0E = 0 means that the male and female populations are identical, E=1E = 1 is the natural level of dimorphism, and E=2E = 2 is twice the natural level of dimorphism. The model was then run and the error in the level of dimorphism was calculated.

Effect of sample size on estimated dimorphism level Figure 3: Effect of sample size on estimated dimorphism level

There is a tendency to overestimate dimorphism when it is small, but converges to the correct value once it reaches a large enough size. Note that a significant contributor to the error in dimorphism estimation is the inaccuracy of our naive sex predictor. By plotting the accuracy of the dimorphism estimation against the computed accuracy in sex prediction, we can see a clear relationship.

Relationship between the error in dimorphism estimation and the accuracy of sex prediction Figure 4: Relationship between the error in dimorphism estimation and the accuracy of sex prediction

Overall, the results of the simulation show that this method is reliable, with caveats that at low sample sizes or levels or dimorphism we need to be especially skeptical of the results.

Results

Once the method had been validated, the same method was used to calculate the posterior distribution for sexual dimorphism for three species of dinosaur: Maiasaura peeblesorum, Psittacosaurus lujiatunensis, and Tyrannosaurus rex. The data used were compiled in (Saitta et al. 2020) and were originally from Woodward et al. 2015 (Maiasaura peeblesorum), Erickson, Makovicky, Inouye, et al. 2015 (Psittacosaurus lujiatunensis), and Erickson, Makovicky, Currie, et al. 2004, Horner and Padian 2004, and Lee and Werning 2008 (Tyrannosaurus rex).

In each case, a logistic growth curve for mass against age with a constant standard deviation was used.

logistic(age,L,q,k)=L1+eq+k∗age\text{logistic}(\text{age}, L, q, k) = \frac{L}{1 + e^{q+k \ast \text{age}}}

Since LL is the asymptotic mass, the level of dimorphism is Ll−LsL_l - L_s, where here ll represents the larger sex and ss represents the smaller one. For each species, priors were chosen again centred around the population level parameters and with a wide enough standard deviation to cover all of the existing data.

For each species, the prior 95% quantiles for randomly generated populations, along with the actual data, are displayed, along with a selection of posterior growth curves.

Maiasaura peeblesorum Prior and Posterior Estimates
Psittacosaurus lujiatunensis Prior and Posterior Estimates
Tyrannosaurus rex Prior and Posterior Estimates
Figure 5: Prior and Posterior Estimates

The posterior curves are biologically plausible and in line with both the data and our understanding of physiology; none of the values are immediately out of the realm of possibility.

Of course, we can also compare the levels of dimorphism as the larger sex’s size as a percent of the smaller one (e.g. a dimorphism level of 50% would mean that the larger sex was 50% larger than the smaller).

Dimorphism Levels Figure 6: Dimorphism Levels

From this, we can see strong evidence for the Maiasaura peeblesorum and Psittacosaurus lujiatunensis dimorphism, while the Tyrannosaurus rex displays either no or a very small amount of dimorphism.

Discussion

As shown above, the modification of the method in (Saitta et al. 2020) using Bayesian regression to calculate the posterior distributions for sex-specific growth curves and the level of sexual dimorphism is a generally reliable and useful way to quantify dimorphism in extinct populations. While it shares some weakness with the original method, it also has several notable advantages. One is that by carrying uncertainty through all parts of modelling, the final result contains more information about the process (McElreath 2020). Another is that by incorporating prior knowledge, we can constrain our final result to be biologically plausible.

Despite these advantages, there remain many areas in which this method could be improved. The most critical, as noted in (Saitta et al. 2020), is sex determination. As shown by Figure 4, sex determination accuracy was a strong predictor of the error in the level of dimorphism. There are two complementary paths which should be taken towards alleviating this issue. The first is the incorporation of the uncertainty in sex determination into the posterior curves generated. At the moment we behave as though the sex that we predict is correct, with no uncertainty, when the reality is that this is a significant source of error. This would result in widening the posterior distributions for the levels of dimorphism, but would not address the underlying issue. For that, we should look at alternate means of predicting sex, such as cluster analysis involving traits other than the one that we are looking at.

Determining levels of sexual dimorphism when both sex and size determination are difficult is naturally a difficult problem, and the method here, while offering an incremental improvement over existing methods, is certainly no panacea. Further efforts, both on the statistical and palaeontological side of the equation, will be required in order to improve our understanding of sexual dimorphism within extinct taxa.

For all data and code used here, along with additional analyses, see https://github.com/EricRobertCampbell/aps-2024-poster-sexual-dimorphism.

References