vmc is a visualization grammar designed to make it easy to generate informative visualizations for model checking. The vmc grammar assumes a basic workflow for creating model checks. First, the model predictions or model features need to be extracted as a distribution of data from a model object, i.e. data tidying. Then the user must specify an uncertainty representation to describe the selected distribution(s). They must also specify the presentation of the observed data. The user can choose among multiple comparative layouts to structure their comparison between observed data and model predictions.
Installation
You can install the development version of vmc from GitHub with:
# install.packages("devtools")
devtools::install_github("MUCollective/vmc")Usage
Model
First, we fit a beta regression model to predict the penguins’ bill ratio (i.e., bill depth / bill length) based on penguins’ flipper length.
penguins <- penguins %>%
drop_na(sex) %>%
mutate(is_gentoo = species == "Gentoo") %>%
mutate(bill_ratio = bill_depth_mm / bill_length_mm)
model_beta <- brm(
bf(bill_ratio ~ flipper_length_mm,
phi ~ flipper_length_mm),
family = Beta(),
init = "0",
data = penguins,
prior = c(prior(normal(0, 1), class = "b"),
prior(exponential(1), class = "b", dpar = "phi", lb = 0)),
file = "models/model.rds"
)The results…
model_beta
#> Family: beta
#> Links: mu = logit; phi = log
#> Formula: bill_ratio ~ flipper_length_mm
#> phi ~ flipper_length_mm
#> Data: penguins (Number of observations: 333)
#> Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#> total post-warmup draws = 4000
#>
#> Population-Level Effects:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
#> Intercept 3.06 0.16 2.75 3.37 1.00 4082
#> phi_Intercept -1.35 1.30 -3.93 1.22 1.00 1844
#> flipper_length_mm -0.02 0.00 -0.02 -0.02 1.00 4210
#> phi_flipper_length_mm 0.03 0.01 0.02 0.04 1.00 1854
#> Tail_ESS
#> Intercept 2918
#> phi_Intercept 1479
#> flipper_length_mm 2768
#> phi_flipper_length_mm 1625
#>
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).Check on the shape of distrition
You can start model checking by check the shape of the distribution of predicted bill ratio and observed bill ratio. Using vmc, you can create a posterior predictive check by only one mcplot() (and also using coord_flip() to flip the response variable to the x-axis).
model_beta %>%
mcplot() +
mc_gglayer(coord_flip())
#> Warning: No shared levels found between `names(values)` of the manual scale and the
#> data's fill values.
You can not only check on the shape of posterior predictive distribution, but also check on other posterior distributions generated by model_beta using the mc_draw component. For example, there are two parameter mu and phi besides the response variable in the beta regression model we are using. Here, we are check on the posterior distribution of mu.
model_beta %>%
mcplot() +
mc_draw("mu") +
mc_observation_transformation(mean, group = vars(flipper_length_mm)) +
mc_gglayer(coord_flip())
#> Warning: No shared levels found between `names(values)` of the manual scale and the
#> data's fill values.
Check on the marginal effects
You may find the misalignment between the shape of predictive distribution and observed distribution. You can scrutinize the misalignment further by checking on the marginal effects of the predictor using mc_condition_on.
model_beta %>%
mcplot() +
mc_draw(ndraws = 1) +
mc_condition_on(x = vars(flipper_length_mm))
#> Warning: No shared levels found between `names(values)` of the manual scale and the
#> data's fill values.
To reveal the uncertainty information inside of the model predictions, you can use either HOPs or other uncertainty representations.
vmc supports HOPs as a method to group samples (see vignette("uncertainty-representation")).
model_beta %>%
mcplot() +
mc_draw(ndraws = 50) +
mc_model_point(group_sample = "hops") +
mc_condition_on(x = vars(flipper_length_mm))
#> Warning: No shared levels found between `names(values)` of the manual scale and the
#> data's fill values.
#> `nframes` and `fps` adjusted to match transitionYou may find that some of the points of model predictions and observed data are overlapped together, which may cause visual clutter. To fix this problem, you can change to another comparative layout using mc_layout_*.
model_beta %>%
mcplot() +
mc_draw(ndraws = 50) +
mc_model_point(group_sample = "hops") +
mc_condition_on(x = vars(flipper_length_mm)) +
mc_layout_juxtaposition()
#> Warning: No shared levels found between `names(values)` of the manual scale and the
#> data's fill values.
#> No shared levels found between `names(values)` of the manual scale and the
#> data's fill values.
#> No shared levels found between `names(values)` of the manual scale and the
#> data's fill values.
#> No shared levels found between `names(values)` of the manual scale and the
#> data's fill values.
#> `nframes` and `fps` adjusted to match transition
#> # A tibble: 50 × 7
#> format width height colorspace matte filesize density
#> <chr> <int> <int> <chr> <lgl> <int> <chr>
#> 1 gif 800 400 sRGB FALSE 0 28x28
#> 2 gif 800 400 sRGB FALSE 0 28x28
#> 3 gif 800 400 sRGB FALSE 0 28x28
#> 4 gif 800 400 sRGB FALSE 0 28x28
#> 5 gif 800 400 sRGB FALSE 0 28x28
#> 6 gif 800 400 sRGB FALSE 0 28x28
#> 7 gif 800 400 sRGB FALSE 0 28x28
#> 8 gif 800 400 sRGB FALSE 0 28x28
#> 9 gif 800 400 sRGB FALSE 0 28x28
#> 10 gif 800 400 sRGB FALSE 0 28x28
#> # ℹ 40 more rows
Other uncertainty representations are supported by the visual representation components (mc_model_* and mc_obs_*).
model_beta %>%
mcplot() +
mc_draw(ndraws = 1000) +
mc_model_lineribbon() +
mc_condition_on(x = vars(flipper_length_mm))
Facet to find the source of multimodal
You may find there seem to be some clusters among the points of observed bill ratio. To find the source of this multimodal, you can change the conditional variable to facet the model checking on rows or columns using mc_condition_on.
model_beta %>%
mcplot(observation = penguins) +
mc_draw(ndraws = 1000) +
mc_model_lineribbon() +
mc_condition_on(x = vars(flipper_length_mm),
row = vars(sex))
model_beta %>%
mcplot(observation = penguins) +
mc_draw(ndraws = 1000) +
mc_model_lineribbon() +
mc_condition_on(x = vars(flipper_length_mm),
row = vars(island))
model_beta %>%
mcplot(observation = penguins) +
mc_draw(ndraws = 1000) +
mc_model_lineribbon() +
mc_condition_on(x = vars(flipper_length_mm),
row = vars(species))
New model
When faceting on the species of penguins, you can see the points of observed bill ratio are clustered into three groups. This finding can inform an update on the model to include species as an independent variable.
model_species <- brm(
bf(bill_ratio ~ flipper_length_mm * species,
phi ~ flipper_length_mm * species),
family = Beta(),
init = "0",
data = penguins,
prior = c(prior(normal(0, 0.5), class = "b"),
prior(normal(0, 0.15), class = "b", dpar = "phi")),
cores = 4,
iter = 24000,
warmup = 14000,
file = "models/model_species.rds"
)The result…
model_species
#> Family: beta
#> Links: mu = logit; phi = log
#> Formula: bill_ratio ~ flipper_length_mm * species
#> phi ~ flipper_length_mm * species
#> Data: penguins (Number of observations: 333)
#> Draws: 4 chains, each with iter = 24000; warmup = 14000; thin = 1;
#> total post-warmup draws = 40000
#>
#> Population-Level Effects:
#> Estimate Est.Error l-95% CI u-95% CI
#> Intercept -0.40 0.26 -0.91 0.12
#> phi_Intercept 7.76 2.29 3.25 12.24
#> flipper_length_mm 0.00 0.00 -0.00 0.00
#> speciesChinstrap -0.13 0.32 -0.76 0.50
#> speciesGentoo -0.45 0.31 -1.05 0.15
#> flipper_length_mm:speciesChinstrap -0.00 0.00 -0.00 0.00
#> flipper_length_mm:speciesGentoo -0.00 0.00 -0.00 0.00
#> phi_flipper_length_mm -0.01 0.01 -0.04 0.01
#> phi_speciesChinstrap -0.01 0.15 -0.31 0.28
#> phi_speciesGentoo 0.00 0.15 -0.29 0.30
#> phi_flipper_length_mm:speciesChinstrap 0.01 0.00 0.00 0.01
#> phi_flipper_length_mm:speciesGentoo 0.01 0.00 0.00 0.01
#> Rhat Bulk_ESS Tail_ESS
#> Intercept 1.00 21184 25569
#> phi_Intercept 1.00 30353 28724
#> flipper_length_mm 1.00 21186 25744
#> speciesChinstrap 1.00 21067 24050
#> speciesGentoo 1.00 21687 24792
#> flipper_length_mm:speciesChinstrap 1.00 20879 23602
#> flipper_length_mm:speciesGentoo 1.00 20853 23819
#> phi_flipper_length_mm 1.00 30039 28653
#> phi_speciesChinstrap 1.00 35917 28882
#> phi_speciesGentoo 1.00 36834 27541
#> phi_flipper_length_mm:speciesChinstrap 1.00 38256 30833
#> phi_flipper_length_mm:speciesGentoo 1.00 27289 27179
#>
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).You can check on this new model on the marginal effects of flipper length and species in one plot by using mc_condition_on.
model_species %>%
mcplot() +
mc_draw(ndraws = 1000) +
mc_model_lineribbon() +
mc_condition_on(x = vars(flipper_length_mm),
color = vars(species))
By this model checking, the new model is validated to capture the multimodal in the observed data. Then you can check more on the posterior distribution of the model conditioned on species using vmc.
model_species %>%
mcplot() +
mc_model_slab(side = "left") +
mc_obs_dots() +
mc_condition_on(x = vars(species))
Check on model assumptions
To check on some particular assumptions of the model, you can pick different mark types of the visual representation. For example, to check the calibration of the model, you can choose to use intervals to represent the model distribution.
model_species %>%
mcplot() +
mc_model_interval() +
mc_obs_interval() +
mc_condition_on(x = vars(species)) +
mc_layout_nested()
#> Warning: No shared levels found between `names(values)` of the manual scale and the
#> data's fill values.
vmc also supports to check on other assumptions, like heteroskedasticity and the normality of residuals.
model_species %>%
mcplot() +
mc_draw(ndraws = 1) +
mc_condition_on(x = vars(flipper_length_mm),
color = vars(species)) +
mc_layout_encoding("residual") +
mc_gglayer(geom_hline(yintercept = 0))
model_species %>%
mcplot() +
mc_draw(ndraws = 1) +
mc_condition_on(color = vars(species)) +
mc_layout_encoding("qq") +
mc_gglayer(geom_abline())