Simulating model for several interaction types

[achauffou]

At the end of last week, I have prepared models that can handle multiple interactions types. I have performed simulations on the model all_binom_03 (its Stan code can be found there). The model can be summarized as follows:

logit(p[site, sp1, sp2]) = alpha + beta[site] + gamma[sp1] + gamma[sp2] + lambda[site] * SS[site, sp1, sp2]
gamma[sp] = zgamma[sp] * sigma_gamma[sp_fun_group[sp]]
lambda[site] = zlambda[site] * sigma_lambda[site_int_type[site]] + lambda_bar[site_int_type[site]]
...

Each species belongs to a functional group and each site to an interaction type, and the sigma_gamma and sigma_lambda depend on the functional group, resp. interaction type

I have performed three different simulations:

1 interaction type, 20 sites, 52 species, species at all sites of their interaction type, removing no data point: param_post_true_alpha.pdf param_post_true_beta.pdf param_post_true_gamma.pdf param_post_true_lambda.pdf param_post_true_lambda_bar.pdf param_post_true_sigma_beta.pdf param_post_true_sigma_gamma.pdf param_post_true_sigma_lambda.pdf These results look really nice, nothing alarming.

3 interaction types, 40 sites, 60 species, species at all sites of their interaction type, removing no data point: param_post_true_alpha.pdf param_post_true_beta.pdf param_post_true_gamma.pdf param_post_true_lambda.pdf param_post_true_lambda_bar.pdf param_post_true_sigma_beta.pdf param_post_true_sigma_gamma.pdf param_post_true_sigma_lambda.pdf It does not look completely off, but some parameters are not estimated that well (sigma_gamma, lambda_bar).

3 interaction types, 90 sites, 900 species, species are not present at all site, removing empty columns/rows: param_post_true_alpha.pdf param_post_true_beta.pdf param_post_true_gamma.pdf param_post_true_lambda.pdf param_post_true_lambda_bar.pdf param_post_true_sigma_beta.pdf param_post_true_sigma_gamma.pdf param_post_true_sigma_lambda.pdf The model has very small confidence intervals for alpha, beta, sigmas, but some parameters are completely off (alpha, sigma_beta, sigma_gamma). There must be something wrong, the models is very confident in estimations that are quite off.

I am currently investigating the issue as well as starting to think about how I will write the manuscript. Let me know if by any chance you have an idea of what is going wrong in here.

[bernibra]

Two things (so far):

• Alphas should also be interaction-dependent. That is, alpha should be "alpha[sp_fun_group[sp]]" (using your notation). This way you allow different interaction types to be more or less likely to exist than others.
• I thought we decided the lambdas to be constant? That made a lot of sense to me, and I think it simplifies the model quite a bit. If you end up making them constant though, just like with alpha, it should be interaction-dependent. That is, "lambda[sp_fun_group[sp]]" (using your notation).

[achauffou]

You are right I forgot to make the alphas interaction-dependent. Regarding your second point, I have prepared all models that we did for the case of pollination interactions. But it is true I should maybe only use the ones with constant interaction-dependent lambdas since we are not going to talk about the others when speaking of multiple-interactions models. However, interestingly for pollination models the WAIC of the model with pooled lambdas is slightly lower than with constant lambda (but the different in WAIC is very very small in comparison to the other models). I will keep you updated.

[bernibra]

The reason I feel that the single lambdas might be better is simply that it will be easier to interpret. That said, I'm not too surprised that the WAIC is lower. You should keep all, and see where the discussion takes us.

[achauffou]

I fully agree it definitely makes things simpler

[achauffou]

Just a quick update on my current struggles with these simulations...

• The simulation with a single lambda works fine, has good diagnostics and seems to deliver results that are appropriate. That is encouraging.
• The simulation with pooled lambdas seems to deliver appropriate results, but has 1768/8000 transitions that reach the max tree depth. I am not sure yet whether I should increase the max allowed tree depth to try and get rid of those, or if it indicates a more complex issue.
• Both simulations with/without pooled lambda completely fail when I try to generate data that have missing data points (to account for the fact that not all species are seen everywhere). There are 6000 (resp. 8000) out of 8000 transitions reaching max tree depth and large R hat values for all parameters. I am still investigating the issue.

Oddly enough, the two problems that I have here did not happen yesterday when I was working with less simulated data. So for now I suspect that I have made an error somewhere in the code.

[bernibra]

Do you want to meet to talk about this? I am not concerned about hitting the max tree depth, it can happen often and it does not necessarily indicate problems with the model but only difficulty to sample. Moreover, I think we should focus on single lambdas for your master thesis because, the more I think about it, I do not necessarily see why we would expect lambdas to vary per site (since we already have an intercept that varies and captures sampling differences across sites). Third problem you mention is a bit harder for me to comment on. I don't remember how we were generating the missing data, but it sounds like you might have some species that don't have data to work with and the model is struggling.

Anyway, let me know if you want to meet!

[achauffou]

It could be a good idea to meet, but maybe I should first do a few things to ensure that the issue is real (rather than meet only to realize later that I made a mistake in declaring a variable or else).

Regarding the pooled lambdas, I agree that it might not be the biggest concern now (and will likely drop the pooled lambdas models to focus on single lambda ones). Initially, I would expect that the influence of climatic suitability could vary depending on the local conditions and, technically, I think the model could capture such effect due to the difference between slope/intercept. But I guess the amount of data we have and uncertainty on suitability does not reasonably allow to infer such a small effect considering how difficult to capture it is. So I agree that it is more reasonable to use only the single-lambda models.

Regarding the missing data simulations, I worked on ensuring that all species/sites are in the generated data yesterday, but I think it is not foolproof yet. I will keep investigating the issue this afternoon, and if I do not find out where it comes from, I will let you know so we can meet and talk about it.

[bernibra]

Sounds good. I am currently on holiday, so I am happy to wait until next week. I was just offering in case you felt stuck in this, but it's good that you want to double-check the data-generating process before the meeting.

[achauffou]

Great, then we can wait and see next week if we need to meet. I can work on other things in parallel if I end up stuck anyway (starting to draft the manuscript, collecting data on invasiveness). Enjoy your holiday!

[achauffou]

Just so you know, I have found an error that caused the issue yesterday. The simulation outcomes now looks good, and I think I will run the model on the real data soon.

[bernibra]

Fantastic news! Great work. Looking forward to seeing the results for the real data.