Last updated: 2026-01-06
Checks: 7 0
Knit directory: fiveMinuteStats/analysis/
This reproducible R Markdown analysis was created with workflowr (version 1.7.1). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(12345) was run prior to running the
code in the R Markdown file. Setting a seed ensures that any results
that rely on randomness, e.g. subsampling or permutations, are
reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
The results in this page were generated with repository version 2d6d41e. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
Note that you need to be careful to ensure that all relevant files for
the analysis have been committed to Git prior to generating the results
(you can use wflow_publish or
wflow_git_commit). workflowr only checks the R Markdown
file, but you know if there are other scripts or data files that it
depends on. Below is the status of the Git repository when the results
were generated:
Unstaged changes:
Modified: Makefile
Staged changes:
Modified: analysis/likelihood_ratio_simple_continuous_data.Rmd
Modified: analysis/likelihood_ratio_simple_models.Rmd
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
These are the previous versions of the repository in which changes were
made to the R Markdown (analysis/bayes_multiclass.Rmd) and
HTML (docs/bayes_multiclass.html) files. If you’ve
configured a remote Git repository (see ?wflow_git_remote),
click on the hyperlinks in the table below to view the files as they
were in that past version.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| html | 2d6d41e | Peter Carbonetto | 2026-01-06 | Ran wflow_publish("analysis/bayes_multiclass.Rmd"). |
| Rmd | a361e11 | Peter Carbonetto | 2026-01-06 | A bunch of small edits to the bayes_multiclass vignette. |
| Rmd | b705a3f | Peter Carbonetto | 2025-12-29 | Created pdf version of bayes_multiclass vignette. |
| html | b705a3f | Peter Carbonetto | 2025-12-29 | Created pdf version of bayes_multiclass vignette. |
| Rmd | 4b35829 | Peter Carbonetto | 2025-12-29 | Created pdf version of bayes_multiclass vignette. |
| html | 4b35829 | Peter Carbonetto | 2025-12-29 | Created pdf version of bayes_multiclass vignette. |
| html | 5f62ee6 | Matthew Stephens | 2019-03-31 | Build site. |
| Rmd | 0cd28bd | Matthew Stephens | 2019-03-31 | workflowr::wflow_publish(all = TRUE) |
| html | 34bcc51 | John Blischak | 2017-03-06 | Build site. |
| Rmd | 5fbc8b5 | John Blischak | 2017-03-06 | Update workflowr project with wflow_update (version 0.4.0). |
| Rmd | 391ba3c | John Blischak | 2017-03-06 | Remove front and end matter of non-standard templates. |
| html | 8e61683 | Marcus Davy | 2017-03-03 | rendered html using wflow_build(all=TRUE) |
| html | 5d0fa13 | Marcus Davy | 2017-03-02 | wflow_build() rendered html files |
| Rmd | d674141 | Marcus Davy | 2017-02-26 | typos, refs |
| html | c3b365a | John Blischak | 2017-01-02 | Build site. |
| Rmd | 67a8575 | John Blischak | 2017-01-02 | Use external chunk to set knitr chunk options. |
| Rmd | 5ec12c7 | John Blischak | 2017-01-02 | Use session-info chunk. |
| Rmd | 8bcd7f9 | mfrasco | 2016-04-04 | changed references to four groups to five groups |
| Rmd | 9156677 | stephens999 | 2016-01-19 | add multiclass example |
See here for a PDF version of this vignette.
Be familiar with the computing the posterior probability on classes for 2 classes.
Suppose we have an observation \(X\), which we believe to have come from one of \(K\) models, \(M_1 \dots, M_K\). Suppose we can compute the likelihood for any models. This document lays out the computation of the posterior probability that the model came from model \(K\), and emphasizes that the result depends only on the likelihood ratios. This is a straightforward extension of the 2-class calculations.
In a previous example, we considered the question of whether a tusk came from one of two classes—a savanna elephant or a forest elephant—based on its DNA. In practice, we might be interested in finer-scale distinctions than this. For example, forest elephants from West Africa actually differ genetically from those in Central Africa. And Savanna elephants from North Africa differ from those in the East and South. (Actually, elephant genetics within each subspecies varies roughly continuously across the continent, and any division into discrete groups can be seen as an approximation.)
So what if now we have allele frequencies for “North Savanna”, “South Savanna”, “East Savanna”, “West Forest” and “Central Forest” groups? How do we decide which group a tusk likely came from? Now we have five models, but the calculation is the same for \(K\) models, so we might as well do it for \(K\) models. Here is the general outline.
Suppose we are presented with a series of observations \(x_1, \dots, x_n\), each of which are generated from a model \(M_k\), for some \(k \in {1, \dots, K}\). Let \(Z_i \in {1, \dots, K}\) indicate which model the \(i\)th observation came from. Bayes Theorem says that \[ \Pr(Z_i = k \mid x_i) = \frac{\Pr(x_i \mid Z_i = k) \, \Pr(Z_i=k)}{\Pr(x_i)}. \]
Recall that, by the law of total probability, \[ \Pr(x_i) = \sum_{k'=1}^K \Pr(x_i, Z_i=k') = \sum_{k'=1}^K \Pr(x_i \mid Z_i=k') \, \Pr(Z_i=k'). \] Also note that \(\Pr(x_i \mid Z_i = k)\) is the “likelihood” for model \(k\) given data \(x_i\), which we write as \(L_{ik}\). And we use \(\pi_k\) as shorthand for \(\Pr(Z_i=k)\). Putting this together gives \[ \Pr(Z_i = k \mid x_i) = \frac{L_{ik} \pi_k}{\sum_{k'=1}^K L_{ik'} \pi_{k'}}. \] Note that the denominator \(\Pr(x_i)=\sum_{k'=1}^K L_{ik'} \pi_{k'}\) is the same for all \(k\). So an equivalent way of laying out this calculation is to write \[ \Pr(Z_i = k \mid x_i) \propto L_{ik} \pi_k \] and to note that the constant of proportionality is determined by the fact that probabilities must sum to 1. This way of applying Bayes theorem is very common and convenient in practice, so you should get used to it. In words, this formula can be said \[ \text{posterior} \propto \text{likelihood x prior}. \]
Here is an example of the calculation in practice. The five rows of
the matrix ref_freqs represent the allele frequencies in
five groups: North Savanna, South Savanna, East Savanna, West Forest and
Central Forest. The calculation presented here assumes that the
population of tusks we are looking at is equally drawn from all four
groups, so \(\pi =
(0.2,0.2,0.2,0.2,0.2)\), but it would of course be easy to change
to any other value of \(\pi\).
x <- c(1,0,1,0,0,1)
ref_freqs <- rbind(
c(0.39,0.14,0.22,0.12,0.03,0.38),
c(0.41,0.10,0.18,0.12,0.02,0.28),
c(0.40,0.11,0.22,0.11,0.01,0.3),
c(0.75,0.25,0.11,0.18,0.25,0.25),
c(0.85,0.15,0.11,0.16,0.21,0.26)
)Function to compute the posterior probabilities given \(L\) and \(\pi\):
normalize <- function (x)
x/sum(x)
posterior_prob <- function (L_vec, pi_vec)
normalize(L_vec * pi_vec)Now another function to compute the likelihood:
L <- function (f, x)
prod(f^x*(1-f)^(1-x))Now use these functions to compute the likelihoods and posterior probabilities:
L_vec <- apply(ref_freqs,1,L,x = x)
L_vec
# [1] 0.023934466 0.016038570 0.020702326 0.009513281 0.013712299
posterior_prob(L_vec,c(0.2,0.2,0.2,0.2,0.2))
# [1] 0.2852705 0.1911608 0.2467472 0.1133871 0.1634344Remember that when comparing two models, only the likelihood ratios matter, not the actual likelihoods. In fact the same is true when comparing \(K\) models, as we can see by examining the calculation above. Specifically, imagine multiplying all the likelihoods by some positive constant \(c\), and notice that this would not change the final answer, because of the normalization step.
Notice that, just as with the 2-model case, the calculation involves weighing the relative support from the data for each model (from the likelihood function) against the “prior” plausibility of each model (from the vector \(\pi\)).
In practice we might not know \(\pi\). And although in such a case it might seem natural to assume that all the values of \(\pi\) are equal, one has to be careful to note that this is still an assumption, and such assumptions may have unforeseen implications. For example, in this case, this assumption implies that 60% of the tusks are from savanna elephants and 40% from forest elephants, not 50-50 (because three of our five groups are savanna). The difference between 60-40 and 50-50 is probably not a big deal in most applications, but imagine that we had 20 different savanna groups and 2 forest groups. Would we still be happy to assume that every group was equally common (and so savanna tusks are 10 times as common as forest tusks)? The answer would depend on the context, but quite possibly not.
sessionInfo()
# R version 4.3.3 (2024-02-29)
# Platform: aarch64-apple-darwin20 (64-bit)
# Running under: macOS 15.7.1
#
# Matrix products: default
# BLAS: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib
# LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
#
# locale:
# [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#
# time zone: America/Chicago
# tzcode source: internal
#
# attached base packages:
# [1] stats graphics grDevices utils datasets methods base
#
# loaded via a namespace (and not attached):
# [1] vctrs_0.6.5 cli_3.6.5 knitr_1.50 rlang_1.1.6
# [5] xfun_0.52 stringi_1.8.7 promises_1.3.3 jsonlite_2.0.0
# [9] workflowr_1.7.1 glue_1.8.0 rprojroot_2.0.4 git2r_0.33.0
# [13] htmltools_0.5.8.1 httpuv_1.6.14 sass_0.4.10 rmarkdown_2.29
# [17] evaluate_1.0.4 jquerylib_0.1.4 tibble_3.3.0 fastmap_1.2.0
# [21] yaml_2.3.10 lifecycle_1.0.4 whisker_0.4.1 stringr_1.5.1
# [25] compiler_4.3.3 fs_1.6.6 Rcpp_1.1.0 pkgconfig_2.0.3
# [29] later_1.4.2 digest_0.6.37 R6_2.6.1 pillar_1.11.0
# [33] magrittr_2.0.3 bslib_0.9.0 tools_4.3.3 cachem_1.1.0