vignettes/repeated.Rmd
repeated.RmdA common feature of comparative datasets is repeated observations. In many analyses, we only have one observation per taxon. But often there will be more than one observation for each taxon, such as when we have observed multiple individuals of the same species. In these cases, it can be useful to include all of these observations in the model and estimate the residual variation that is not due to the coevolutionary process.
When repeated observations are present, we add taxon-level varying intercepts \(\alpha\) to the model as follows:
\[ \begin{split} y_{i,j} &\sim f(\mu_{i,j}, \phi_j) \\ g(\mu_{i,j}) &= \boldsymbol{\eta}_{i,j} + \alpha_{i,j} \\ \alpha_{i,j} &\sim N(0, \Sigma) \end{split} \]
where \(y_{i,j}\) is the observed value for trait \(j\) from taxa \(i\), \(f()\) is a likelihood function with expected value \(\mu\) and possible additional parameters \(\phi\), \(g()\) is a link function converting the expected value to the latent continuous space of the evolutionary process, \(\boldsymbol{\eta}_{i,j}\) is the latent trait value of trait \(j\) for taxa \(i\), \(\alpha_{i,j}\) is the taxon-level varying intercept for trait \(j\) for taxa \(i\), and \(\Sigma\) is a \(j\) x \(j\) variance-covariance matrix for the varying intercepts. The varying effects represent within-species/population variance in the latent trait values.
We show this using an example dataset from de Villemereuil & Nakagawa (2014). Suppose we have measured two continuous variables (\(x\) and \(y\)) for 20 species, with five observations for each species.
head(repeated$data) species x y
1 sp_1 11.223724 107.41919
2 sp_1 9.805934 109.16403
3 sp_1 10.308423 91.88672
4 sp_1 8.355349 121.54341
5 sp_1 11.854510 105.31638
6 sp_2 4.314015 64.99859We can fit the dynamic coevolutionary model to this dataset.
fit <-
coev_fit(
data = repeated$data,
variables = list(
x = "normal",
y = "normal"
),
id = "species",
tree = repeated$phylogeny,
# additional arguments for cmdstanr
parallel_chains = 4,
iter_warmup = 2000,
iter_sampling = 2000,
refresh = 0,
seed = 1
)Running MCMC with 4 parallel chains...
Chain 3 finished in 272.3 seconds.
Chain 4 finished in 366.5 seconds.
Chain 2 finished in 367.6 seconds.
Chain 1 finished in 373.3 seconds.
All 4 chains finished successfully.
Mean chain execution time: 344.9 seconds.
Total execution time: 373.4 seconds.
summary(fit)Variables: x = normal
y = normal
Data: repeated$data (Number of observations: 100)
Phylogeny: repeated$phylogeny (Number of trees: 1)
Draws: 4 chains, each with iter = 2000; warmup = 2000; thin = 1
total post-warmup draws = 8000
Autoregressive selection effects:
Estimate Est.Error 2.5% 97.5% Rhat Bulk_ESS Tail_ESS
x -1.58 0.77 -3.17 -0.21 1.00 5111 3604
y -1.34 0.71 -2.85 -0.15 1.00 5223 3187
Cross selection effects:
Estimate Est.Error 2.5% 97.5% Rhat Bulk_ESS Tail_ESS
x ⟶ y 0.42 0.86 -1.27 2.10 1.00 7501 6161
y ⟶ x 0.04 0.85 -1.65 1.71 1.00 7967 6517
Drift parameters:
Estimate Est.Error 2.5% 97.5% Rhat Bulk_ESS Tail_ESS
sd(x) 2.93 0.32 2.37 3.60 1.00 5491 5809
sd(y) 2.58 0.29 2.05 3.21 1.00 5128 6351
cor(x,y) 0.87 0.06 0.71 0.96 1.00 3770 4970
Continuous time intercept parameters:
Estimate Est.Error 2.5% 97.5% Rhat Bulk_ESS Tail_ESS
x -0.06 0.78 -1.61 1.46 1.00 8432 5963
y 0.04 0.76 -1.48 1.53 1.00 9003 5867
Residual parameters:
Estimate Est.Error 2.5% 97.5% Rhat Bulk_ESS Tail_ESS
sd(x) 0.28 0.02 0.23 0.32 1.00 9704 5265
sd(y) 0.29 0.02 0.25 0.34 1.00 10692 5832
cor(x,y) -0.18 0.10 -0.38 0.02 1.00 8699 5553Warning: There were 30 divergent transitions after warmup.
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmupIn the model output, we can see that coev_fit() has
detected the presence of repeated observations and has consequently
modelled residual standard deviation and residual correlation parameters
for \(x\) and \(y\).
If we wanted to remove the additional varying effects structure when
repeated observations are present, we could set
estimate_residual = FALSE when fitting the model.