Access summary statistics (e.g. fit quantiles) data

[mschubert]

Thank you for this package! I started using it a couple of days ago, and I'm impressed how simple and intuitive it makes MCMC model fitting.

If I want to get the data behind the fit quantiles shown in the plot, how would you recommend to do that?

I'm exploring reimplementing parts of the plot.mcpfit function and extracting the data from the ggplot object.

Does it maybe make sense in the medium term to split these summary statistics into separate functions for a user to be able to access too?

[lindeloev]

Thanks a lot for the encouraging words!

All the MCMC data is in fit$samples which is an mcmc.list object. This is a widely used format so a lot of packages could be put to work on that and mcp should not duplicate it. In particular, tidybayes (perhaps in concert with ggplot2), coda, and bayesplot would be useful for all sorts of plotting and diagnostics. You can see a few lines on how to get started here: https://github.com/lindeloev/mcp#do-much-more

But I'm eager to hear if you have a use case that would benefit from mcp having some "intermediate" functions?

And yes, splitting out some of the functionality from plot.mcpfit would definitely be good for the internal structure of mcp. This is on the to-do list for making a proper predict function that works for out-of-sample data.

[lindeloev]

BTW, are you aware that summary.mcpfit returns a data.frame with quantiles? :-)

result = summary(fit)
print(result$lower[1])  # lower quantile of the first parameter

[mschubert]

Thank you, this is very helpful! A predict method would indeed solve this issue for me (as I'm interested in the prediced y values, not the parameters in this case).

But I'm eager to hear if you have a use case that would benefit from mcp having some "intermediate" functions?

I'm trying to extract the best fitted value and its confidence intervals on the domain of the original data.

I have read up a bit on tidybayes and used your plotting function to come up with fitted quantiles like this (maybe this helps someone):

pars_pop = paste0(fit$pars$population, collapse = "|")
samples = tidybayes::spread_draws(fit$mcmc_post, !!rlang::sym(pars_pop), regex=TRUE) %>%
    tidybayes::sample_draws(1000) %>%
    tidyr::expand_grid(x=seq(min(fit$data$x), max(fit$data$x), length.out=30)) %>%
    dplyr::mutate(low = rlang::exec(fit$simulate, !!!., type="fitted",
                                    quantile=0.05, add_attr=FALSE)) %>%
                  high = rlang::exec(fit$simulate, !!!., type="fitted",
                                    quantile=0.95, add_attr=FALSE)) %>%
    dplyr::group_by(x) %>%
    dplyr::summarize(low = mean(low),
                     high = mean(high))

This does, however, seem to fit only the median and not the high/low quantiles:

ggplot(samples, aes(x=x)) +
    geom_point(aes(y=low), color="green") +
    geom_line(aes(y=high), color="red")

whereas if I use your plotting function I get:

plot(fit, q_fit=c(0.05,0.95))

I guess simulate does not support fitted quantiles, but I have to predict and calculate them myself?

[lindeloev]

fit$simulate(..., quantile = XXX) is a work-in-progress that I haven't finished yet so quantile has no effect as of now and that's why you end up with identical line for the two quantiles. Indeed, it will eventually do what you expect it to. Sorry about that! Once it's finished, I plan on writing a vignette about the fit$simulate() function because it will be a workhorse underlying much of mcp post-processing , including the upcoming predict() function.

In the meantime, try replacing the last part of your first code chunk with:

dplyr::mutate(pred= rlang::exec(fit$simulate, !!!., type="fitted", add_attr=FALSE)) %>%
dplyr::group_by(x) %>%
dplyr::summarize(
    low = quantile(pred, 0.05),
    high = quantile(pred, 0.95))

(untested, so let me know if this solves it)

[mschubert]

That works perfectly, thank you!

[lindeloev]

mcp 0.3 has now been published on CRAN. You can use predict(fit) and fitted(fit) to get quantiles. See more in this article and in the documentation for these functions.