## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", eval = FALSE)

## -----------------------------------------------------------------------------
# library(bqmm)
# data(Orthodont, package = "nlme")
# 
# fit <- bqmm(distance ~ age + (1 | Subject),
#             data = Orthodont,
#             tau  = 0.5)          # the conditional median
# 
# summary(fit)

## -----------------------------------------------------------------------------
# fixef(fit)        # population-level coefficients at tau = 0.5
# ranef(fit)        # subject-specific deviations
# VarCorr(fit)      # random-effect SDs (and correlations, if any)
# coef(fit)         # fixed effects
# predict(fit)      # fitted conditional medians

## -----------------------------------------------------------------------------
# fit_q <- bqmm(distance ~ age + (1 | Subject),
#               data = Orthodont,
#               tau  = c(0.1, 0.25, 0.5, 0.75, 0.9))
# 
# coef(fit_q)   # a tau-by-coefficient matrix
# plot(fit_q)   # coefficient-versus-tau paths

## ----echo = FALSE, eval = TRUE, out.width = "70%", fig.cap = "Coefficient-versus-quantile path: the estimated effect of a predictor at each quantile, with uncertainty. A non-flat path is distributional information a mean model discards."----
knitr::include_graphics("figures/coef-path.png")

## -----------------------------------------------------------------------------
# predict(fit_q, noncrossing = "rearrange")

## -----------------------------------------------------------------------------
# my_prior <- bqmm_prior(
#   beta_sd     = 5,    # SD of the normal prior on fixed effects
#   sigma_scale = 1,    # half-normal scale for the ALD scale sigma
#   re_scale    = 2,    # half-normal scale for random-effect SDs
#   lkj         = 2     # LKJ shape (correlated REs only)
# )
# fit <- bqmm(distance ~ age + (1 | Subject), Orthodont,
#             tau = 0.5, prior = my_prior)

## -----------------------------------------------------------------------------
# vcov(fit, adjusted = TRUE)    # corrected (default)
# vcov(fit, adjusted = FALSE)   # naive posterior covariance
# confint(fit, adjusted = TRUE)
# summary(fit)                  # uses the adjusted intervals

## ----echo = FALSE, eval = TRUE, out.width = "70%", fig.cap = "Frequentist coverage of nominal-95% intervals across simulated designs: the naive posterior under-covers; the Yang--Wang--He--adjusted intervals are at or above nominal."----
knitr::include_graphics("figures/coverage.png")

## -----------------------------------------------------------------------------
# fit_c <- bqmm(distance ~ age + (1 + age | Subject),
#               data = Orthodont, tau = 0.5,
#               cov  = "unstructured")
# 
# VarCorr(fit_c)                       # SDs plus...
# attr(VarCorr(fit_c), "correlation")  # the posterior-median correlation matrix

## -----------------------------------------------------------------------------
# fit <- bqmm(distance ~ age + (1 | Subject), Orthodont, tau = 0.5,
#             chains = 4, iter = 4000,
#             control = list(adapt_delta = 0.99, max_treedepth = 12))

## -----------------------------------------------------------------------------
# library(posterior)
# summarise_draws(as_draws(fit))            # R-hat, ESS per parameter
# 
# library(bayesplot)
# mcmc_trace(as_draws(fit), regex_pars = "b_")

## -----------------------------------------------------------------------------
# yrep <- posterior_predict(fit)            # draws x observations
# bayesplot::ppc_dens_overlay(Orthodont$distance, yrep[1:50, ])

## ----echo = FALSE, eval = TRUE, out.width = "70%", fig.cap = "Posterior predictive check: the observed outcome density (dark) against draws from the fitted model (light). Systematic mismatch flags a misfit the quantile of interest may not capture."----
knitr::include_graphics("figures/ppcheck.png")