Allow multiPhylo object in tree argument?

[ScottClaessens]

Currently, the package only allows single phylo objects. Should we allow the user to enter a multiPhylo object? The model could iterate over all (or a reduced number of) phylogenetic trees and report the combined estimates.

Or, alternatively, we could just leave this to the user to do manually (by e.g. looping over the coev_fit() function themselves).

[erik-ringen]

Now that the package is further along, I have some thoughts about how to implement this. We would need to segment the tree each time, such that that are N_tree x N_seg, node_seq, parent, and ts. z_drift would also become [N_tree, N_seg - 1, J].

To marginalize over the sampled trees in the model block, store the likelihood from each tree like this:

for (i in 1:N_obs){
vector[N_tree] lp;
for (n in 1:N_tree) lp[n] = some_lpdf(y[i] | eta[n, tip_id[i]) + log(1/N_tree);
target += log_sum_exp(lp);
}

[ScottClaessens]

I've made some changes to the Stan code on the multi-phylo branch. I'm not sure yet exactly how to modify the generated quantities block, but I'll have another look at this soon.

[ScottClaessens]

Okay, I've finished implementing multiPhylo objects in the Stan code, data list, and all post-processing functions. You can find the changes in the multi-phylo branch.

I have tested the authority model in the README against the same model fitted with two identical phylogenies in a multiPhylo object. The results are basically identical to the README results, which is comforting since there should be no phylogenetic uncertainty in this case:

Variables: political_authority = ordered_logistic 
           religious_authority = ordered_logistic 
     Data: authority$data (Number of observations: 97)
    Draws: 4 chains, each with iter = 1000; warmup = 1000; thin = 1
           total post-warmup draws = 4000

Autoregressive selection effects:
                    Estimate Est.Error  2.5% 97.5% Rhat Bulk_ESS Tail_ESS
political_authority    -0.62      0.52 -1.89 -0.02 1.00     2308     1622
religious_authority    -0.68      0.53 -1.97 -0.03 1.00     2398     1873

Cross selection effects:
                                          Estimate Est.Error 2.5% 97.5% Rhat Bulk_ESS
political_authority ⟶ religious_authority     2.62      1.23 0.26  5.08 1.01     1093
religious_authority ⟶ political_authority     2.49      1.30 0.02  5.18 1.00      868
                                          Tail_ESS
political_authority ⟶ religious_authority     1718
religious_authority ⟶ political_authority     1655

Drift scale parameters:
                    Estimate Est.Error 2.5% 97.5% Rhat Bulk_ESS Tail_ESS
political_authority     1.49      0.80 0.12  3.19 1.00     1141     1170
religious_authority     1.31      0.78 0.08  2.99 1.00     1102     1119

Continuous time intercept parameters:
                    Estimate Est.Error  2.5% 97.5% Rhat Bulk_ESS Tail_ESS
political_authority     0.19      0.91 -1.57  1.93 1.00     4314     2443
religious_authority     0.21      0.93 -1.60  2.04 1.00     4297     2734

Ordinal cutpoint parameters:
                       Estimate Est.Error  2.5% 97.5% Rhat Bulk_ESS Tail_ESS
political_authority[1]    -1.21      0.87 -2.93  0.62 1.00     2246     1815
political_authority[2]    -0.47      0.85 -2.11  1.24 1.00     2332     1949
political_authority[3]     1.69      0.94  0.00  3.62 1.00     1546     1742
religious_authority[1]    -1.54      0.87 -3.25  0.18 1.00     2512     2800
religious_authority[2]    -0.87      0.85 -2.57  0.81 1.00     2585     2881
religious_authority[3]     1.59      0.96 -0.15  3.59 1.00     2125     2309