Phylogenetic meta-analysis: model optimization and WAIC calculation

[Marina Golivets]

Dear All,

I am conducting a global meta-analysis on pairwise plant species interactions in order to understand whether invasive species are generally better competitors than native species. 

Each data point in the dataset represents an outcome of competition between two different species, with each species being classified according to its biogeographic origin as either native or nonnative (invasive). The two species origins give four combinations of interacting species pairs, all of which are represented in the dataset: native-native, native-nonnative, nonnative-native, and nonnative-nonnative. The overall goal of the analysis is to conclude on whether the intensity of competition between those four groups significantly differs.

In my model, I've been trying to incorporate multiple potential sources of nonindependence among effect sizes (Hedges' d), the main three of which are:

1) the study effect (effect sizes coming from the same experiment/research group are most probably more similar that those from different experiments);

2) sampling nonindependence, which arises when the same control was used to estimate multiple effect sizes;

3) phylogenetic relatedness among competitor species and phylogenetic relatedness among target species (closely related competitors might have similar effects on target species and, similarly, closely related target species may respond in a similar way to competition).

I got the model run, and it returns plausible results. I am not sure, however, how to specify pointwise log-likelihood in Stan that I would need to perform model evaluation and comparison. 

I've attached the data, Stan model, and all the R code needed to perform the analysis. I would enormously appreciate any suggestions regarding the model itself (mistakes, ways to optimize it, etc.) and, even more importantly, any help on the model evaluation approach that would work best in this case.

Kind regards,

Marina

[Aki Vehtari]

Hi,

WAIC is a method to approximate predictive performance, and thus you need to define which prediction task is interesting for you. Instead of WAIC it is easier to think about that by considering cross-validation (leave-one-out cross-validation being asymptotically equivalent with WAIC, but having generally better finite sample performance). If you consider cross-validation, you need to think at which level you want to make the predictions. Could you tell whether you would be interested in predictions for new (or left-out)

 - single observations (for same studies, sampling groups etc.)

 - studies

 - sampling groups

 - biogeogr origin groups

or

 - species?

(If I did understood correctly your model, then I found these different ways to group the data).

Easiest and most robust way to do the predictive performance estimate, would then be to use k-fold-CV by leaving out one or several "things" at the desired level (observation, study, sampling group, biogeogr origin group, species).

A computationally faster way is to use leave-one-thing-out with PSIS-LOO approximation.

You could also easily compute log likelihood terms required by PSIS-LOO for sampling groups by adding following to the generated quantities block

for (i in 1:K){

    log_lik[i] = multi_normal_lpdf(segment(d, pos, L[i]) | segment(theta, pos, L[i]), block (V_d[i], 1, 1, L[i], L[i]));

but then it depends how large L[i]'s are whether importance sampling approach works for this.

Luckily loo package will give you warnings if this doesn't produce a well behaving approximation (note that is PSIS-LOO fails then WAIC fails even more badly).

Thus to continue you need to first think what you want to predict.

Aki

[Michael Betancourt]

+1

[Marina Golivets]

Aki,

Thank you so much for your help! 

I am primarily interested in predictions at the level of biogeographic origin. 

Marina

[Aki Vehtari]

Hmm... I tried to look more of your code, but now I'm confused

Looking at these lines in your code

---

   int <lower=1> N; // num obs

   int <lower=1> G; // num biogeogr origin groups

   matrix[N, G] X; // design matrix

...   

   vector[G] beta; // group predictors

...

   for (i in 1:N) theta[i] = mu_group[i] + eff_study[jj[i]] + eff_comp_phylo[mm_comp[i]] + eff_targ_phylo[mm_targ[i]] + eff_comp_species[mm_comp[i]] + eff_targ_species[mm_targ[i]] + sigma_obs * eta_obs[i];

---

each theta can depend on all biogeogr origin groups, and I can't see what you want to predict at biogeogr origin level. Probably needs more time to understand this...

Aki

[Marina Golivets]

X is a matrix of ones and zeroes, i.e. if observation i belonged to the biogeographic origin group 2, the row X[i, ] would be 0 1 0 0. 

I used matrix notation because I was planning on adding covariates in the model and I thought it would be easier to modify the model later. It is, however,  not necessary; instead, I could provide indexing for a group in the data block:

...

int gg[N]; // biogeographic origin group index

...

mu_group ~ mu + sigma_mu*eta;

   for (i in 1:N) theta[i] = mu_group[gg[i]] + eff_study[jj[i]] + eff_comp_phylo[mm_comp[i]] + eff_targ_phylo[mm_targ[i]] + eff_comp_species[mm_comp[i]] + eff_targ_species[mm_targ[i]] + sigma_obs * eta_obs[i];

I am trying to predict overall mean effect for each group.

Marina

[Bob Carpenter]

Multiplying by 0 or 1 is much less efficient in Stan than
doing a select. That is, if

int<lower=0, upper=1> cond;
real a;

rather than writing

cond * a

it's better to write

cond ? a : 0

or even better yet, to just skip all the calculations where the
zero would propagate. So if you're adding in an intercept, you'd
just have something like:

if (cond) predictor = predictor + a;

The reason is that all the arithmetic operations get stored in the
expression graph we travers to propagate the chain rule for derivatives.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Aki Vehtari]

If each K term in here

--- 

    for (i in 1:K){

        segment(d, pos, L[i]) ~ multi_normal(segment(theta, pos, L[i]), block (V_d[i], 1, 1, L[i], L[i]));

---

depends only on one biogeographic origin group, then when you compute log_liks in generated quantities

---

for (i in 1:K){

    log_lik[i] = multi_normal_lpdf(segment(d, pos, L[i]) | segment(theta, pos, L[i]), block (V_d[i], 1, 1, L[i], L[i]));

---

you can sum log_lik's for each biogeographic origin group to get G log_lik's to be given to loo package.

And as I said previously PSIS-LOO might be able to handle this, but if not it will tell you and then you need to make G-fold-CV or k-fold-CV by yourself.

Aki

[Marina Golivets]

Thank you, Aki!

Unfortunately, effect sizes within sampling groups might be coming from different biogeographic origin groups. 

Here is a short explanation why it is the case: 

Each sampling group is a collection of effect sizes that were calculated using a single control. In this study, a control was the biomass of a target species when grown in monoculture, and then a treatment was the biomass of the same target species when grown in a mixture with a competitor. Sampling dependencies arose from the experiments where species-specific effects of multiple competitor species (among which there could be both native and non-native species) were tested on the same target species. 

Marina