Update `coev_get_delta_theta()` and `coev_plot_delta_theta()` functions

[ScottClaessens]

It would be useful to update these functions to use the new coev_calculate_theta() function provided by @erik-ringen.

[ScottClaessens]

I have updated these functions in a4a8e79 and 514abf6. Both functions now rely on coev_calculate_theta() under the hood.

It's worth noting that in the updated function versions, $\Delta\theta_{z}$ is calculated using posterior medians of the median trait values, rather than looping over the posterior. See e.g. this line. This is because coev_calculate_theta() takes a fitted model and returns a full posterior, so it's not possible to loop over the posterior samples in the same way as before.

@erik-ringen do you think this is okay? Or is there important posterior uncertainty in the median trait values that we are missing with this approach?

[ScottClaessens]

Closing this for now, we can come to this point if need be.

[erik-ringen]

I have updated these functions in a4a8e79 and 514abf6. Both functions now rely on coev_calculate_theta() under the hood.

It's worth noting that in the updated function versions, Δ θ z is calculated using posterior medians of the median trait values, rather than looping over the posterior. See e.g. this line. This is because coev_calculate_theta() takes a fitted model and returns a full posterior, so it's not possible to loop over the posterior samples in the same way as before.

@erik-ringen do you think this is okay? Or is there important posterior uncertainty in the median trait values that we are missing with this approach?

Re-opening this because I think it is important to integrate uncertainty in the medians and MADs. My intuition is this will stabilize the estimation of the delta-theta-z values, which I notice now often have really extreme tail values.

If we extract the posterior into a rethinking-style list with the now built-in function extract_samples(), then it should be easy, unless I am missing something.

post <- extract_samples(object) eta <- post$eta[ ,1:object$stan_data$N_tips, ] med_eta <- apply(eta, c(1,3), median) mad_eta <- apply(eta, c(1,3), mad)